*Article* **GPR and Digital Survey for the Diagnosis and the 3D Representation of the Battle of Issus Mosaic from the House of the Faun, Pompeii (Naples, Italy)**

**Marilena Cozzolino 1,\*, Antonio De Simone 2, Vincenzo Gentile 1, Paolo Mauriello <sup>1</sup> and Amanda Piezzo <sup>3</sup>**


**Abstract:** The application of non-invasive geophysical techniques and digital surveys to explore cultural heritage is becoming a very important research field. The capability to detect inner and superficial changes in the inspected surfaces allows for imaging spatial inhomogeneity and material features and planning targeted conservation and restoration interventions. In this work, the results of a research project carried out on the famous Battle of Issus Mosaic, also known as the "Alexander Mosaic", are presented. It is a masterpiece of ancient art that was found in 1831 in the House of Faun, the most luxurious and spacious house in Pompeii. It is notable for its size (3.41 × 5.82 m), the quality of workmanship and the subject that represents the culminating phase of the battle between Alexander Magno's army and the Persian one of Darius. In 1916, it was moved inside the National Archaeological Museum of Naples, where the original horizontal location was changed with a vertical arrangement supported by an inner wooden structure, whose exact manufacture is unclear. Today, the mosaic is affected by important instability phenomena highlighted by the appearance of the significant detachment of tiles, superficial lesions and swelling of the surface. Given the important need to preserve it, a high-detail diagnostic study was realized through a digital survey and non-invasive geophysical surveys using ground-penetrating radar (GPR). The investigation was repeated after two years, in 2018 and 2020, with the aim of verifying the evolution of degradation. The work provided a high-resolution estimate of the state of the health of the mosaic and allowed for obtaining a three-dimensional reconstruction of the internal mosaic structure, including the formulation of hypotheses on the engineering supporting works of the twentieth century; this provides an essential tool for the imminent conservation project, which also implies restoring the original horizontal position.

**Keywords:** mosaic of alexander; GPR; digital survey; pre-conservation diagnosis

#### **1. Introduction**

In recent years, the use of non-invasive geophysical and geomatic techniques has assumed an increasingly important role in the field of cultural heritage, especially by supporting conservation and restoration projects. Geophysical methods are very useful for assessing the presence of underground structures in preventive archeology at different scales [1–6]; addressing conservation and stability issues of architectural monuments by inspecting soil foundations; assessing the mechanical properties of structural elements [7–13]; and exploring internal and superficial structures of precious and delicate targets, such as statues, wall paintings and mosaics [7,14–18]. Among these techniques, ground-penetrating radar (GPR) is widely used thanks to the miniaturization of the instrumentation, the high

**Citation:** Cozzolino, M.;

De Simone, A.; Gentile, V.; Mauriello, P.; Piezzo, A. GPR and Digital Survey for the Diagnosis and the 3D Representation of the Battle of Issus Mosaic from the House of the Faun, Pompeii (Naples, Italy). *Appl. Sci.* **2022**, *12*, 6965. https://doi.org/ 10.3390/app12146965

Academic Editor: Tung-Ching Su

Received: 7 June 2022 Accepted: 6 July 2022 Published: 9 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

investigation resolution and the minimal impact on the analyzed surfaces. For example, Masini et al. [7] presented GPR prospecting on three different constructive elements that are typical of historical buildings (a wall, a masonry pillar and a marble column), allowing for the characterization of the masonry, the detection of cracks and the imaging of metallic reinforcement bars. The obtained information was relevant for providing better knowledge of the history of the monuments and their current internal state to be evaluated for in any possible restoration intervention. Matias et al. [14] used GPR combined with seismic transmission tomography to provide important results in an investigation of columns and walls of a 14th-century UNESCO monument by giving information on the quality and spatial distribution of the materials used in the construction of the monument. Manataki et al. [15] applied different GPR systems with frequencies of 1600 MHz, 500 MHz and 250 MHz to study the mosaics of Delos island, which holds a significant body of ancient Greek art of the Hellenistic period that is on UNESCO's World Heritage List. The 1600 MHz system allowed for identifying the boundaries of the mosaic layers, as well as problematic areas, such as bulges and high levels of moisture that may cause deterioration.

For the same topics, geophysical methods can also be applied jointly with geomatic techniques that, thanks to the advances in data collection, processing and visualization, represent an important source of knowledge regarding diagnostics and documentation [19–26]. Cozzolino et al. [20] conducted a 3D metric survey through photogrammetry and groundpenetrating radar (GPR) tests applied to the study of the trapezophoros with two griffins attacking a doe of Ascoli Satriano, which is a masterpiece of ancient art that needs to be protected. The work provided information on both visible and hidden defects, such as numerous cracks that affect the sculpture. Arias et al. [21] showed the usefulness of a multidisciplinary approach to heritage documentation involving close-range photogrammetry and ground-penetrating radar techniques, as well as the development of finite-elementbased structural models. The study was focused on the documentation of a medieval bridge concerning the geometric shape, the building material, and the current damage and its causes. Danese et al. [25] showed a spatial-analysis-based protocol for the interpretation of data coming from different non-invasive tests to improve the extraction process of the pattern's decay. The case study was a frescoed wall of a Gymnasium in Pompeii, which was investigated with the following non-invasive techniques: structure-from-motion photogrammetry (SfM), ground-penetrating radar and multitemporal infrared thermography. This approach enabled the extraction of decay patterns to construct a 3D model that constituted the deformation map of the painting analysis methodology and to establish the first step for the restoration of an important multilevel characterization of a fresco that is useful for the protection and mitigation of its deterioration risk.

In this study, a high-detail diagnostic investigation was realized through a digital survey and non-invasive geophysical survey using ground-penetrating radar (GPR) on the Battle of Issus Mosaic, also known as the "Alexander Mosaic". It was found on 24 October 1831 in the exedra of the peristyle of the House of the Faun (Regio VI, 12, 2) in Pompeii, which was the ancient city destroyed by the eruption of Vesuvius in 79 AD (Figures 1 and 2) [27]. In 2018, the National Archaeological Museum of Naples (MANN), where the mosaic is exposed, set up a working group, mainly composed of personnel within the institute, which is collaborating with the Central Institute for Restoration, as well as the University of Molise and the Center for Research on Archaeometry and Conservation Science (University of Naples Federico II and University of Sannio) for surveys and diagnostic investigations. The study determining the state of conservation of the mosaic was aimed at supporting the restoration project and the change of its placement. To this end, a high-detail survey was realized by the University of Molise in 2018, and it was repeated in 2020 to verify the evolution of degradation in compliance with the following workflow:

(1) The creation of a high-resolution 3D model and orthophotos of the external surfaces through a photogrammetry digital survey in order to highlight the decay and even the type of decay, which is not perceivable by direct sight.


**Figure 1.** Location of Pompeii on the map of Italy (**a**), plan of the city with the position of the House of Faun (**b**) and the exedra of the peristyle where the mosaic was found (**c**).

**Figure 2.** A close-up picture of the Battle of Issus Mosaic.

#### **2. Test Site: The Battle of Issus Mosaic**

The Battle of Issus Mosaic represents the triumph of Alexander the Great over Darius III of Persia in 333 BC during the Battle of Issus in Turkey. The masterpiece was created at the end of the second century BC in the opus vermiculatum technique with about 1,800,000 millimetric tesserae. In the left portion, Alexander the Great is depicted with his horse Bucephalus. Medusa is represented on the cuirass with her wavy hair. On the right, there is Dario on a chariot as he tries to launch an assault with his men while the coachman is already whipping the horses. We also note Dario Oxyathres, who sacrifices himself to allow his brother Darius III to save himself by letting himself be pierced by the Macedonian leader. The spears, the crowding of men and horses, the fallen horse and the Persian soldier in the foreground who looks at himself in agony in a mirror bring to mind the dramatic moment of the battle.

From the moment of its discovery, there was a debate between the court, the academicians, the archaeologists, the architects and the mosaicists regarding its state of conservation and the choice of keeping it at Pompeii in its original location or detaching and transferring it to the Royal Bourbon Museum. The conservation events of the Alexander Mosaic from the House of the Faun to its current location are described in detail in [28]: "At the end, after many uncertainties, the King resolved to transfer the Mosaic in Naples and on 16 November 1843, under pouring rain, the crate containing the Mosaic of Alexander was placed on a railway wagon coming from Maddaloni, near Caserta, pulled by sixteen oxen and escorted by armed soldiers, and after a grueling journey of nine days finally came in the Royal Museum where it was exhibited, not without controversy, in the ground floor of the Gallery, in the so-called Room of the Balbi. The Mosaic remained in this place until 1916 when it was moved to the western mezzanine where Vittorio Spinazzola, Director of the National Archaeological Museum of Naples between 1910 and 1924, picked up the mosaics in a new independent collection and where it still is exposed. A detailed documentation on photographic plates, made between 1916 and 1917, is preserved in the Historical Photo Archive of Special Superintendence for the Archaeological Heritage of Naples and Pompeii. The photographs gather information relative to all phases of the moving and the final installation". However, the type of structure created at the base of the mosaic to support it in a vertical position is not well understood as there are no shots of it.

The events briefly explained profoundly affected the state of conservation of the work, which is currently affected by important instability phenomena highlighted by the appearance of significant detachments of tiles (the occurrence is emphasized by the vibrations transmitted by car and railway traffic outside the museum), superficial lesions and swelling of the surface. In addition, it must be considered that the mosaic was designed and built to be placed on the floor, and consequently, the detachment from the site of discovery produced a significant alteration. In fact, the almost complete, albeit brief, removal of the rudus (the layer of lime, sand and aggregates, such as gravel or pebbles, which formed a concrete of consistent thickness) and the transport from the lower to the upper floor significantly affected the nucleus (the layer of mortar mixed with fragments of bricks that served as a support for the floor made of tiles). This procedure unnaturally transformed a floor mosaic into a wall mosaic [29].

#### **3. Material and Methods**

#### *3.1. Photogrammetric Digital Survey*

The production of the three-dimensional model of the mosaic surface was carried out through two photogrammetric surveys carried out in 2018 and 2020. The frames were acquired along horizontal bands (ensuring an average overlap of 80% both horizontally and vertically) at a distance of about 0.5 m from the mosaic surface using a Nikon D80 reflex camera (CCD sensor (23.6 × 15.8 mm) with 12.2 million pixels and a fixed focal length of 24 mm). The camera was set in aperture priority mode (f/9 value) and with an ISO sensitivity equal to 100. In order to reference the 3D model in Cartesian space in real metric units, a survey with a total station was also created using 163 clearly recognizable elements on the mosaic surface as markers. Therefore, a reference system was set up with an origin (P = 0,0,0) located in the lower-left corner of the mosaic and all the project data were implemented in an information system managed with Quantum GIS software (Version 3.16.11 Hannover, Open Source Geospatial Fundation, Beaverton, OR, USA).

Data processing was carried out using the Photoscan software (Agisoft Metashape Pro, Version 1.6.4, Agisoft LLC, St. Petersburg, Russia) through the following operations:


#### *3.2. Ground-Penetrating Radar (GPR)*

GPR is based on the diffusion of electromagnetic pulses into the soil and the recording of those re-radiated by buried targets characterized by sufficient dimensions and electromagnetic properties different from those of the surrounding ground. The quantities that are measured are the time required for the wave to travel the path from the transmitting antenna to a discontinuity and return to the surface (double time or two-way time) and the amplitude of the reflected wave. The double travel time depends on the speed with which the wave propagates within the material and provides information on the depth at which the reflectors are located. However, the amplitude, which represents how much energy returns to the surface after reflection, depends on the initial energy of the sent wave, the quantity that is dissipated along the way and the contrast of the electromagnetic properties of the materials that comprise the surface of the reflection. A complete description of the method is available in textbooks, such as [31–35].

The factors that influence the performance of the system in terms of the detectability of existing targets are the electromagnetic properties of the propagation medium, which determine the depth of investigation that can be reached, which varies from point to point. Since the attenuation of the means is a function of the radiated frequency, the use of highfrequency antennas generally allows for enhancing the resolution power, but at the expense of the penetration depth of the signals.

The georadar surveys were carried out on the entire surface of the mosaic in its actual vertical position. An IDS georadar RIS-K2 with a 3000 MHz high-resolution antenna was used for the data acquisition, whose features were considered sufficient and suitable for an overall assessment of the conservation status of the mosaic with respect to the type of target to be investigated and the thicknesses of surfaces to be analyzed.

Technically, data acquisition took place on 114 lines, from top to bottom every 5 cm, where instrumental readings were executed in continuous mode (Figure 3). In addition, several horizontal profiles and some profiles on the sides of the mosaic were acquired.

**Figure 3.** Grid data acquisition (**a**) with an IDS georadar (**b**).

Raw data were processed using GPR-SLICE 7.0 software (Screening Eagle Tecnologies Ag, Zurich, Switzerland) [36] using standard methodological approaches. In the first step, data and trace editing were executed, inserting information such as the temporal and spatial sampling intervals (time window, 20 ns; samples/scan, 512; scan/mark, 25; unit/marker, 1). Data were recorded and processed as 16-bit data and were converted by subtracting out the DC drift (wobble) in the data and, at the same time, adding a gain with time. A time-zero correction was determined to designate the starting point of the wave and the center frequency of the antenna was matched. Then, a bandpass filter and background removal were respectively applied to reduce noise from the oscillating components that had a regular frequency cycle in the frequency domain and to remove striation noises that occurred at the same time. Processed radargrams were subsequently corrected with an automatic gain control (ACG) function [37] that was applied to each trace based on the difference between the mean amplitude of the signal in the time window and the maximum amplitude of the trace. In the final step, the resulting filtered radargrams were inserted into a three-dimensional matrix from which sections were obtained at a particular double-time interval (measured in nanoseconds). Considering the complex layering and the reduced length of the profiles, we preferred to not apply the migration filter and we avoided presenting results converting time to depth using a mean value with the possibility to obtain an arbitrary and inaccurate estimate.

#### **4. Results and Discussion**

Figure 4 shows the DEMs obtained in 2018 and 2020 and Figure 5 shows some details of the last one transparently overlaid onto the processed orthophoto. The DEM analysis in 2018 (Figure 3a,b) highlighted depressions and reliefs relative to a hypothetical plane passing through the point of origin of the reference system (P = 0,0,0) in a spatial range of 3.66 cm. The reliefs were mainly located on the frame of the panel where the swellings responsible for the detachment of some mosaic tiles were clearly highlighted. The depressions, on the other hand, affected the central and right portions of the mosaic. The presence of an internal iron frame at the edges of the panel could, on the one hand, have favored the onset of bubbles and swelling as a consequence of the natural oxidation processes of the metal; on the other hand, it could have resisted deformations due to volumetric variations of the mortars in response to changes in temperature and humidity and/or deformations triggered by kinematic stresses (natural or artificial). This could explain the existence of depressions detected in the central part and the onset of lesions, especially on the edges (most resistant points) and the points of greatest weakness of the panel (central and right part). In order to identify the surface variations from 2018 to 2020, the point clouds of the surveys were metrically compared through the CouldCompare software (Verison 2.11 beta, Telecom ParisTech and the R&D division of EDF, Villeurbanne, France). Figure 6 shows the absolute distances between the point cloud obtained by subtracting the 2020 data from

those of 2018. On a global scale, a uniform swelling of the central part of the mosaic was estimated to be in the range of about 4 mm.

**Figure 4.** DEMs obtained in 2018 (**a**) and 2020 (**c**) and the same drawings transparently placed on the orthophoto (**b**,**d**). The values in the legend are expressed in meters.

**Figure 5.** Details of the orthophoto obtained in 2020 (**a**–**c**) and DEMs transparently overlaid on them (**d**–**f**).

**Figure 6.** Absolute distance obtained by subtracting the 2020 data from 2018 data.

Figure 7 shows the images obtained in 2018 from the visual analysis of the orthophoto in which the areas affected by restorations, the detachment of tesserae and fractures are highlighted. The indication of the fractures was just quantitative and therefore they were not classified by taking into account the elevation, thickness and dimensions. The visual analysis of the orthophoto made it possible to detect in detail every single fracture and the detachment of tesserae present on the surface. The latter phenomenon, caused by multiple factors, was aggravated by the vertical position of the mosaic, which, compared with the original horizontal position, determined the displacement on the ground and the irreversible loss of the individual tiles or portions of the mosaic. The fractures had two main alignments: one vertical/horizontal, which mainly affected the edges of the mosaic, and one diagonal to the mosaic, which concerned the central/right portion of the panel. The latter was particularly relevant and fell into the areas of greatest depression recorded in the digital surface elevation model. The presence of restored and consolidated areas on the left side may have increased the resistance of the surfaces, avoiding the onset of superficial lesions.

**Figure 7.** Restorations, detachments of tesserae and fractures reported on the orthophoto (**a**) and the DEM (**b**).

Regarding the GPR results, anomalies due to the presence of the iron support reinforcement present under the frame stood out in the different time slices (Figure 8). The degrees of amplitude variation in the time-slices were assigned a color scale that was chosen in order to show sufficient contrast to make the anomalies easily recognizable: light green corresponded to low amplitudes, while red corresponded to high amplitudes. The anomalies visible in these representations depicted the spatial distribution of the amplitudes of the reflections at specific depths within the grid. Within the sections, low amplitude variations expressed small reflections that indicated the presence of homogeneous material. High

amplitudes, on the other hand, denoted significant discontinuities in the investigated surface. In the most superficial part, we saw some anomalies of high amplitude that denoted a variation in the consistency of the materials (Figure 8b), which were located mainly on the left side of the mosaic. This aspect could be related to the presence of depressions and diagonal fractures in the central area and right side of the panel. The anomaly on the right side, near the border, was located at the point where there was swelling.

**Figure 8.** Details of the orthophoto obtained in 2020 (**a**) and slices relating to the time slices of 0.35–0.88 ns (**b**), 1.88–2.4 ns (**c**) and 3.34–3.87 ns (**d**).

In order to understand the inner structure, in the following part, the analysis of photographic documentation, made between 1916 and 1917, which records all phases of the moving and the final vertical installation, was jointly analyzed with GPR radargrams and slices. Following the images, it is known that the mosaic was housed on some blocks of rock. In Figure 9a, it is possible to note the original stratigraphy by recognizing, from top to bottom, the nucleus, the rudus and the screed covering the stone beams. Subsequently, the section was cleaned, the screed was removed and the stone structure was brought to light. Some wedges separated the rudus from the blocks of rock (Figure 9b). Then, the surface of the section was leveled and the rudus was cleaned of residual lime for the realization of the support frame (Figure 9c). In the next phase, equally spaced holes were drilled in the rudus, three of which are visible in Figure 9d.

From the analysis of the georadar results, it was assumed that these holes were drilled to anchor the rudus to the frame under construction by means of pins. Figure 10b shows the acquired radargram at the edge of the frame, as indicated by the green arrow on the mosaic orthophoto (Figure 10d). With red circles, the positions of the holes visible in the photos are indicated. At those points, hyperbolas are attributable to the presence of iron pins. However, there are other anomalies (indicated with yellow arrows) of the same nature. At the edges, the close hyperbolas reflected the positions of the nails used to fix the iron L-bars (blue arrows) installed in the following stages. Figure 10c shows the radargram acquired at the edge of the original frame, as indicated by the magenta arrow on the mosaic. The anomalies persisted and another one was highlighted on the left (red arrow). There were no more anomalies in the upper band. Then, the pins penetrated the mosaic by at least

15 cm (a few cm in the original one). The radargram acquired in the upper band, again at the edge of the original mosaic, highlighted the same anomalies found below (Figure 10a). Here, the pins are still visible from the top of the frame.

**Figure 9.** Photo no. 4138-1916-EX 306 (**a**), photo no. 4131-1916-EX 299 (**b**), photo no. 4129-1917-EX 297 (**c**) and photo no. 4130-1916-EX 298 (**d**).

**Figure 10.** Location of the anchor pins (indicated by red circles) and further anomalies on three horizontal radargrams (**b**–**d**) located at the positions of the colored arrows (**a**).

The same situation is found on the sides where two hyperbolas are visible (Figure 11), in addition to the pins visible with the naked eye.

Figure 12 shows the digital model of the processes followed before applying the frame and removing the stone.

**Figure 11.** Orthophoto of the left and right sides of the mosaic and radargrams labeled 01 and 02. Red dots indicate visible pins.

**Figure 12.** Digital model with the locations of the blocks of rock and the anchor pins.

The external frame was probably made using wooden crosspieces measuring 25 × 13 cm. Figure 13a shows a detail of Figure 9d (photo no. 4130-1916-EX 298), where in the background, a wooden crosspiece appears with regular carvings. They were assumed to be about 10 cm wide and placed at intervals of about 22 cm. Figure 13b shows the digital model of the carved frame. It was assumed that this element is related to the underlying part of the longitudinal crosspiece that was specially prepared to house the elements of the internal filling.

**Figure 13.** Photo no. 4130-1916-EX 298 (**a**) and digital model of the frame with indication of the notches using red arrows (**b**).

The anchoring of the wooden crosspieces with the mosaic is described in the previous section. To reinforce the structure, four angular iron bars were added (Figure 14). Holes were also drilled to position the protective covering of the front part of the mosaic. A 2 cm board was placed under the wooden crosspieces, which were placed under the stone crosspieces (Figure 14a). This table was probably removed in the next phases.

**Figure 14.** Photo no. 4132-1916-EX 300 (**a**), a digital model with the anchoring of the frame and addition of iron angles (**b**) and a cross-section (**c**).

A protective sheet was added to the surface, which was blocked with a board on the edge by means of nails in the pre-drilled points in the previous step (Figure 15a). Some sides were added (10 cm wide and 5 cm high) in between where the wide boards were to be placed (Figure 15b). The gaps were subsequently sealed with mortar (Figure 15c). After

which, the lifting of the recessed mosaic began (Figure 14d), which was placed in a vertical position on a support (Figure 15e).

**Figure 15.** Photo no. 4133-1916-EX 301 (**a**), photo no. 4134-1916-EX 302 (**b**), photo no. 4135-1916-EX 303 (**c**), photo no. 4157-1917-EX 325 (**d**) and photo no. 4173-1917-EX 341 (**e**).

On the back, from a historical photo, an alternation of narrow dark boards and wide light boards, 10 cm and 22 cm wide, respectively, can be noted (Figure 16a). In the GPR slice relating to the time window 1.88–2.4 ns (Figures 8b and 16b), a series of vertical anomalies were highlighted that began directly in contact with the rudus. It was plausible that the wooden structure placed at the rear of the mosaic traced this structure in some way. In particular, it is likely that the narrow boards were actually beams or murals 10 cm wide (with a maximum length of 14 cm, depending on the thickness of the shaving made under the rudus), which rested directly under the mosaic. As can be seen from many of the radargrams, some profiles seemed to be characterized by the presence of metal elements as brackets to block the vertical wooden beams from the rudus. An example is shown in Figure 17, where three identifiable hyperbolas emerged clearly with the rudus, the filling lime and the closing table. Hyperbolas 1 and 2 (Figure 17a) often penetrated the mosaic and probably represented a system for anchoring the possible wooden beam (vertical georadar anomalies) directly to the mosaic. At approximately hyperbola 3, a sort of reinforcement was evident under the vertical structures (Figure 18). The spaces were probably filled with plaster and closed everything with the boards that were visible on the back. The slice (Figure 17b), which was related to the time window between 3.34 ns and 3.87 ns, was enhanced to visualize the three anomalies jointly, even if detected at different depths.

**Figure 16.** Historical photo (photo no. 4172-1917-EX 340) straightened and positioned on the back of the mosaic (**a**) and time slice relating to the time window 1.88–2.4 ns (**b**).

**Figure 17.** Main anomalies, labeled 1, 2 and 3, indicated on the radargram (acquired as indicated by the magenta arrow) (**a**) and on the slice relating to the time window 3.34–3.87 ns (**b**).

**Figure 18.** Anchor system digital model.

Therefore, after the removal of the blocks of rock, the vertical beams would have been housed (directly or indirectly) thanks to the notches located at the base of the longitudinal beams of the frame. In this system, it is possible to make four hypotheses based on the presumed depth of the notch, the location of the vertical beams and the closing plank. The hypotheses in question took into account the position and shape of the acquired georadar anomalies:

**Hypothesis 1.** *(Figure 19): Support boards (25* × *12 cm) were inserted at regular intervals under the rudus. The boards were closed in the lower part with boards (10* × *2 cm) that protrude by 2 cm outside. The empty space between the support boards was filled with plaster. We consider this hypothesis the least probable, as it does not explain the presence of the georadar anomalies previously discussed.*

**Figure 19.** Hypothesis 1 digital model: locations of support boards under the rudus (**a**), complete support system (**b**), cross-section (**c**) and longitudinal section (**d**).

**Hypothesis 2.** *(Figure 20): The vertical beams inserted at regular intervals were worked to be housed and nailed to the frame. Between one beam and the other, closing tables were inserted through brackets. The gaps left empty between the beams were filled with plaster.*

**Figure 20.** Hypothesis 2: locations of vertical beams (**a**), locations of anchoring brackets (**b**), locations of closing boards (**c**), complete support system (**d**), cross-section (**e**) and longitudinal section (**f**).

**Hypothesis 3.** *(Figure 21): This differs from hypothesis 2 in that the incision was made to be 12 cm in order to completely contain the vertical beam that rests directly under the rudus and protrudes 2 cm outside.*

**Figure 21.** Hypothesis 3 digital model: positions of vertical beams (**a**), locations of anchoring brackets (**b**), arrangement of closing tables (**c**), complete support system (**d**), cross-section (**e**) and longitudinal section (**f**).

**Hypothesis 4.** *(Figure 22): The planks were anchored to the vertical beams by means of brackets according to a "T" system. Closing tables were inserted between one plank and the other. The gaps left empty between the beams were filled with plaster. This is the hypothesis that best matches the shape and position of the georadar anomalies.*

**Figure 22.** Hypothesis 4 digital model: vertical beam positions (**a**), anchoring bracket locations (**b**), board arrangements (**c**), external board positions (**d**), complete support system (**e**), cross-section (**f**) and longitudinal section (**g**).

#### **5. Conclusions**

In this study, the issue connected to the conservation of the Alexander Mosaic was investigated in an organic and exhaustive way, starting from the analysis of the available knowledge framework. From an operational point of view, the adopted research methodology provided for the recognition of all available documentation, an accurate survey of the current state from the morphological point of view and non-invasive diagnostic investigations. This allowed for obtaining a further preliminary and in-depth knowledge phase of the asset to also provide a framework for comparison with what emerged from previous studies and, therefore, data on any changes in the state of conservation of the asset. In detail, in addition to a slight progression of the degradation phenomena already recorded in 2018, regarding the lesions, anomalies in the support emerged from the diagnostic results, which suggested the presence of discontinuities in the background mortar. This could have probably been due to the presence of different materials, with particular reference to the mosaic embedding system put in place during the transfer in 1916, consisting of a wooden frame placed along the perimeter, with interposed wooden joists reinforced with metal elements. However, to date, insufficiently exhaustive data has emerged regarding the state of conservation of the support that is not visible and cannot be inspected directly and the actual state will be fully detected only after moving the mosaic from its current position, making the back accessible. Currently, the first phase of the restoration can now be considered completed, which was limited to interventions to make the mosaic surface safe, which was necessary in the face of the detected degradation phenomena and limited only to the elimination of the conditions that can generate critical issues during movement (loose tiles, presence of continuity solutions between the layers, etc.). Subsequently, the mosaic will be overturned, following which it will be possible to access the support to

promptly check its conservation status, fine-tune the subsequent conservation interventions and possibly provide for the construction of a new collection system should the existing one prove to be no longer suitable. Subsequently, the mosaic will be placed in a horizontal position that will allow for the restoration of the mosaic surface and the final rearrangement to be carried out. The restoration operations, which will also affect the rear parts, will be preceded by further diagnostic analyzes once the mosaic is overturned. In this way, at the end of the project, a complete model will be obtained, which will integrate all the results and will represent a precious source of knowledge on this important artifact.

**Author Contributions:** Conceptualization, M.C., P.M., A.D.S. and A.P.; Data curation, M.C., P.M. and V.G.; Methodology M.C., P.M. and V.G.; Validation, M.C., P.M., A.D.S., A.P. and V.G.; Writing original draft, M.C., P.M. and V.G.; Writing—review and editing, M.C., P.M. and V.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** Special thanks to Paolo Giulierini for having involved us in this important project and for having been a constant guide, always demonstrating the greatest passion and dedication to the protection of this very important masterpiece. We are also very thankful to the company "Boviar. Integrated systems for diagnostics and monitoring", in particular Filippo Bovio, for the technical-scientific support during the data acquisition and processing phases.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Tung-Ching Su \*, Tsung-Chiang Wu, Ming-Hung Wun and Cheng-Wei Wang**

Department of Civil Engineering and Engineering Management, National Quemoy University, Kinmen 89250, Taiwan; tsung\_chiang@nqu.edu.tw (T.-C.W.); ji394su3love20@gmail.com (M.-H.W.); joswawagiga@gmail.com (C.-W.W.) **\*** Correspondence: spcyj@nqu.edu.tw

**Abstract:** Many studies in the literature have presented multiple remote sensing techniques for defect inspection of paintings. At present, however, papers on defect inspection and restoration of oriental architectural arts—such as door god paintings—are still rare. If an aged and damaged door god painting needs a restoration, then following the style and treatment skill of the original artist as much as possible is important for the restoration. Unfortunately, it is usually difficult to access the original artists for some of the aged door god paintings. This paper considers the texture features of auspicious patterns of armors on warrior door gods as useful information to recognize styles of door god paintings by unknown artists. First, a two-level two-dimensional discrete wavelet transform coupled with co-occurrence matrix calculation was adopted to analyze the texture features, based on the descriptors of angular second moment (ASM), entropy (ENT), contrast (CON), homogeneity (HOM), dissimilarity (DIS), correlation (COR), and cluster tendency (CLU), in the four orientations of 0◦ (horizontal), 45◦ (vertical), and 90◦ and 135◦ (double diagonal). Second, a two-tailed *t*-test based on the analyzed texture features was introduced into the hypothesis testing for demonstrating the master and apprentice relationships between the surveyed artists, and for recognizing the door god painting styles of unknown artists as well. The experimental results show that the proposed method effectively describes the texture features of the auspicious patterns of the surveyed door god paintings, and is able to determine the useful co-occurrence features for recognizing unknown artists' painting styles.

**Keywords:** door god paintings; texture features; discrete wavelet transform; two-tailed *t*-test; painting style recognition

#### **1. Introduction**

Door god paintings are frequently seen on the gates of traditional residences or temples, and belong to one kind of Chinese architectural artworks (see Figure 1 [1]). Because of Chinese people's worship of animism in ancient times, some of the door god characters are imaginary. Today, there are several pairs of door gods representing faith in home safety and evil avoidance. The noun "door god" can be seen in the period as early as before the Qin Dynasty [2]. Essentially, the door god characters can be separated into two categories: warrior, and civil official. The greatest difference between warriors and civil officials is that a set of armor is necessary for a warrior, but unnecessary for a civil official. Shen Tu (Chinese transliteration: 神荼) and Yu Lei (Chinese transliteration: 鬱壘) are the earliest warrior door gods and prevailed in the Han Dynasty [3]. In addition to Shen Tu and Yu Lei, Chin Shu Pao (Chinese transliteration: 秦叔寶) and Yu Chih Kung (Chinese transliteration: 尉遲恭) are also frequently seen warrior door gods in the temples of Taiwan. Several auspicious patterns—such as lock chain, turtle back, fish scales, flowers, clouds, and the Chinese characters of "回" and "卍" (see Figure 2)—have been adopted to decorate warriors' armor. For each of the above auspicious patterns, the displayed texture features

**Citation:** Su, T.-C.; Wu, T.-C.; Wun, M.-H.; Wang, C.-W. Style Recognition of Door God Paintings by Hypothesis Testing for Texture Features of Painting Patterns. *Appl. Sci.* **2022**, *12*, 2637. https://doi.org/10.3390/ app12052637

Academic Editor: Asterios Bakolas

Received: 16 February 2022 Accepted: 1 March 2022 Published: 3 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

may be different due to the creativity of different artists. The different texture features mean that every artist has his/her own unique painting style via his/her special treatment or skill with door god paintings [4].

Any famous door god painting generally has a long history, but that long history may result in the impossibility of identifying the original artist. Additionally, door god paintings use residence or temple gates as their support materials, and so the placed locations of door god paintings belong to a semi-outdoor environment, where door god paintings are apt to suffer as a result of climate factors—such as temperature, relative humidity, and wind—or human factors, such as ritual activities and inadvertent collisions. Hence, the conservation and maintenance of door god paintings are more challenging than for other paintings or artworks collected in a well-controlled, indoor environment, such as an art gallery or museum. When restoring an aged and damaged door god painting, the treatment skill and the procedure must follow its original style as best as possible in order to preserve the unique style of the door god painting. In case of the original artist being unknown, the original style cannot be correctly followed by the painting restoration.

**Figure 1.** A pair of door god paintings: (**a**) Yu Chih Kung (Chinese transliteration: 尉遲恭): warrior on the left door; (**b**) Chin Shu Pao (Chinese transliteration: 秦叔寶): warrior on the right door [1].

**Figure 2.** Example images of auspicious patterns on paintings of warriors' armor: (**a**) lock chain pattern; (**b**) turtle back pattern; (**c**) fish scales pattern; (**d**) flowers pattern; (**e**) clouds pattern; (**f**) Chinese character pattern: "回"; (**g**) Chinese character pattern: "卍".

#### **2. Research Aim**

The types of deterioration of door god paintings usually include lifting, cracks, shrinkage, atomization, chalking, etc. As far as we know, the current literature discussing restorations of door god paintings is quite scarce. Related image processing methods have also demonstrated their effectiveness in artwork defect detection by providing useful information for artwork restoration [5–7]. This paper's main contributions are as follows: First, to apply an image-processing-based texture analysis method to extract the texture features from the acquired auspicious pattern images. Second, hypothesis testing is implemented on the texture feature dataset to see the similarity of the auspicious patterns between the surveyed artists, and to demonstrate the master and apprentice relationship as well. Finally, hypothesis testing also is applied to the texture feature dataset to see to which painting style(s) among the surveyed artists the door god painting style(s) of the unknown artist(s) are the closest. In this way, a door god painting restoration can retain its originality as much as possible.

#### **3. Backgrounds**

#### *3.1. Artists of Door God Paintings in Taiwan*

Since around the 1960s, there have been certain famous artists who created door god paintings for traditional residences or temples in Taiwan. However, most of these famous artists are now deceased, and so their created door god paintings are regularly restored by their apprentices, some of whom are even members of the famous artists' family. For instance, the two famous artists of door god paintings—Yu-Feng Chen (Chinese transliteration: 陳玉峰, 1900–1964) and Chun-Yuan Pan (Chinese transliteration: 潘春源, 1891–1972)—were both apprentices to an artist from Quanzhou, Fujian, China, and most of their door god paintings are distributed over temples in southern Taiwan. Shou-Yi Chen (Chinese transliteration: 陳壽彝, 1934–2012, the eldest son of Yu-Feng Chen) and his cousin, Tsao-Ju Tsai (Chinese transliteration: 蔡草如, 1919–2007), were both apprentices of door god painting to Yu-Feng Chen. Among the door god paintings of Yu-Feng Chen, Shou-Yi Chen, and Tsao-Ju Tsai, there is a special painting style belonging to their family genre. Similarly, Chun-Yuan Pan, Li-Shui Pan (Chinese transliteration: 潘麗水, 1914–1995, the son of Chun-Yuan Pan), and Yueh-Hsiung Pan (Chinese transliteration: 潘岳雄, 1943–present, the eldest son of Li-Shui Pan) also created door god paintings with a special painting style belonging to their family genre. Consequently, a mentorship genealogy concerning the creation of door god paintings in southern Taiwan can be established. As for artists in other areas of Taiwan, Hung [2] established a mentorship genealogy for a famous artist in Hsinchu, Taiwan, and compiled a list of the door god paintings created by the artist and his apprentices.

There have been a few cases of door god painting restoration in Taiwan in the past. In 1968, S.-Y. Chen cooperated with two other artists to refurbish aged door god paintings of the Qingshui Zushi (Divine Ancestor) Temple in New Taipei City, Taiwan [4]. Another artist, Lien-Cheng Hsu (Chinese transliteration: 許連成, 1919–2002), created door god paintings for a historic temple in northern Taiwan in 1975, and Hsu's door god paintings were refurbished in 1989 by Chia-Cheng Liu (Chinese transliteration: 劉家正, 1955–present). In fact, C.-C. Liu refurbished many of Hsu's door god paintings, but Liu's painting style belongs to the family genre of Pan rather than Hsu, and so there is a little controversy about the restoration approach.

The original artists of the door god paintings in Taiwan can almost always be identified but, unfortunately, in Kinmen and Penghu—which are the outlying islands of Taiwan—the original artists of door god paintings are mostly difficult to ascertain. Along with the immigration of ancestors from southern Fujian (China) to Taiwan, Kinmen and Penghu were selected as the relay stops. Therefore, there certainly should be some kind of master and apprentice relationships between the artists in southern Fujian and Taiwan, including the outlying islands.

#### *3.2. Texture Feature Analysis for Paintings*

An image showing a frequency of tonal change, such as zebra crossings and knitting, can be regarded as one with texture features [8]. Similarly, an auspicious pattern image also shows the frequency of tonal change. Several studies in the literature have discussed the applications of texture feature analysis to painting restoration [9–11], painting classification by artistic genre [12–14], painting style recognition [15,16], and raw material investigation [11,17]. Cai and Siegel [9] explored two texture features—energy (or homogeneity) and entropy (or disorderliness)—to model the visual appearance of paintings before and after surface varnishing. Based on a calculation of the gray-level co-occurrence matrix (GLCM), energy and entropy demonstrated that they are sensitive to the effect of varnishing—especially as entropy shows the increase in a painting's contrast by varnishing. The statistic calculation of GLCMs is a common approach, and has been used to extract texture features for the classification of traditional Chinese paintings based on the painters [18].

Wavelet-based analysis, where multiscale and multiorientation image decomposition is performed, has been applied to a collection of high-resolution digital scans of drawings or paintings to describe the painting characteristics [19]. Cetinic and Grgic [20] introduced various classifiers into automated painter recognition based on the texture features extracted by GLCM and discrete wavelet transforms (DWTs). Undoubtedly, the statistical properties of wavelet coefficients have demonstrated that they are successful in the stylistic analysis of paintings [21,22]. Referring to related research [23,24], the process of texture feature extraction first conducts image decomposition using DWT to obtain sub-band images that are multiscale and multiorientation. Secondly, GLCMs are generated from sub-band images, and finally the texture descriptors—such as energy, entropy, contrast, correlation, and others—are calculated from GLCMs.

#### **4. Research Material**

Based on the investigations of Lee [25] and Kang [26] on door god paintings in Taiwan, some of their investigated temples and other temples in Taiwan were selected as the study sites, and smartphones (model ASUS\_Z01HDA or Sony\_J9110) were used to acquire the auspicious pattern images. In total, 52 temples were visited and 453 auspicious pattern images were acquired (see the zip folder in the Supplementary Materials). Table 1 lists the number of acquired auspicious pattern images, which involve fish scales, lock chain, turtle back, and "回" patterns by the 31 artists. Among the four auspicious patterns, lock chain and fish scales are the most common patterns adopted in door god paintings, and so they have a much greater number of acquired images than the others. Moreover, the stronger the artist's reputation, the greater the number of auspicious pattern images acquired.

In addition to the 453 auspicious pattern images of the surveyed artists, we also acquired 12 auspicious pattern images of unknown artists from temples in Kinmen for the study of painting style recognition. In order to avoid blurred imaging of the door god paintings, support of the smartphones by a tripod instrument instead of a handheld approach was necessary. Considering the different sizes of the door god paintings, the appropriate instrument height and object distance had to be determined by trial and error during the imaging process. Most of the trial and error results in the field indicate that the appropriate instrument height and object distance are ~1.5 m and 1–2 m, respectively. Because the door god paintings are located in semi-outdoor environments, sometimes the imaging process suffers as a result of the structures of the temples, the azimuths of sunshine, or the weather, leading to an inappropriate illuminated image. An inappropriate illumination will result in an image with low quality in contrast and brightness. The f-stop numbers (apertures) of the ASUS\_Z01HDA and Sony\_J9110 are *f*/1.7 and *f*/1.6 (*f*: lens focal length), respectively. Thus, an automated adjustment of lens focal length was adopted in order to acquire the auspicious pattern images with appropriate contrast and brightness.


**Table 1.** Numbers of auspicious pattern images, with their corresponding artists.

#### **5. Methodology**

#### *5.1. Texture Feature Extraction*

Figure 3 shows a diagram of the texture feature extraction proposed in this paper. Before performing the 2-level 2-dimensional (2D) DWT, each RGB-wise auspicious pattern image is transformed into a grayscale one. In other words, the color information is not taken into consideration for our DWT operation. After the 2-level 2D DWT, an approximation image and three detailed images in horizontal (0◦), vertical (90◦), and diagonal (45◦ or 135◦) orientations are derived. Based on the three detailed images, the GLCMs in the 4 orientations are calculated. Hereafter, each GLCM calculates the co-occurrence features, consisting of angular second moment (ASM), entropy (ENT), contrast (CON), homogeneity (HOM), dissimilarity (DIS), correlation (COR), and cluster tendency (CLU). In this research, the co-occurrence features in a certain orientation can be expressed as a 1 × 7 vector, and so the co-occurrence features in the 4 orientations are finally integrated into a 4 × 7 matrix.

**Figure 3.** Scheme of texture feature extraction from auspicious pattern imagery.

#### 5.1.1. Two-Level 2D DWT

Wavelet transform (WT) is a linear transform developed from Fourier transform, where the basic functions are sinusoids, but the wavelet functions vary with frequency and limit duration, thus offering better resolutions along the frequency scale [27,28]. In the beginning, the development of WT was for signal processing, but not for image processing, until Daubechies and Mallat provided the discretization of WT and established the connection between WT and the multiresolution theory, respectively. This paper regards an auspicious pattern image as the change in a discrete signal along a 2D scale. A 2D DWT can decompose the auspicious pattern image into many scales, which range from the roughest scale to the finest [29–31]. Through a decomposition of 2D DWT, which is implemented by consecutive low-pass (L) and high-pass (H) filtering through one-dimensional convolution, the auspicious pattern image I(*m*, *n*) can be divided into sub-band images, including an approximation image (LL) and three detailed images in horizontal (HL), vertical (LH), and diagonal (HH) orientations [28]. The parameters of *m* and *n* denote the number of pixels of the image in the row and column directions, respectively.

The approximation image (LL) can be further decomposed in the next level to obtain the images of LL, HL, LH, and HH, with the sizes of (*m*/2, *n*/2). For the decomposition of each level, the alternative convolutions of the approximation image and the low-pass or high-pass filters in the column or row directions are operated using a downsampling function by 2. Thus, the 2D DWT finally produces the pyramid representations of the sub-band images, which occupy the same amount of storage as the original auspicious pattern image [32].

Several studies in the literature indicate that the 2-level 2D DWT (i.e., 2D DWT in the second level) can robustly extract the details of the image texture features from the original images [23,24,28,33]. After decomposition in the second level, the 7 sub-band images, including 6 detailed images (HL1, LH1, HH1, HL2, LH2, and HH2) and one approximation image (LL2), are obtained. The indices of "1" and "2" signify the first and second levels, respectively. Among the above 6 detailed images, HL2, LH2, and HH2 were chosen to calculate the paper's GLCMs.

#### 5.1.2. Calculation of the Gray-Level Co-Occurrence Matrix (GLCM)

A GLCM is a one-dimensional square matrix, where each element in the row (*i*) and column (*j*) directions records a relative occurrence frequency (*Mij*) of a pair of pixels with the same element value, separated by a certain pixel distance (*D*) in one orientation (θ) [23,28]. Figure 4 shows an illustration of GLCM production at different orientations and pixel distances for a pixel of interest in the square matrix, along with an example of statistics for *Mij* at 90◦ (vertical direction) and one pixel distance. Figure 4b shows the GLCM of an example detailed image, where the image size is 7-by-7 (but not limited to square) and the number of gray levels is 8. In the GLCM, the element (3, 0) has a value of 1, because there is only one instance in the example detailed image where the vertically adjacent pixels (2, 1) and (1, 1) have gray-level values of 3 and 0, respectively.

**Figure 4.** Illustration of GLCM production: (**a**) A pair of pixels with the same element value separated by the different pixel distances in 4 orientations. (**b**) Relative occurrence frequency recording of gray levels of digital images, using a GLCM with *D* = 1 and θ = 90◦.

**E**

Except for the four necessary orientations of 0◦, 45◦, 90◦, and 135◦, several studies have indicated that the different pixel distances (*D*) have a great influence on the accuracy of GLCM-based texture description, and demonstrated that one pixel distance (*D* = 1) can lead to better accuracy of the feature extractions [23,34]. This paper introduces an adjacent distance of bordering on pixels instead of striding over pixels into the calculation of the GLCMs.

#### 5.1.3. Calculation of Co-Occurrence Features

The GLCMs obtained from the detailed images of HL2 and LH2 further calculate the co-occurrence features in 0◦ and 90◦, respectively. Both of the GLCMs in 45◦ and 135◦ obtained from the detailed image of HH2 further calculate the co-occurrence features in the double-diagonal orientations. The 7 aforementioned co-occurrence features can be obtained using the following texture descriptors:

Angular Second Moment (ASM): The uniformity of distribution of the gray level in an auspicious pattern image can be represented thus [35]. The texture descriptor of ASM is expressed as follows:

$$\sum\_{i=1}^{n} \sum\_{j=1}^{n} P\_{ij}^{2} \tag{1}$$

The value of ASM ranges from 1/*n*<sup>2</sup> to 1. When the value is 1, it means a constant image. Here, *Pij* is calculated as *Mij* ∑*n <sup>i</sup>*=<sup>1</sup> <sup>∑</sup>*<sup>n</sup> <sup>j</sup>*=<sup>1</sup> *Mij* .

Entropy (ENT): A statistical measure of randomness determines the textural interference in an auspicious pattern image [24,35]. The texture descriptor of ENT is obtained as follows:

$$-\sum\_{i=1}^{n}\sum\_{j=1}^{n}P\_{ij}\cdot\log\_{10}^{P\_{ij}}\tag{2}$$

The larger the ENT value, the higher the textural complexity of the auspicious pattern. Thus, there should be a highly negative correlation between ASM and ENT.

Contrast (CON): The local variations in an auspicious pattern image can be measured by this descriptor, which is calculated as follows:

$$\sum\_{i=1}^{n} \sum\_{j=1}^{n} \left(i - j\right)^{2} \cdot P\_{ij} \tag{3}$$

The higher the CON, the higher the image contrast will be. Hence, a CON of 0 means a constant image.

Homogeneity (HOM): HOM can also be called an inverse differential moment, which measures the similarity between the distributions of elements in a GLCM and those in the diagonal GLCM [24]. The textural descriptor is defined as follows:

$$\sum\_{i=1}^{n} \sum\_{j=1}^{n} \frac{P\_{ij}}{1 + (i - j)^2} \tag{4}$$

The representation of HOM is contrary to that of CON—the higher the HOM, the lower the image contrast (or the more the image homogeneity) will be.

Dissimilarity (DIS): The degree of dissimilarity of the gray levels in an auspicious pattern image is now measured. DIS is very sensitive to the arrangement of gray-level values or tones in an imagery space, and is expressed as follows:

$$\sum\_{i=1}^{n} \sum\_{j=1}^{n} |i - j| P\_{ij} \tag{5}$$

A higher DIS value means a greater dissimilarity of gray levels in the image.

Correlation (COR): The spatial dependencies between the image pixels can be defined [24] to indicate the textural directionalities of the auspicious patterns [36] by this textural descriptor, which is calculated as follows:

$$\frac{\sum\_{i=1}^{n} \sum\_{j=1}^{n} (i \cdot j) P\_{ij} - \mu\_x \mu\_y}{\sigma\_x \sigma\_y} \tag{6}$$

In other words, COR is a correlation coefficient of GLCM. The higher the COR is, the greater is the spatial dependency of gray levels in a certain direction.

Cluster tendency (CLU): The degree of the textural clustering in an auspicious pattern image is measured and can be obtained as:

$$\sum\_{i=1}^{n} \sum\_{j=1}^{n} \left( i - \mu\_{\mathcal{X}} + j - \mu\_{\mathcal{Y}} \right)^{2} \cdot P\_{\hat{\mathcal{Y}}j} \tag{7}$$

It should be noted that a lower CLU value means higher textural clustering. Moreover, there is an interrelationship between CLU and COR. When the values of CLU and COR are large, there will be an indefinite textural directionality or a wide distribution of gray levels. Here, *μx*, *μy*, *σx*, and *σ<sup>y</sup>* are calculated as follows:

$$\begin{cases} \begin{aligned} \mu\_{\boldsymbol{x}} &= \sum\_{i=1}^{n} \sum\_{j=1}^{n} i \cdot P\_{ij} \\ \mu\_{\boldsymbol{y}} &= \sum\_{i=1}^{n} \sum\_{j=1}^{n} j \cdot P\_{ij} \\ \sigma\_{\boldsymbol{x}} &= \sqrt{\sum\_{i=1}^{n} \sum\_{j=1}^{n} (i - \mu\_{\boldsymbol{x}})^2 \cdot P\_{ij}} \\ \sigma\_{\boldsymbol{y}} &= \sqrt{\sum\_{i=1}^{n} \sum\_{j=1}^{n} (j - \mu\_{\boldsymbol{y}})^2 \cdot P\_{ij}} \end{aligned} \tag{8}$$

#### *5.2. Statistical Testing*

This paper adopted a two-tailed *t*-test to test the similarity between any two auspicious patterns, the similarity of an auspicious pattern between any two surveyed artists, or the similarity of an auspicious pattern between the surveyed and unknown artists. The hypotheses of the *t*-tests are described as follows:

#### 5.2.1. Testing for Similarity between Any Two Auspicious Patterns

This paper aims to examine the similarities and differences between the auspicious patterns in door god paintings in Taiwan by applying hypothesis testing to the obtained co-occurrence features. As shown in to Figure 3, we calculated the co-occurrence feature values for each auspicious pattern image, which were recorded as a 4-by-7 matrix. Thus, the obtained co-occurrence feature values in a certain orientation can be arranged as an *n*-by-7 matrix, where *n* is the number of auspicious pattern images. In the fish scales pattern, for example, the data size of its obtained co-occurrence feature values in one orientation is 245-by-7. The arrangement of the 245-by-7 matrix also means that the hypothesis testing does not consider the different painting styles between the artists. For each orientation, any 2 of the 4 auspicious patterns performed a two-tailed *t*-test, where the null (H0) and alternative (H1) hypotheses were expressed as follows:

$$\begin{cases} \text{ } \mathcal{H}\_0: D\_i = D\_j\\ \text{ } \mathcal{H}\_1: D\_i \neq D\_j \end{cases} \tag{9}$$

and *Di* and *Dj* are the *n*-by-7 data matrices of the extracted co-occurrence features of auspicious patterns *i* and *j*, respectively. Thus, there is a total of C4 <sup>2</sup> = 6 combinations for the two-tailed *t*-test.

#### 5.2.2. Testing for Similarity of an Auspicious Pattern between Any Two Surveyed Artists

This paper demonstrates whether any two surveyed artists have a master and apprentice relationship by applying hypothesis testing to the extracted co-occurrence features of the auspicious patterns of their door god paintings. For some artists, unfortunately, there are extremely few images accessibly acquired, resulting in a non-normal distribution for the number of tested samples. In order to control the risk of accepting a false hypothesis, both α and β—which are the probabilities of rejecting and accepting the null hypothesis, respectively—must be considered. For this paper, the values of α/2 = 0.05 and β = 0.1 were introduced into the calculation of the minimum sample size (N) from the National Institute of Standards and Technology [37]:

$$\mathbf{N} \ge \left(\mathbf{t}\_{1-\frac{\alpha}{2}} + \mathbf{t}\_{1-\beta}\right)^2 \left(\frac{\mathbf{s}}{\delta}\right)^2 \tag{10}$$

where s and *δ* denote the standard deviations of the samples and the population, respectively. In this paper, s is equal to *<sup>δ</sup>*, so Equation (10) can be rewritten as t1<sup>−</sup> <sup>α</sup> <sup>2</sup> + t1−<sup>β</sup> 2 . Considering the degrees of freedom as 30, the values of t1<sup>−</sup> <sup>α</sup> <sup>2</sup> and t1−<sup>β</sup> are 1.6973 and 1.3104, respectively; thus, the *N* value approximates to 9.

According to the above estimation for the minimum size of samples, this paper considers the surveyed artists, from whom we can acquire more than 10 images (including 10 images) for some auspicious patterns, using a two-tailed *t*-test. Therefore, from Table 1, only the fish scales and lock chain patterns satisfy the above condition, and there is a total of C*na* <sup>2</sup> (*na*: number of considered artists) combinations for the two-tailed *t*-test. In this paper, the *na* values are 7 and 6 for the fish scales and lock chain patterns, respectively. The cooccurrence features of the fish scales and lock chain patterns in an orientation are arranged as an *m*-by-7 data matrix, where *m* (*m* ≥ 10) is the number of acquired auspicious pattern images for a certain considered artist. Based on Equation (9), here *Di* and *Dj* represent the *m*-by-7 data matrices of artists *i* and *j*, respectively.

5.2.3. Testing for Similarity of an Auspicious Pattern between Surveyed and Unknown Artists

The two-tailed *t*-test was applied to the extracted co-occurrence features of the fish scales or lock chain patterns of the unknown artists' door god paintings to determine whether there was a similar painting style between the surveyed and unknown artists. Among the acquired 12 auspicious pattern images of the unknown artist(s), the number of fish scales and lock chain pattern images was half each. Hence, the sample of just 6 images is able to control the risk of accepting a false hypothesis by the two-tailed *t*-test, with α/2 = 0.05, β = 0.05, and *δ* = 1.5s. In the two-tailed *t*-test, an *l*-by-7 (*l* = 6 in this paper) data matrix was built for recording the co-occurrence features of the fish scales or lock chain patterns in an orientation, where *l* is the number of acquired fish scales or lock chain pattern images for the unknown artist(s). Based on Equation (9), here *Di* and *Dj* represent the *m*-by-7 and *l*-by-7 data matrices of surveyed artist *i* and unknown artist *j*, respectively.

#### **6. Results and Discussion**

#### *6.1. Texture Features of Surveyed Auspicious Patterns*

The four auspicious patterns in Table 1, including fish scales, lock chain, turtle back, and "回", were surveyed and calculated for their co-occurrence features (texture features), as shown in Figure 5. Here, each co-occurrence feature plotted in Figure 5 is a mean of all of the calculated co-occurrence features of the acquired images for the auspicious pattern. Figure 5 shows that for any auspicious pattern the co-occurrence features in the double-diagonal orientations are approximately equal. This result demonstrates that all four of the auspicious patterns have texture symmetry in the double-diagonal orientations. The co-occurrence features of the four auspicious patterns are discussed as follows:

**Figure 5.** Texture features of auspicious pattern examples in different orientations: (**a**) ASM. (**b**) CON. (**c**) HOM. (**d**) ENT. (**e**) DIS. (**f**) COR. (**g**) CLU.

In ASM the fish scales and turtle back patterns obtained the highest and lowest cooccurrence feature values, respectively. This result demonstrates that the fish scales pattern images are more constant than the others. For any auspicious pattern, the distribution of gray levels in the horizontal orientation is more uniform than those in the other orientations.

In CON, Figure 5b illustrates that the imagery contrasts of the auspicious patterns in the double-diagonal orientations are higher than those in the horizontal or vertical orientations. Except for the "回" pattern, the co-occurrence feature values in the vertical orientation are lower than those in the other orientations, where the co-occurrence feature values of the turtle back pattern are the highest. Thus, compared with the other patterns, the turtle back pattern usually has the highest imagery contrast. Figure 5b,c illustrate that the concave shapes of the obtained curves in HOM run contrary to the convex ones in CON. Due to the slightly higher HOM values of the fish scales pattern, the textural descriptor of HOM demonstrates that the fish scales pattern has a lower contrast than the other patterns. The lower imagery contrast also means a greater imagery constant.

In addition to CON and HOM, the shapes of the obtained curves between ENT and ASM are also contrary. The ENT value of the "回" pattern in the vertical orientation is clearly much larger than that of the other patterns in the different orientations. This result shows that, compared with the other patterns, the "回" pattern in the vertical orientation has the highest textural complexity. Excluding the vertical orientation, the turtle back pattern has the highest textural complexity. However, there are similar shapes of the obtained curves between CLU and COR. In COR, the values in the horizontal and vertical orientations are far higher than those in the double-diagonal orientations. Among the auspicious patterns, the "回" pattern has an extreme difference in the COR values between the horizontal or vertical orientations and the double-diagonal ones. Thus, this paper demonstrates that greater spatial dependency of gray levels in the auspicious pattern images means higher CLU values, i.e., higher textural clustering.

In DIS, Figure 5e shows that the auspicious pattern images have higher dissimilarity of the gray levels in the double-diagonal orientations than in the horizontal or vertical orientations. It is noticeable that the shapes of the obtained curves between CON and DIS are similar (see Figure 5b,e). Hence, textural analysis demonstrates that there is an interrelationship between CON and DIS for the four auspicious patterns. A higher contrast of an auspicious pattern image means a higher dissimilarity of the gray levels in the image.

#### *6.2. Similarity between Any Two Auspicious Patterns*

A two-tailed *t*-test was applied to the calculated co-occurrence features to see the similarity between any two of the four auspicious patterns, and the results of the six combinations are shown in Table 2. The results confirm that there is the most textural difference between the fish scales and lock chain patterns, and the most textural similarity between the lock chain and turtle back patterns. This result indicates that it is easy to distinguish the fish scales patterns from lock chain patterns, but very difficult to distinguish the lock chain patterns from turtle back patterns. In spite of that, the textures of the fish scales and lock chain patterns were sometimes confused when considering some of the textural descriptors in certain orientations, e.g., CON in the vertical orientation.

Table 2 also indicates that the fish scales and turtle back patterns in the horizontal and vertical orientations almost have a significant texture similarity. However, applying the textural descriptors of ASM, HOM, ENT, and DIS to the texture features in the doublediagonal orientations confirms that there is a significant difference between the fish scales and turtle back patterns. In other words, the four above textural descriptors in the doublediagonal orientations are useful in recognition of the fish scales and turtle back patterns. Finally, most of the textural descriptors in any orientation are inadequate to distinguish the "回" pattern from the others unless one considers the textural descriptor of COR.


**Table 2.** Two-tailed *t*-test results for similarity between any two auspicious patterns.

H0: supports the hypothesis that the texture features between the two auspicious patterns are similar; H1: supports the hypothesis that the texture features between the two auspicious patterns are dissimilar. Confidence level α: 0.05.

#### *6.3. Similarity of an Auspicious Pattern between Any Two Surveyed Artists*

Figure 6 shows a diagram of master and apprentice relationships for the surveyed artists, who acquired more than 10 auspicious pattern images for the fish scales or turtle back patterns. Some of the artists shown in Figure 6 were not surveyed for this paper, but they are/were the teachers of the surveyed artists, and are necessary to display in the diagram of the master and apprentice relationships. Figure 6 illustrates that the surveyed artists can be approximately separated into two family genres, i.e., Y.-F. Chen and C.-Y. Pan. Indeed, C.-C. Liu not only apprenticed to his uncle (belonging to the family genre of Pan), but also learned from S.-Y. Chen (belonging to the family genre of Chen) and the other artists by observing their treatment skills with the door god paintings. Thus, the theory is that the painting style of C.-C. Liu should involve those of the two family genres. On the other hand, F.-T. Kuo (Chinese transliteration: 郭佛賜) and Y.-C. Chen (Chinese transliteration: 陳陽春) do not have any master and apprentice relationships with the other surveyed artists, including with one another. Therefore, the theory is that there should be a great difference in the painting styles of F.-T. Ku or Y.-C. Chen from the others, including one another.

The two-tailed *t*-test results for the similarities between any two surveyed artists' fish scales and lock chain patterns in the four orientations are recorded in the symmetric matrices as shown in Tables 3 and 4, respectively. Excluding the diagonal elements of the matrices, each element uses a vector to record the two-tailed *t*-test result. In a vector of [], the elements express the acceptances of the hypotheses based on the co-occurrence features of ASM, CON, HOM, ENT, DIS, COR, and CLU in sequence, and the symbols of "0" and "1" indicate confirmation of the H0 and H1 hypotheses, respectively.

**Figure 6.** A diagram of master and apprentice relationships between the surveyed artists (black coloring) and their mentors (gray coloring).

**Table 3.** Two-tailed *t*-test results for similarity between any two surveyed artists' fish scales patterns in the orientations.


Confidence level α: 0.1. The elements in a vector of [] indicate the hypotheses supported based on the cooccurrence features of ASM, CON, HOM, ENT, DIS, COR, and CLU, in sequence, and "0" and "1" denote H0 and H1, respectively.


**Table 4.** Two-tailed *t*-test results for similarity between any two surveyed artists' lock chain patterns in the orientations.

Confidence level α: 0.1. The elements in a vector of [] indicate the hypotheses supported based on the cooccurrence features of ASM, CON, HOM, ENT, DIS, COR, and CLU, in sequence, and "0" and "1" denote H0 and H1, respectively.

Based on the diagram of the master and apprentice relationships in Figure 6, the two-tailed *t*-test results in Tables 3 and 4 are discussed. Most of the two-tailed *t*-test results successfully indicate the differences in the painting styles between the two main family genres of Chen and Pan. For instance, Table 3 demonstrates that the painting style of S.-Y. Chen is significantly different from those of L.-S. Pan, L.-C. Tsai, and Y.-H. Pan. Moreover, Table 4 shows that the painting style of L.-S. Pan is robustly different from those of S.-Y. Chen and W.-N. Chuang. However, there is an unreasonable result of recognizing the painting styles between the two main family genres of Chen and Pan as W.-N. Chuang versus T.-F. Su (see Table 4); W.-N. Chang and T.-F. Su belong to the family genres of Chen and Pan, respectively, but the two-tailed *t*-test result indicates that they have a robust similarity in painting style.

In the same family genre, this paper notes that the painting styles of the direct relatives of master and apprentice—such as L.-S. Pan versus L.-C. Tsai, T.-F. Su, or Y.-H. Pan—should be extremely similar. Unfortunately, Tables 3 and 4 cannot effectively demonstrate that a similar painting style exists among the above direct relatives of master and apprentice. For the collateral relatives of apprentices, Table 4 shows in the family genre of Chen that the painting styles of S.-Y. Chen and W.-N. Chuang have higher textural similarity in the double-diagonal orientations than in the horizontal and vertical orientations. The above characteristic also exists in the family genre of Pan, such as with C.-C. Liu versus T.-F. Su, but excludes L.-C. Tsai versus C.-C. Liu or Y.-H. Pan.

For the fish scales pattern, the co-occurrence features of ASM, HOM, ENT, and DIS have higher applicability to demonstrate the master and apprentice relationships in Figure 6 than the others. Based on the four above co-occurrence features, however, L.-C. Tsai and Y.-C. Chen, relative to the other surveyed artists in Table 3, seem to obtain fewer reasonable demonstrations of master and apprentice relationships between them and the other artists. For the lock chain pattern, except for COR, the other co-occurrence features have approximate applicability to demonstrate the master and apprentice relationships. Furthermore, among the surveyed artists in Table 4, T.-F. Su shows the fewest reasonable demonstrations of master and apprentice relationships between him and the other artists. Consequently, in this paper the three artists L.-C. Tsai, Y.-C. Chen, and T.-F. Su were excluded from the following two-tailed *t*-test to observe the similarity of the fish scales and lock chain patterns between the surveyed and unknown artists. Moreover, ASM, HOM, ENT, and DIS were considered in the two-tailed *t*-test based on the fish scales pattern, and in addition to the four above co-occurrence features, CON and CLU were considered in the test based on the lock chain pattern.

#### *6.4. Similarity of Auspicious Patterns between Surveyed and Unknown Artists*

The two-tailed *t*-test was also used to assist in the painting style recognition of the unknown artists' door god paintings in Kinmen. If a robust similarity of an auspicious pattern between the surveyed and unknown artists is obtained, this suggests that the door god paintings of the unknown artist(s) could be restored by referring to the painting style of the surveyed artist. Tables 5 and 6 show the two-tailed *t*-test results for the similarities of the fish scales and lock chain patterns between the surveyed and unknown artists, respectively, and demonstrate that F.-T. Kuo and W.-N. Chuang have the most similar painting styles to the unknown artists' fish scales and lock chain patterns, respectively.

To further observe the painting styles of each surveyed artist, we performed an interquartile range calculation for the co-occurrence features extracted from the acquired auspicious pattern images. Therefore, one co-occurrence feature in the four orientations can calculate four medians via the interquartile range analysis, and this paper adopts a mean value of the four medians of the co-occurrence feature to describe the painting style. The charts shown as Figure 7 plot the means of the medians of the extracted co-occurrence feature values for the surveyed artists' fish scales and lock chain patterns, respectively. By discussing Figure 7, the painting styles of the fish scales and lock chain patterns of the unknown artists in Kinmen could be surmised.


**Table 5.** Two-tailed *t*-test results for similarity between surveyed and unknown artists' fish scales patterns in the orientations.

Confidence level α: 0.1. The elements in a vector of [] indicate the hypotheses supported based on the co-occurrence features of ASM, HOM, ENT, and DIS, in sequence, and "0" and "1" denote H0 and H1, respectively.

**Table 6.** Two-tailed *t*-test results for similarity between surveyed and unknown artists' lock chain patterns in the orientations.


Confidence level α: 0.1. The elements in a vector of [] indicate the hypotheses supported based on the cooccurrence features of ASM, CON, HOM, ENT, DIS, and CLU, in sequence, and "0" and "1" denote H0 and H1, respectively.

**Figure 7.** Means of medians of co-occurrence feature values in four orientations extracted from acquired auspicious pattern images of surveyed artists: fish scales pattern (**upper chart**); lock chain pattern (**lower chart**).

The upper chart of Figure 7 shows that F.-T. Kuo's fish scales patterns have lower ASM and HOM and higher ENT and DIS than those of the other surveyed artists. The lower ASM and HOM mean that the acquired fish scales pattern images of F.-T. Kuo exhibit lower images constancy, leading to a higher image contrast. Based on the higher ENT and DIS, this paper finds that F.-T. Kuo's fish scales patterns (or the unknown artists' fish scales patterns in Kinmen) have higher textural complexity and higher hue gradient than those of the other surveyed artists.

In addition to ASM, HOM, ENT, and DIS, CON and CLU can help describe the painting style of W.-N. Chuang's lock chain patterns (or the unknown artists' lock chain patterns in Kinmen). The lower chart of Figure 7 shows that W.-N. Chuang's lock chain patterns also have lower ASM and HOM and higher ENT and DIS than those of the other surveyed artists. Because the representation of HOM is contrary to that of CON, the lower chart illustrates that W.-N. Chuang's lock chain patterns have higher CON than those of the other surveyed artists. The lower chart also shows that W.-N. Chuang's lock chain patterns have the highest CLU among the surveyed artists. Thus, this result indirectly explains why the unknown artists' lock chain patterns in Kinmen have higher textural clustering than the surveyed artists.

From the above discussion, this paper finds that there should be a similar door god painting style between F.-T. Kuo and W.-N. Chuang, because both of their auspicious patterns have the same characteristics of lower ASM and HOM and higher ENT and DIS. Thus, the door god paintings in Kinmen could be restored by referring to the painting styles of F.-T. Kuo or W.-N. Chuang.

#### **7. Conclusions**

Based on a two-level two-dimensional DWT and gray-level co-occurrence matrix calculation, this research proposes a texture feature extraction procedure and, coupled with a two-tailed *t*-test to analyze the door god painting styles, demonstrates the master and apprentice relationships between the surveyed artists so as to recognize unknown artists' painting styles. In total, 52 temples in Taiwan were visited to acquire 453 images for the auspicious patterns—including fish scales, lock chain, turtle back, and "回"—on the armor of door god warriors created by the surveyed artists. Additionally, 12 auspicious pattern images of unknown artists' door god paintings in Kinmen, which is an outlying island of Taiwan, were also acquired.

In the door god painting style analysis, the result indicates that all four auspicious patterns have texture symmetry in the double-diagonal orientations. The texture of the fish scales patterns is more constant (i.e., lower contrast) than that of the other patterns. For any auspicious pattern, the textural contrast in the double-diagonal orientations is higher than those in the horizontal or vertical orientations and, among the four auspicious patterns, the turtle back pattern usually has the highest textural contrast. Except in the vertical orientation, the turtle back pattern has the highest textural complexity. In the vertical orientation, the highest textural complexity exists in the "回" pattern. Generally, there is the most textural difference between the fish scales and lock chain patterns, and the most textural similarity between the lock chain and turtle back patterns. The analysis herein also successfully demonstrates the interrelationships between the seven co-occurrence features.

The results of our analysis effectively demonstrate the differences in the painting styles between the family genres of Chen and Pan. In a family genre, however, the results of our analysis fail to indicate that artists with a direct relationship of master and apprentice are expected to have a similar painting style. As for the collateral relatives of apprentices, our findings confirm that the artists' auspicious patterns seem to have higher textural similarity in the double-diagonal orientations than in the horizontal and vertical orientations. Through the two-tailed *t*-test, based on the fish scales pattern, we found that the co-occurrence features of ASM, HOM, ENT, and DIS are useful to demonstrate the apprentice relationships between the five artists of C.-C. Liu et al. Based on the lock chain pattern, except for COR, the other six co-occurrence features are useful to demonstrate the apprentice relationships between the four artists of C.-C. Liu et al.

Among the surveyed artists, F.-T. Kuo and W.-N. Chuang have the most similar painting styles to the unknown artists' fish scales and lock chain patterns, respectively. Thus, the door god paintings in Kinmen could be restored by referring to the painting styles of F.-T. Kuo or W.-N. Chuang. According to the analyzed painting styles of Kuo and Chuang, this paper surmises the painting styles of unknown artists in Kinmen, indicating that the unknown artists' fish scales patterns exhibit higher textural contrast, textural complexity, and hue gradient than those of the surveyed artists. Furthermore, the unknown artists' lock chain patterns also have higher textural contrast and clustering than those of the surveyed artists.

The hypothesis testing method can currently determine the useful co-occurrence features for recognizing the unknown artists' painting styles. In the future, the useful co-occurrence features will be input into neural networks in order to recognize the style of every single door god painting. We also believe that deep learning convolutional neural networks based on the human visual system would be useful in recognizing the style of every single door god painting.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //drive.google.com/file/d/1ax76s22qelbdEmsi5ILhCKe3HlrvhAoP/view?usp=sharing (accessed on 1 December 2021).

**Author Contributions:** Conceptualization, T.-C.S. and T.-C.W.; methodology, T.-C.S.; software, T.-C.S. and M.-H.W.; validation, T.-C.S. and T.-C.W.; formal analysis, M.-H.W. and C.-W.W.; investigation, T.-C.S., M.-H.W. and C.-W.W.; data curation, M.-H.W.; writing—original draft preparation, T.-C.S.; writing—review and editing, T.-C.S. and T.-C.W.; project administration, T.-C.S.; funding acquisition, T.-C.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Technology, Taiwan, grant number MOST 108-2410-H-507-013-MY2. The APC was funded by the Ministry of Science and Technology, Taiwan.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Haoyang Jiao 1,2,3, Fayuan Li 1,2,3,\*, Hong Wei 1,2,3 and Wei Liu 1,2,3**


**Abstract:** Shoulder lines can best depict the morphological characteristics of the Loess Plateau. Moreover, a shoulder line depicts the external appearance of spatial differentiation of loess landforms and the internal mechanism of loess landform evolution. The efficient and accurate extraction of shoulder lines can help to deepen the re-understanding of the morphological structure and differentiation of loess landforms. However, the problem of shoulder line continuity in the extraction process has not been effectively solved. Therefore, based on high-resolution satellite images and digital elevation model (DEM) data, this study introduced the regional growing algorithm to further correct edge detection results, thereby achieving complementary advantages and improving the accuracy and continuity of shoulder line extraction. First, based on satellite images, the edge detection method was used to extract the original shoulder lines. Subsequently, by introducing the regional growing algorithm, the peaks and the outlet point extracted with the DEM were used as the growth points of the positive and negative (P-N) terrains to grow in four-neighborhood fields until they reached a P-N terrain boundary or a slope threshold. Finally, the P-N terrains extracted by the regional growing method were used to correct the edge detection results, and the "burr" was removed using a morphological image-processing method to obtain the shoulder lines. The experimental results showed that the method proposed in this paper can accurately and effectively complete the extraction of shoulder lines. Furthermore, the applicability of this method is better and opens new ideas for quantitative research on loess landforms.

**Keywords:** satellite images; shoulder line; edge detection; regional growing algorithm; positive and negative (P-N) terrain

#### **1. Introduction**

One of the most significant recent discussions has been the study of spatial differentiation in the Loess Plateau [1–5]. Shoulder lines are the topographic structural lines that best depict loess landform characteristics. The shape, grade, spatial distribution, development trend, and other characteristics of shoulder lines reflect the regional variations in loess landforms, as shown in Figure 1. However, due to the complexity of the Loess Plateau, the effective extraction of shoulder lines is influenced by a variety of factors, such as the accuracy of the source data, the local discontinuity of the lines, and the applicability of the method [6–10]. Consequently, a more efficient and accurate shoulder line extraction method is still needed. The traditional method of extracting shoulder lines is to use the contour lines of topographic maps or aerial images to directly delineate them [11–15]. Although the precision is excellent, the efficiency of this method is low, and the workload is heavy.

**Citation:** Jiao, H.; Li, F.; Wei, H.; Liu, W. An Improved Shoulder Line Extraction Method Fusing Edge Detection and Regional Growing Algorithm. *Appl. Sci.* **2022**, *12*, 12662. https://doi.org/10.3390/ app122412662

Academic Editor: Tung-Ching Su

Received: 26 October 2022 Accepted: 7 December 2022 Published: 10 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** An illustration of shoulder line of the Loess Plateau.

In recent years, with the support of RS and GIS technologies, the characteristic terrain elements in a landform can be automatically extracted and analyzed, and the extraction of shoulder lines has become one of the hot spots for study [16–21]. Researchers have studied the automatic extraction of shoulder lines. The extraction methods can be classified into three types: The first method is the geometric morphological method. Terrain parameters, such as slope, slope aspect, and curvature, are used for regular judgment to obtain shoulder lines based on the given rule of shoulder lines in the geometric form. For example, Evans used DEM data as the basis for extracting shoulder lines by combining areas of low difference between mean elevation and high positive plan curvature [22]. Lei obtained shoulder lines by comparing the slopes of upstream and downstream grids on the same slopes [23]. Xiao proposed a method for extracting shoulder lines based on slope aspect orientation according to the slope-turning characteristics of loess landforms [24]. This method is simpler, faster, and easier to implement compared with the traditional approach. However, the actual effect obtained is not ideal. The second method is the hydro-geomorphologic method. The constraint and drawing of shoulder lines are achieved as a result of the interaction between the hydrological process and the landform and are based on the principle of a flow analysis of the terrain surface. For example, Lv applied a surface-flow digital simulation method and a contour vertical-line-tracking method to achieve shoulder line extraction based on the principle of a terrain surface flow analysis [25]. Liu used water flow paths for distributed water flow computation to extract shoulder lines [26]. Yang used the direction of confluence to determine the streamline of a slope and drew a shoulder line based on the inflection point of the streamline [27]. The shoulder lines obtained by this method have a certain geological significance compared with the geometric morphological method, and the problem of shoulder line continuity is overcome to a certain extent. However, the digital simulation of surface water is influenced by several factors, including water flow direction and catchment threshold, which increase the uncertainty of a slope's flow path, making judging points on a characteristic shoulder line difficult, lacking in precision, and necessitating substantial calculations. The third method is based on an image-processing method. The edge detection method achieves the purpose of detecting sudden edge change by comparing the changing characteristics in the image brightness value and detecting the sudden change point of the image brightness. For example, Vrieling used the supervised classification approach for the maximum likelihood classifier to classify gullies and non-gullies [28]. Yan used image binarization and various edge detection operators to extract shoulder lines [29]. Wang extracted a shoulder line by combining the P-N opening of the terrain and the threshold segmentation of a difference image [30]. Compared to the previous two methods, edge detection operators can quickly and effectively extract such changing features, and some operators can even extract weak mutation features to better

reflect details; this provides a method for indistinct shoulder line extraction. However, the shoulder lines extracted with edge detection methods tend to be less closed, and the line segments are fragmented and do not match the actual shoulder lines. In conclusion, the above methods share the same problem in that they require a means to effectively balance continuity and accuracy. Drawing on the advantages of the good applicability of the edge detection method, finding ways to further improve the integrity of shoulder lines is an urgent problem to be solved. The regional growing algorithm, a fundamental feature in image segmentation research, can produce continuous, closed regions with relatively little domain information [31,32]. The regional growing method has advantages for solving the continuity problem of shoulder lines based on DEM data [33]. However, the resolution of DEM data has a great influence on the extraction results, and it is difficult to obtain highprecision DEM data, especially in large areas. Hence, a shoulder line extraction method that combines the edge detection method and regional growing algorithm is proposed here to comprehensively improve the continuity and accuracy of shoulder lines. Based on highresolution satellite images, this study uses the advantages of edge detection in highlighting details to complete the extraction of an original shoulder line. At the same time, based on DEM data, the regional growing algorithm is used to realize the determination of P-N terrain, to overcome the discontinuities of the original shoulder line, and to complete the accurate extraction of the shoulder line.

#### **2. Materials and Methods**

#### *2.1. Study Areas*

The study area is one of the key national governance regions for soil and water conservation and is located in Yijun and Luochuan counties, Shaanxi Province, China, with latitudes between 35◦25 N and 34◦42 N and longitudes between 109◦20 E and 109◦37 E, as shown in Figure 2. The topography of this place is the Loess Plateau. In this area, the elevation ranges from 85 m to 3719 m. The watershed unit is usually used as the basic landform analysis unit due to its relative internal homogeneity; hence, in this study, the extraction of shoulder lines was also presented with small watersheds as the sample area.

**Figure 2.** Location of the study area: (**a**,**b**) location of Loess Plateau; (**c**,**c1**–**c4**) hillshade map of the study area.

#### *2.2. Data*

In this study, high-resolution satellite images were used to produce higher-density spectral and textural information. Imagery of3mresolution downloaded from Planet Explorer (https://www.planet.com/explorer/, (accessed on 1 September 2021)) was used. The image data were mainly captured between June and September 2021. During this period, the images had less cloud content, and crops were harvested on the land surface; therefore, the boundaries of objects with loess morphologies were more obvious, which aided image edge detection and shoulder line extraction. The Advanced Land-Observing Satellite (ALOS) Digital Elevation Model with a spatial resolution of 12.5 m was used to calculate the peak point and outlet point of the sample area. The DEM data were downloaded from the National Aeronautics and Space Administration (NASA, https://search.asf.alaska.edu/#/, (accessed on 1 October 2016)). The horizontal and vertical accuracies of the elevation data could reach 12.5 m. These were the most accurate data from open-source DEM data. More details on the data are shown in Table 1.

**Table 1.** Details of data used in this study.


#### *2.3. Edge Detection*

In image-processing technology, commonly used edge detection operators include the Roberts operator, Sobel operator, Prewitt operator, Laplace operator [34], and Canny operator [35]. These operators have different edge detection capabilities for images with different characteristics. In an image, shallow gullies appear as linear features with weak feature information. The Canny operator has the characteristics of good anti-interference ability and accurate positioning and can effectively identify and locate edges on slope data. Compared with other edge detection operators, the Canny operator can seek the best solutions for anti-noise solutions and precise positioning. Therefore, this study selected the Canny operator to complete detection of the experimental sample area. To make the shoulder line features more obvious and to avoid image noise interference, the images needed to be preprocessed.

#### 2.3.1. Image Grayscale

Performing grayscale operations on images can reduce the amount of computation needed. The maximum value of three image components is taken as the result of image grayscale processing, as shown in Formula (1):

$$\text{Gray}(i, j) = \max[R(i, j), G(i, j), B(i, j)]\tag{1}$$

where *Gray*(*i*, *j*) represents the gray value of an image, and *R*(*i*, *j*), *G*(*i*, *j*), *B*(*i*, *j*) represent the three components of the image.

#### 2.3.2. Binary Image

A target image has a large difference in gray value from its background image, and it is partitioned according to the gray value. The gray value of a target image is marked as 0, and the gray value of a background image is marked as 1. If *F*(*x*, *y*) is the gray value of a pixel in the image, the transformation function of the gray threshold *Th* is as follows:

$$F(\mathbf{x}, \mathbf{y}) = \begin{cases} 1, F(\mathbf{x}, \mathbf{y}) > Th \\ 0, F(\mathbf{x}, \mathbf{y}) \le Th \end{cases} \tag{2}$$

#### 2.3.3. Canny Edge Detection

To perform edge detection on the preprocessed target images, the Canny edge detection method included the following four steps: The first step was to smooth the image. We constructed a filter with a Gaussian function, performed convolution filtering on the images, removed noise, and obtained smooth images. The second step was to calculate the gradient magnitude and gradient direction of the images. The gradient magnitude and gradient direction of the smoothed images were calculated by the finite difference in the first-order partial derivatives. The third step was to perform non-maximum suppression on the gradient amplitude. To determine the edge, it was necessary to keep the point with the largest local gradient and to suppress the non-maximum value, that is, we set the non-local maximum value point to zero in order to obtain a refined edge. If the gradient value of the neighborhood center point was larger than the value of the two adjacent points along the gradient direction, the current neighborhood center point was determined as a possible edge point. Otherwise, it was assigned a value of zero, and the pixel point was judged as a non-edge point. The fourth step was double threshold detection. After applying nonmaximal suppression, the remaining edge pixels provided a more accurate representation of the true edges in the images. However, some edge pixels were still caused by noise and color variations. To account for these spurious responses, edge pixels must be filtered out with weak gradient values, as well as edge pixels with high gradient values. This was achieved by choosing high and low thresholds. If the gradient value of an edge pixel was higher than the high threshold, it was marked as a strong edge pixel. If the gradient value of an edge pixel was less than the high threshold and greater than the low threshold, it was marked as a weak edge pixel. If the gradient value of an edge pixel was less than the low threshold, it was suppressed.

#### *2.4. Regional Growing Algorithm*

The regional growing algorithm is a process of merging pixels or sub-regions into a larger area according to similarity criteria [36]. This algorithm is based on the theory that regions start with a group of growing points and that the same or similar adjacent pixels merge into new growing points. This process is continuously repeated until there are no more points to merge. The three steps are as follows: (1) choosing the appropriate growing points, (2) determining the similarity criteria, and (3) establishing the stop rules.

#### 2.4.1. Identifying Growing Points for P-N Terrain

Shoulder lines lie on the boundaries of P-N terrain. Accordingly, extracting the P-N terrain boundary is considered a premise of shoulder line extraction. Positive terrain is an area that is higher than the adjacent region or that is located in a tectonic uplift region. The P-N terrain method can be used to classify positive terrain. Errors in the positive terrain can be classified into two categories, as follows: (1) depressions, where the condition is caused by artificial modification or slight topographic relief in small regions, and (2) flatland, where a small difference can be observed between the original elevation and the elevation after smoothing when using a filter window slide on a nearly flat DEM. Slight elevation changes can affect the results. This phenomenon is especially evident in the Loess Plateau area, which means that peaks located in positive terrain must be correctly classified. Hence, peaks should be chosen as the growing points of positive terrain and should grow until there are no more points of the same type to merge. The above-mentioned positive terrain was extracted using the P-N terrain method. The size of the analysis window depended on the fragmentation of the landform. If the landform had more fragments, then we tended to choose a smaller analysis window. Negative terrain is an area that is lower than the adjacent region or an area that is located at a tectonic down-lift region. The test area was a complete watershed, and the negative terrain was connected. Accordingly, one growing point for negative terrain was enough. The outlet was the lowest elevation point and was located in negative terrain. On this basis, the outlet point was chosen as the growing point

of the negative terrain and grew until there were no more points of the same type to merge. The above-mentioned negative terrain was extracted using the P-N terrain method.

#### 2.4.2. Growth Criteria

We assumed that the pixel with a gray value of 8 was the initial growth point, as shown in Figure 1, where the numbers in the figure indicate the gray levels of pixels. The growing criterion for the four neighborhoods of 8 was that the calculated result should be in the range between −1 and 1 when the growing point was less than the neighbor point. If no difference could be found between the adjacent growing results, then the growing process stopped. Figure 3 illustrates four stages of the regional growing algorithm. Figure 3a depicts an original image, and the numbers are gray values. Figure 3b–d represents the growth process.

**Figure 3.** Illustration of regional growth: (**a**–**d**) represent the four steps of the regional growing algorithm, respectively.

The regional growing algorithm can avoid most misclassification areas. Still, for those regions that are connected to the correct area, it remains difficult to accurately identify them. This problem results in the inaccurate position of the shoulder line, and it is serious in the Loess Plateau area. To improve this concept, this study introduced slope gradient, and defined the following rules: if the slope gradient of the negative terrain was smaller than the given threshold and the negative terrain was adjacent to the positive terrain, then this negative terrain could be considered as positive terrain, as shown in Figure 4. The threshold of the slope gradient depended on the type of loess landform. This study tested different thresholds, compared the test results with the manual visual interpretation results, and chose 7◦ as the final slope threshold.

**Figure 4.** Slope's influence on the positive and negative terrain classification results: (**a**) represents the classification map of positive and negative terrain; (**b**) represents the results after slope correction.

#### *2.5. Burr Removal*

The resulting shoulder line had some parasitic components in line corners, as shown in Figure 5. These parasitic components are called burrs. These burrs can lead to uncertain positions of the shoulder line. When the shoulder line was transformed from a grid to a vector, this defect became more prominent.

**Figure 5.** Contrast between before (**a**) and after (**b**) burr removal.

This study used a morphological image-processing method to remove burrs. The basic principle of this method was that fixed structure elements were used to detect and thin the endpoint (Equation (3)). Structure elements were used to detect and remove burrs. The algorithm was as follows: (1) The threshold of endpoint thinning was determined. We observed the image; if the length of a burr was less than 3 pixels, the threshold was set to 5 pixels to guarantee the reliability of the experimental results. (2) The structure elements were determined. The diagonal pixels could be ignored because the growing method only used four neighbor pixels. This study used the structure elements as shown in Figure 6. In these elements, "×" represents ignored pixels, "1" depicts shoulder line pixels, and "0" denotes background pixels. In an analysis window, if the image value agreed with any structure element in Figure 6, then "0" was set as the center pixel value (Equation (4)). The whole image was scanned using structure elements that could only remove one endpoint pixel at once. The above operation was repeated five times to ensure that the burrs were completely removed. The difference between before and after the above-mentioned steps is shown in Figure 5. (3) The correct end points were recovered. Given that this method processed all the end points, the correct end points were also removed, as were the burrs. The next step was to recover the correct end points. Structure elements were used to scan the whole image and remove the end points again. The burrs were already removed after the previous five scans. Accordingly, the pixels removed at this time were the correct end points. Thereafter, we performed four neighbor expansions and obtained the intersection with the original shoulder line (Equation (5)). The expansion time was the same as the thinning time. The complete process of this study method is shown in Figure 7.

$$X\_1 = A \odot \{B\} \tag{3}$$

where *A* is the original shoulder line, *B* is the structure elements, is the thinning, and *X*<sup>1</sup> is the thinning result.

$$X\_2 = \bigcup\_{k=1}^{8} \left( X\_1 \odot B^k \right) \tag{4}$$

where *<sup>X</sup>*<sup>2</sup> is the end point, *<sup>B</sup><sup>k</sup>* (*<sup>k</sup>* = 1, 2, 3, 4) are the structure elements, and is the hit-miss transformation.

$$X\_3 = (X\_2 \oplus H) \cap A \tag{5}$$

where *X*<sup>3</sup> is the shoulder line, *H* is a four-neighborhood structure element, *A* is the original shoulder line, and ⊕ is the expansion.



**Figure 6.** Structure elements.

**Figure 7.** The complete process of our experiment.

#### *2.6. Accuracy Assessments*

To objectively evaluate and verify the superiority of the method in this paper, the following indicators were used to evaluate the edge detection method, the regional growing algorithm, and the method in this paper. Class pixel accuracy (CPA) is the percentage of correct pixels out of all the extracted result pixels; the closer the value is to 100%, the better. Pixel accuracy (PA) is the percentage of correctly extracted pixels in an image, that is, the proportion of correctly extracted pixels out of the total pixels; the larger the value, the better. The dice similarity coefficient (Dice) is a measure of set similarity, indicating the ratio of the area where two objects intersect the total area; the value range is (0, 1), and the effect is best when it is 1. Intersection over union (IOU) is the overlapping area between the extraction results and the real value divided by the joint area between the extraction results and the real value; the value is between (0, 1), and the larger the value, the better the effect. The formulas for each indicator are as follows:

$$\text{CPA} = \frac{TP}{\left(TP + FP\right)} \times 100\% \tag{6}$$

$$PA = \frac{(TP + TN)}{(TP + TN + FP + FN)} \times 100\% \tag{7}$$

$$\text{Dice}(A, B) = \frac{2|A \cap B|}{|A| + |B|} \tag{8}$$

$$IOU = \frac{|A \cap B|}{|A \cup B|} \tag{9}$$

where *A* is the method extraction result; *B* is the manual visual interpretation result; *TP* = (*A* ∩ *B*) is the method correctly extracting the region; *FN* = *A* − (*A* ∩ *B*) is the target area missed by the method; *TN* = *I* − *A* is the real background area; and *I* is the set of image pixels.

#### **3. Results**

#### *3.1. Parameter Settings*

When using edge detection to detect high-precision remote-sensing images, to identify the edges of gullies more accurately, it was necessary to set appropriate parameters for the Canny operator function. The Canny function had two parameters: The first parameter represented the first threshold, and the calculated boundary points were greater than this threshold to be the real boundary. The second parameter represented the second threshold; the calculated boundary points below this threshold were discarded. Based on the understanding of the terrain features of the study sample area, the parameters were adjusted, and the final optimal first and second parameters were 50 and 210, respectively, when using the regional growing algorithm to generate positive and negative terrain based on DEM data. Based on the comparison of extraction results that used different window sizes, the window size of 13 × 13 was the best size for shoulder line extraction [37]. Positive and negative terrain growth points were set as peaks and outlet points, respectively.

#### *3.2. Results Analysis*

In this study, the Canny operator was used to complete the detection of the gully range on high-resolution satellite images, and the extraction results were superimposed on the image data, as shown in Figure 8. Based on the DEM data, the regional growing algorithm was used to complete the extraction of positive and negative terrain for each sample area. This study chose peaks and outlets as P-N terrain growing points, respectively. The peak and outlet points could be detected through neighborhood analysis and the watershed boundary method. In the process of growing, the same values of pixels were merged until there were no points to merge. The generated results were also superimposed on hillshade data, as shown in Figure 9. Overall, compared with high-resolution satellite images and hillshade data, the extraction results of edge detection and the regional growing algorithm were better.

**Figure 8.** The result of extracting Negative terrain by edge detection method. (**a**–**d**) represent 4 sample areas.

**Figure 9.** The results of extracting P-N terrain using the regional growing algorithm: (**a**–**d**) represent 4 sample areas.

To compare the performance of the regional growing and edge detection methods for the extraction of the shoulder line, the real shoulder line needed to be defined. In this paper, manual visual interpretation supported by expert knowledge was employed to extract the real shoulder line. The shoulder line extracted using the manual visual interpretation method based on high-resolution satellite images was used as the evaluation criterion. The extraction results of the regional growing and edge detection methods were evaluated by comparing the negative terrain areas of these two methods with the real shoulder line results. It can be seen from Table 2 that, for the four sample areas, the regional growing algorithm was generally better than the edge detection method when the negative terrain

areas were used as the evaluation indicator. For example, sample area a, the difference between the results obtained by the regional growing algorithm and the standard results was within 12 ha, and the error was within 5%; the difference between the results obtained by the edge detection method and the standard results was within 41 ha, and the error was within 13%. These results also showed that the results of shoulder line extraction are credible when studying the characteristic indicators related to area, such as gully erosion. However, in the details, we can find that there were still differences between the two methods. The specific analysis was as follows. It can be found from Figure 10 that the regional growing algorithm was more advantageous in the detailed expression of the shoulder line, such as the turning point of the shoulder line. However, it can be found from Figure 11 that, although the edge detection method could reflect more details, the extraction results were fragmented, resulting in inaccurate positioning of the shoulder line in some places, and the generated shoulder line was discontinuous. Although DEM data do not contain more information than high-resolution satellite images, they become an advantage for the positive and negative terrain generated by the regional growing method. This is because the core of the regional growing algorithm is the determination of the growing point. We could accurately find the lowest point (the outlet) and the highest point (the peak) as the positive and negative topographic growth points based on the DEM data using the digital terrain analysis method. Furthermore, for shoulder line continuity obtained through extraction, the regional growing algorithm had more advantages.

**Table 2.** Comparison of extracted negative terrain area with standard area differences.


**Figure 10.** Results of edge detection compared to the regional growing algorithm: (**a**–**d**) represent 4 sample areas; (**a1**,**b1**,**c1**,**d1**) represent the results of the edge detection method; (**a2**,**b2**,**c2**,**d2**) represent the results of the regional growing algorithm.

**Figure 11.** Results of the regional growing algorithm compared to edge detection: (**a**–**d**) represent 4 sample areas; (**a1**,**b1**,**c1**,**d1**) represent the results of the edge detection method; (**a2**,**b2**,**c2**,**d2**) represent the results of the regional growing algorithm.

By comparing and analyzing the shoulder line obtained by manual visual interpretation of the shoulder line based on high-resolution satellite images shown in Figure 12, we can see that the shoulder line extracted by edge detection under positive and negative terrain constraints was closer to the artificial shoulder line. Therefore, the advantages of edge detection in detail and the advantages of the regional growing algorithm in continuity were combined. By using the P-N terrain constrained edge detection results obtained by the regional growing method and removing the burrs, we could finally obtain better detailed and continuous shoulder line extraction results, as shown in Figure 13.

**Figure 12.** Comparison of shoulder lines between edge detection, regional growing algorithm, and manual visual interpretation. (**a**) represents the results of edge detection method and regional growing algorithm; (**b**) represents the results of manual visual interpretation; (**c**) represents the overlay of extraction results from different methods.

**Figure 13.** Results of improved shoulder line: (**a**–**d**) represent 4 sample areas.

#### *3.3. Precision Evaluation*

By comparing the CPA, PA, Dice, and IOU values, it can be seen from Table 3 that the extraction method in this study was better than edge detection and the regional growing algorithm. For example, in sample area a, the CPA value of the proposed method was 7.1% and 5.6% higher than those of the edge detection method and the regional growing algorithm, respectively. The largest PA value obtained was for the method proposed in this study, which was closer to 1, and the effect was better than those of the other methods. Similarly, both the Dice and IOU values were closer to 1, which was better than the other methods. Therefore, the reliability and accuracy of the method proposed in this study were verified.



#### **4. Discussion**

*4.1. Comparison of Different Operators*

For edge detection operators for different images, there are differences in edge detection. It is the premise of accurate shoulder line extraction to select an operator that is suitable for the study sample area from among many operators. Therefore, we chose sample area a to discuss the effects of common operators on the shoulder line extraction results, and the experimental results are shown in Figure 14. In this study, the line-related parameters of the shoulder line were used to further compare the detection results of each edge detection operator. These parameters were mainly used to describe the fineness of the line segment extracted by the operator. As shown in Table 4, the overall effect of the line segment extracted by the Canny operator was better than the other operators. We found that the Prewitt operator and the Sobel operator both performed differential and filtering operations on the image and only had some differences in the selection of weights for smoothing. However, the image was blurred to a certain extent, and some edges could not be detected. Therefore, the detection accuracy was relatively low, and this type of operator was deemed more suitable for situations where the gray value of an image edge is relatively obvious. The detection accuracy of the Roberts operator was relatively high, but it was easy to lose part of the edge, which made the detection results incomplete. At the same time, the image was not smoothed, and the noise could not be suppressed; thus, this operator had the best response to images with steep, low noise. The Laplace operator smoothed the image through the Gaussian function and had a relatively obvious effect on noise suppression. However, the original edge could also be smoothed during processing, resulting in some edges that could not be detected. In addition, the noise had a great influence on it; the detected image details were very rich, but at the same time, false edges could appear. If the false edges were reduced, the detection accuracy could also be reduced; many true edges could be lost, and different parameters should be selected for different images. The Canny operator was more accurate than the Laplace operator in detecting edges. Although some edge information could be lost, this operator had the best effect among the above-mentioned edge detection operators and could identify small edges more clearly.

**Figure 14.** (**a**) represents sample area a; (**b**–**f**) represent the extraction results of the Roberts, Prewitt, Sobel, Laplace, and Canny operators, respectively.


**Table 4.** Comparison of evaluation index values of different edge detection operators.

#### *4.2. Applications and Future Research*

The landform types of the Loess Plateau in northern Shaanxi show significant regional differences, and the dominant factors of landform influence vary for different landform types. Loess landform types are mainly divided into tableland, ridges, and hills [38–40], and there were challenges in realizing the fully automatic extraction of the Loess Plateau shoulder line. For the Loess Plateau, it can be seen from the previous accuracy evaluation results that the proposed method had good applicability. To further verify the applicability of this method in extracting shoulder lines, we selected a loess tableland area with an area of 35,975.531 ha. The results are shown in Figure 15. By superimposing the extracted results onto high-precision images, it can be seen that the extracted shoulder lines were

effective and reasonable, which verified that this method had good universality in Loess tableland area.

**Figure 15.** Shoulder lines extracted by applying the method proposed in this study.

In loess ridge and hill regions, the terrain was more complex. As shown in Figure 16, there were many seriously discontinuous shoulder lines due to occasional gravity erosion factors, such as landslides and scattering. Some areas were also affected by artificial landforms, such as terraced fields, dams, etc., resulting in the existence of multilevel shoulder lines (Figure 16a). At the same time, in loess ridge and hill regions, due to the influence of vegetation, the slope of the shoulder line's up- and down-slopes showed little change, and there were invisible shoulder lines (Figure 16b,c). Even if there is no interference from vegetation, in autumn and winter the continuity and visibility of shoulder lines are poor due to the impact of gravity, such as landslides, strays, and collapses, or even human factors (Figure 16d). Therefore, the application of this method in loess ridge and hill areas is also a problem that needs to be discussed and solved in the future. In future research, it is necessary to further analyze the applicability of this method for different data sources and different landform types. To better apply the method to different landform types over a large area, we could try to divide different areas according to the different land use characteristics (soil erosion characteristics and land use directions) of each landform type. On this basis, we could establish quantitative models of different shoulder lines to achieve the high-efficiency and high-precision extraction of shoulder lines.

**Figure 16.** Examples of landforms in loess ridge and hill regions. (**a**–**d**) show the distribution of shoulder lines in loess ridge and hill regions, respectively

#### **5. Conclusions**

Aiming to improve the poor continuity and inaccuracy of extracted shoulder lines, this study proposed an extraction method that fused edge detection and the regional growing algorithm. The experimental results showed that the CPA and PA values of the edge detection method were in the ranges of 77.6~83.7% and 79.1~80.1%, respectively; the CPA and PA values of the regional growing algorithm were in the ranges of 81.1~83.1% and 78.1~84.5%, respectively; the CPA and PA values of the method proposed in this study were in the ranges of 86.7~89.7% and 85.8~91.7%, respectively. Moreover, the Dice and IOU values of the method studied in this paper were closer to 1 than those of the edge detection method and the regional growing algorithm. This method could guarantee shoulder line continuity and integrity. Meanwhile, burr removal reduced errors when the grid shoulder line was transformed into a vector.

Shoulder lines have obvious turning points above and below the line, and the terrain factors (slope, curvature, etc.) also change accordingly. The geomorphic mechanisms of the P-N terrain above and below the line are significantly different. The positive terrain above the line basically maintains the original slope state after loess accumulation, and slope erosion is mainly surface erosion. The negative terrain below the line is dominated by gully erosion and gravity erosion, and various gravity landforms are widely developed. In summary, shoulder lines can be used as an important topographic index for regional soil erosion intensity and landform division. The accurate extraction of shoulder lines can provide a

new perspective for the quantitative study of loess landforms and is very important for the study of Loess Plateau landforms, soil erosion characteristics, and ecological environments.

**Author Contributions:** Conceptualization, F.L. and H.J.; methodology, H.J.; validation, H.J., H.W. and W.L.; writing—original draft preparation, H.J., W.L. and H.W.; writing—review and editing, H.J., H.W. and W.L.; supervision, F.L.; funding acquisition, F.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (grant number: 42271421).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data will be made available on reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Development of an Underground Tunnels Detection Algorithm for Electrical Resistivity Tomography Based on Deep Learning**

**Yin-Chun Hung 1,\*, Yu-Xiang Zhao <sup>2</sup> and Wei-Chen Hung <sup>1</sup>**


**Featured Application: Authors are encouraged to provide a concise description of the specific application or a potential application of the work. This section is not mandatory.**

**Abstract:** Kinmen Island was in a state of combat readiness during the 1950s–1980s. It opened for tourism in 1992, when all troops withdrew from the island. Most military installations, such as bunkers, anti airborne piles, and underground tunnels, became deserted and disordered. The entries to numerous underground bunkers are closed or covered with weeds, creating dangerous spaces on the island. This study evaluates the feasibility of using Electrical Resistivity Tomography (ERT) to detect and discuss the location, size, and depth of underground tunnels. In order to discuss the reliability of the 2D-ERT result, this study built a numerical model to validate the correctness of in situ measured data. In addition, this study employed the artificial intelligence deep learning technique for reprocessing and predicting the ERT image and discussed using an artificial intelligence deep learning algorithm to enhance the image resolution and interpretation. A total of three 2D-ERT survey lines were implemented in this study. The results indicate that the three survey lines clearly show the tunnel location and shape. The numerical simulation results also indicate that using 2D-ERT to survey underground tunnels is highly feasible. Moreover, according to a series of studies in Multilayer Perceptron of deep learning, using deep learning can clearly show the tunnel location and path and effectively enhance the interpretation ability and resolution for 2D-ERT measurement results.

**Keywords:** Electrical Resistivity Tomography (ERT); deep learning; underground tunnel

#### **1. Introduction**

Kinmen is a small island. During the 43 years of military control, the Kinmen government constructed various defence works and military camps, as well as many spiritual landmarks. Military installations, bunkers, and tunnels can be seen across the island. During the military administration, there were at least 1000 barracks, 22 large-scale monuments, 28 memorial pavilions, and 10 large underground Halls in the Kinmen area [1]. As it was in a state of combat readiness for an extended period, 120,000 soldiers were stationed in Kinmen. After the 1990s, the government regained local autonomy and democratic governance, and the original historical sites became the resources of local tourism development. In recent years, with the reduction of military garrison and opening of military spaces, the pace of development has accelerated. It is estimated that over 1000 vacant barracks will be released over the coming years. At least 300 barracks have been or are in the process of being released, and about 50 to 80 barracks are planned for annual release.

Although the vacant barracks are being released, the numerous underground bunkers are still closed or covered with weeds. In recent years, there have been frequent occurrences of collapsing underground bunkers that endanger private property [2]. Public works were frequently halted due to the discovery of underground bunkers during excavation [3] and

**Citation:** Hung, Y.-C.; Zhao, Y.-X.; Hung, W.-C. Development of an Underground Tunnels Detection Algorithm for Electrical Resistivity Tomography Based on Deep Learning. *Appl. Sci.* **2022**, *12*, 639. https://doi.org/10.3390/ app12020639

Academic Editor: Yosoon Choi

Received: 23 November 2021 Accepted: 5 January 2022 Published: 10 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ammunition depots excavated during public works [4]. The original data on these underground bunkers, including their planimetric positions, massing sizes, and underground depths, are difficult to obtain or simply incorrect. These difficulties and inconsistencies result from frequent changes of documents keepers and mistakes in handover during army withdrawal. Therefore, the Kinmen County Government faces a challenge in providing a detailed report on underground bunker investigation necessary when planning public works, thus delaying the planning and design schedule. If an underground bunker is discovered during the construction, the works must be halted. Without such reports, the building schedule cannot resume, which significantly delays the construction progress and affects the efficiency of government administration. In addition, the army faces similar problems of unknown or unclear locations of underground bunkers when returning the land to private owners, thus affecting the military handover schedule. Civilians also worry about whether there are any underground bunkers nearby, endangering their homes and safety.

Over the past decades, despite societal progress, urban disasters occurred frequently. The collapse of underground bunkers across Kinmen is an example of such disasters. Detailed and correct geological data of the site can enhance the effectiveness of urban disaster prevention and reduce the potential damage these disasters may cause. The geotechnical engineering investigation targets the geological conditions of the site. The conventional geotechnical investigation uses mainly the geological drilling method, but the drilling cost is high and provides only a single point of information. Sometimes it needs to be combined with the geophysical exploration method to gain the complete geological picture.

In recent years, nondestructive geophysical technology has been gradually introduced in different underground environmental surveys. In combination with document information, the information of underground "surface and space" can be obtained [5]. In the past, geophysical exploration was used in the geotechnical investigation primarily for seismic detection (refraction, reflection, and surface wave techniques). The elastic wave velocity of strata is closely bound with the material engineering characteristics [6–8], but the seismic detection is likely to be influenced by ambient vibration noise. Moreover, seismic refraction detection cannot positively detect the low-velocity layer under the high-velocity one, and the shallow layer (less than 50 m) reflection seismic detection is challenging to implement. The surface wave seismic detection is mainly 1D and 2D probing technique, and the 3D detection method is still in the preliminary research stage. Ground-penetrating radar is an electromagnetic method similar to seismic reflection detection. Here detecting the distribution of bed boundary and the localisation of depth and the depth of investigation are limited. If a Bistatic antenna is used, it can increase the detection depth capability [9–12]. The ground-penetrating radar has been used in underground structure and archaeological investigations to a great extent, for instance, in underground tunnel mining [13,14], underground pit mining [15,16], and historical building mining [17,18]. However, the ground-penetrating radar only judges the location of possible tunnel structures derived from the discontinuity of radar wave velocity. The shape, size, depth, and distribution range of underground blockhouses are unknown.

Electrical Resistivity Tomography (ERT) is a geophysical method of mapping underground structures using electrodes placed in boreholes or electrical resistivity measurements from the surface. The present electrical resistivity tomography technique can explore 2D and even 3D resistivity distribution, and the depth of investigation can be adjusted easily by the length of the measuring line. The resistivity is highly correlated with the geomaterial and groundwater characteristics, so it has gradually become one of the primary geophysical methods of geotechnical investigation. In the past decade, the ERT has been extensively used in a variety of geotechnical and environmental engineering investigation and monitoring, for instance, in geologic surveys [19,20], fault line surveys [21–25], slipping plane survey and monitoring [26–28], groundwater investigation and pollution monitoring [29,30], reservoir/dam leakage investigation and monitoring [31,32], refuse landfill leakage investigation and monitoring [33–35], and underground tunnel mining [36–38].

Kinmen Island is located between Taiwan and China. The test site is located in an abandoned barracks area in the central region of Kinmen Island, shown in Figure 1. A largescale military defensive infrastructure was built in this area, including multiple underground tunnels. The Shuangru Mountain Tunnel was built in 1949. It is 1142 m long, extending in all directions, and had a significant strategic position and value. After the troops' withdrawal, the Kinmen County Government declared this tunnel an important cultural heritage in 2019, as it had significant historical value. However, the site has been in disrepair; multiple tunnel intersections have collapsed, and it is difficult to determine the tunnel path and true underground location. This study attempts to evaluate the feasibility of using ERT to detect the basic data on the location, size, and depth of the underground tunnel.

**Figure 1.** Test site location and measurement wiring diagram.

The measurement results of ERT in the geotechnical investigation is influenced by soil water content, geological structure, groundwater level, and ambient noise, among other things. These factors affect the sensitivity of ERT measurement and spatial analysis, so this technology still faces problems in spatial resolution capability. Therefore, this study intends to construct a numerical model with the collected data to discuss the reliability of ERT results for detecting underground tunnels.

As all geologic configurations are essentially 3D, theoretically, the 3D resistivity survey is supposed to obtain the most accurate result [39]. At present, the 3D ERT survey has been studied by most scholars actively, but it is not yet as extensively used as 2D ERT. The main reason is that the Conduct 3D ERT test needs a larger area site, 3D ERT has a relatively higher cost, and field tests take a long time. Therefore, the use of the 3D ERT survey is not yet widespread.

In recent years, artificial intelligence has been extensively used in various domains [40]. This study uses the deep learning algorithm of artificial intelligence for reprocessing several 2D-ERT images from a field test. The Multilayer Perceptron (MLP) of the deep learning method is used to approximate resistivity value. More profiles are developed from the original 2D-ERT images, and the 3D resistivity value information is approximated by different profiles to establish the 3D resistivity model. This study aims to enhance the image resolution capability and interpretation capability of deep learning and discusses the feasibility of this type of artificial intelligence algorithm for enhancing 2D-ERT image resolution.

#### **2. Research Method**

#### *2.1. Electrical Resistivity Tomography (ERT)*

In terms of the measuring principle of ERT, the direct current or low-frequency alternating direct current is conducted to the ground through a pair of current electrodes, C1 and C2, to establish a man-made electric field. The electric field is measured using another pair of potential electrodes, P1 and P2, so as to measure the potential difference between P1 and P2, different configurations of electrodes (e.g., Wenner, pole-dipole, Schlumberger, dipole-dipole, and pole-pole), and the movement of electrodes that correspond to different space geometry factors. Thus, apparent resistivity can be obtained. The resistivity measured infield is not the true resistivity of the subsurface structures. Therefore, the apparent resistivity needs to be calculated via an inversion analysis to obtain the approximately real resistivity profile [39].

The direct current resistivity method includes 1D, 2D, and 3D detection methods. The advantage of a 1D survey is a rapid measurement, but the defect is that it cannot consider the transverse resistivity variation, thus affecting the reliability of measurement results and interpretation. In the last few decades, to enhance the accuracy of the electrical prospecting method in result interpretation, the researchers developed the 2D survey from the pseudo depth composed of the VES and profiling results of the 1D survey. As the 2D-ERT has a low cost and short test time, the apparent underground resistivity is obtained by wiring laid, and the resistivity of subsurface soil is obtained by appropriate inverse calculation. In recent years, more research has been devoted to the research on 3D measurement, leading to higher accuracy of 3D detection. However, the effectiveness of 3D detection is still under investigation and is limited to the measurement of space and time [41–45]. Thus, 2D detection is still the most economical and feasible measurement method at present.

To select the test site, several underground tunnel vent holes were found on the surface of a site according to the tunnel literature, indicating a possible presence of an underground tunnel. In order to evaluate the location, size, and depth data of this tunnel, this study attempted ERT detection.

Therefore, three survey lines, L1, L2, and L3, were laid in the Shuangru Mountain Barrack Field, as shown in Figure 1. The overall length of each survey line was 46 m, the electrode spacing was 2 m. The three survey lines were parallel with each other. The distance between L1 and L2 was 5 m, the distance between L2 and L3 was 8 m. The electrode configuration used Wenner Array for data collection infield measurement. This test used the SYSCAL PRO Switch 48 ground resistance instrument of France IRIS. In this experiment, the analysis software Res2dinv (version 3.54 z) [39] developed by Geotomo

was used for inverse analysis. The inverse analysis method used the optimal least square method (L2 norm). To ensure the measurement quality, each survey line was measured repeatedly to ensure the deviation was below 3%.

#### *2.2. Reliability Analysis*

To analyse the reliability of the 2D-ERT result, this study attempts to build a numerical model to validate the correctness of in situ measured data. The numerical model was built using RES2DMOD developed by Geotomo [46]. The in situ geologic configuration was simulated in the built mesh, and an appropriate resistivity value was given. The numerical solution was calculated using the finite element method to obtain apparent resistivity, and the apparent resistivity profile inverse calculation was carried out using the Optimal least square method (Res2dinv [47]). The numerical simulation resistivity profile can be obtained and compared with the in situ actual resistivity profile result.

The numerical model resistivity can be obtained from the collected geologic data, such as geologic maps, drilling data, the resistivity distribution range of geomaterial, and the preliminary detection result of ERT. Different depths, positions, or specific regions of the model are given appropriate resistivity values according to experience and professional judgment.

According to the analysis of the collected geologic data of the site, the geology around the Shuangru Mountain Tunnel mainly comprises two kinds of strata, which are silty sand within 15 m below the earth's surface and silty clay at 15–30 m below the earth's surface. This study referred to the geologic drilling data of the site and took the silty sand as background to build the numerical model of a tunnel and simulate the tectonics stratum of the tunnel.

#### *2.3. Deep Learning*

The MLP is one of the basic types of deep learning architecture. Its architecture follows the neural network system principle, learns, and predicts data. The MLP learns in the perceptron; the weight is changed after each data processing. The weight is adjusted using an algorithm, and the deviation in the training process is reduced to minimize the error in the amount of output and prediction results. The main advantage is the ability to solve complex problems rapidly. The MLP is a feedforward neural network composed of the multilayer structure of linear and nonlinear activation functions. Each layer is composed of basic elements of neurons, and each neuron is fully connected to all others in the previous and subsequent layers. Neurons in the network have a bias value *b* and an activation function *f*. The connection between neurons in different layers is defined by connection weight, *wi* ∈ *R*, *i* ∈ {1, 2, . . . , *n*}. These parameters are updated during the training of MLP. The general activation functions are ReLU, TANH, sigmoid, and softmax. The output value *y* of a neuron is defined as:

$$y = f\left(\sum\_{i=1}^{n} w\_i \mathbf{x}\_i + b\right) \tag{1}$$

where *xi* is the output value of neuron *i* of the previous layer, *i* ∈ {1, 2, . . . , *n*}. The MLP network generally includes one input layer, one output layer, and one or more dense layers, detailed below.


In the training process of MLP, the weight *wi* and offset parameter *b* are updated in every iteration, and the update target is a minimum loss function. After network training, it can be used for classification or function approximation. The architecture is shown in Figure 2 [48–52].

**Figure 2.** MLP architecture diagram.

In this study, we used Adaptive Moment Estimation (Adam) optimization to train MLP and minimize the Mean Squared Error (MSE). Adam optimization was applied to calculate the adaptive training rate of the parameters [53]. This method, in addition to storing the descending mean of the square of the past of the gradient or *vt* and the mean of the descending of the past of the gradient or *mt*, is kept as momentum. This is why the momentum can be seen as a ball sliding on a sloping surface with no friction and, therefore, can be placed at the minimum error level [54]. Both parameters of the mean descending average of the square of the gradient and the average of the descending of the past of the gradient can be calculated as follows:

$$m\_t = \beta\_1 m\_{t-1} + (1 - \beta\_1)g\_t$$

$$v\_t = \beta\_2 v\_{t-1} + (1 - \beta\_2)g\_t^2$$

where *gt* is the gradients at subsequent time steps, *mt* and *vt* are estimates of the first moment (the mean) and the second moment (the uncentered variance) of the gradients, respectively (hence the name of the method). As *mt* and *vt* are initialized as vectors of zeros, the authors of Adam observe that they are biased towards zero, especially during the initial time steps, especially when the decay rates are low (i.e., *β*<sup>1</sup> and *β*<sup>2</sup> are close to 1).

The prediction of unknown data has always been a complex problem, and it is also an important application in neural networks. The high-level network architecture of deep learning can more accurately predict and describe the distribution and changes of unknown data. The ERT profile is an inverted trapezoid; two sides are blank and short of resistivity data. This study used machine learning, data prediction, and computational ability to reduce the error between the actual value and prediction value of MLP for the blank in the field test ERT profile by MLP. The result was obtained after reiteration to convergence, forming a complete rectangular section. More profiles of different positions were developed from MLP, and the in situ 3D resistivity model was built last. In this study, we used only three-layer (discrete) 2D-ERT data to predict that multilayer (continuous) 3D ERT data is itself a complex and challenging problem. We hope that the 3D resistivity image predicted by MLP will enhance image resolution and interpretation capability.

#### **3. Results and Discussion**

*3.1. ERT Field Acquisition and Data Analysis*

#### 3.1.1. L1 Survey Line

The result of the L1 survey line, according toWenner, is shown in Figure 3 (RMS Error = 2.3%). The electrical property of the stratum is approximately divided into two layers, and the resistivity profile descends as the depth increases. There are two higher resistivity regions located at the survey line 10–18 m and 34–42 m, within 1.5 m depth from the surface. Because Kinmen Island had been free of rain for a long time in the measurement period, the topsoil was very dry. The other region shows low resistivity values at 3–18 m and 32–42 m of survey line and within 1.5–7.9 m underground. As the geology of this region is silty sand stratum, the geology is relatively loose and has higher water content. Notably, at 18–25 m of survey line and 2.5–6.5 m of depth, the resistivity value increases suddenly, the resistivity value of this region increases inwards, indicating this may be the tunnel location, according to preliminary analysis. According to Figure 3, the tunnel may be located at 2.5 m underground, with a height of about 2.5 m and a width of about 4 m, indicated by red dotted lines in Figure 3.

#### 3.1.2. L2 Survey Line

Figure 4 shows the result of the L2 survey line (RMS Error = 2.3%). It shows the same trend as Figure 3 because the L2 survey line is close to L1. The stratum is divided into two layers, in which the middle region has relatively high resistivity, while the rest is a stratum of lower resistivity. As the geology of this region is silty sand stratum, the geology is relatively loose and has higher water content. At 18–28 m of survey line and 2.5–5.5 m of depth, the resistivity value increases suddenly, the resistivity value of this region increases inwards, indicating this may be the tunnel location according to preliminary analysis. As shown in Figure, the tunnel may be located at 2.5 m underground, with a height of about 2.5 m and width of about 4 m, indicated by red dotted lines in Figure 4.

#### 3.1.3. L3 Survey Line

The result of the L3 survey line, according to Wenner, is shown in Figure 5 (RMS Error = 1.87%). It shows the same trend as Figure 4 because of the proximity of the L3 survey line to L2. The electrical property of the stratum is approximately divided into two layers. The resistivity profile descends as the depth increases, and there is a region of relatively high resistivity in the middle. At 3–14 m of survey line and within the earth's surface to 1.5 m underground, this region showed higher resistivity. Because Kinmen Island had been free of rain for an extended period during measurement, the topsoil was extremely dry. The other region showed a low resistivity value at 22–42 m of survey line and within the distance of 7.9 m from the earth's surface. As the geology of this region is silty sand stratum, the geology is relatively loose and has higher water content. At 12–20 m of the survey line and 2–5 m of depth, the resistivity value increases suddenly, and the resistivity value of this region increases inwards. This indicates a possible tunnel location, according to preliminary analysis. As shown in Figure 5, the tunnel may be located at 2 m underground, the height is about 2.5 m, and the width is about 4 m, indicated by red dotted lines in Figure 5.

#### **Figure 5.** L3 result.

#### 3.1.4. Comprehensive Interpretation

According to the measurement results of three survey lines, L1, L2, and L3, there is a high resistivity region in the earth resistivity profile. In order to analyse whether the high resistivities are correlated with each other in space, a 2.5D simulated diagram of the three survey lines in relative positions is drawn, as shown in Figure 6. It is observed that if the high resistivity regions of various profiles are selected and connected in line, the selected region and path correspond to the tunnel location and path in the literature.

**Figure 6.** 2.5D simulated diagram of three survey lines.

#### *3.2. Reliability Analysis*

To verify the reliability of the 2D-ERT images, this study constructed a numerical model to validate the correctness of in situ measured data. To simulate the onsite strata condition, this study used the same measurement parameters as in situ measurement, referred to the geologic drilling data of site, and took the silty sand as background to build the numerical model of a tunnel to simulate the stratum tectonics of the tunnel in the stratum. Referring to the site concrete tunnel size, a 3 m high and 5 m wide square tunnel was built at 2 m underground, and the reasonable resistivity of each unit was assumed. The resistivity value of the concrete tunnel was assumed to be 1000 ohm-m, dry air inside a tunnel, the resistivity was assumed to be 2000 ohm-m, the resistivity of the soil layer covering the tunnel was assumed to be 400 ohm-m. The numerical model is shown in Figure 7.

**Figure 7.** Numerical simulation model.

Figure 8 shows the numerical simulation result (RMS Error = 0.46%). As seen, the background silty sand shows low resistivity distribution, and relatively low resistivity values are shown at 2–18 m and 26–42 m of the survey line. Relatively high resistivity distribution is shown at 20–25 m of a survey line and 2–5 m of depth, indicated by red dotted lines in Figure 8. The numerical simulation result in Figure 8 is compared with the in situ test results in Figures 3–5. The four figures are very similar to each other, proving that the stratum tectonics in Figures 3–5 match the stratum tectonics assumed by numerical simulation, proving high reliability. Therefore, the high resistivity region is identified as the tunnel location and path.

**Figure 8.** Numerical simulation result.

#### **4. Deep Learning Prediction**

The 2D ERT test result mainly shows a 2D profile; the tunnel location and path cannot be interpreted directly. To establish a 3D electrical resistivity model, this study used MLP for calculation. The computational ability of learning, data prediction, and reducing the error between the actual value and prediction value of MLP was used to build the in situ 3D electrical resistivity model.

Figure 9a–g shows a series of results of deep learning prediction. Before artificial intelligence calculation, the measured data in this research must be converted into an MLP data format. The L1, L2, and L3 in situ measurement results are inputted to establish the graph of the relation of relative positions of L1, L2, and L3 for the initial model. Figure 9a shows the initial model of the original data of L1, L2, and L3; the result is similar to Figure 6. Table 1 shows the MLP network structure parameters. The MLP is a three-layer architecture, including two hidden and one output layer. Figure 9b shows the training result of the MLP network. The MLP is reiterated during training to reduce the error amount of output and prediction results and performs calculations until the network converges. Through the nonlinear iterative calculation of the MLP, as seen in Figure 9b, the error value is converged slowly after each iterative computation. Figure 9c shows the result of the original data after MLP data prediction. This is the output of the training data, which can be compared with the original data in Figure 9a. It can be seen that the training results can faithfully present the original data characteristics. Figure 9d shows the results of using the trained MLP to predict the unknown regions of L1, L2, and L3. Figure 9d shows the initially blank areas and without onsite measurement data, such as the lower left corner and the lower right corner, which can be predicted by the MLP network to predict the resistance value. Figure 9e shows the result of using the trained MLP to predict successive profiles. When the predicted number of profiles increases, the 3D resistivity model is built, and 14 profiles are predicted in this study. Due to the narrow scope of this study, the 14 layers of data are predicted by three layers of original data of the resistivity values predicted by the MLP network, which is sufficient to establish a 3D model. Figure 9f shows the 3D resistivity model result. The tunnel is a reinforced concrete structure; it is characterized by unlikely electrical conduction and high resistivity, especially in the 2D ERT result. In order to highlight the location of high resistivity, this study concealed the low resistivity colour of the 3D resistivity model in Figure 9f. After repeated tests and the resistivities are arranged, the 60% of the resistance value is taken as a boundary, the colour of resistivity value lower than 60% is hidden, and the colour of resistivity value higher than 60% is displayed. The result, after adjustment, is shown in Figure 9g; the tunnel location and path have been shown clearly.

**Table 1.** MLP network structure parameter list.


From the existing literature, we already know the approximate location and direction of the underground tunnel in reality. The results of the MLP study show that a few (three layers) 2D-ERT data to predict the 3D ERT data can be seen and that the location of the tunnel is consistent with the actual location.

**Figure 9.** *Cont*.

(**c**)

(**d**)

**Figure 9.** *Cont*.

**Figure 9.** *Cont*.

**Figure 9.** *Cont*.

(**g**)

**Figure 9.** Deep learning prediction results (**a**–**g**). (**a**) Original data of L1, L2, and L3; (**b**) Error rate of MLP network training process; (**c**) Result of the original data of L1, L2, and L3 after MLP data prediction; (**d**) Result of unknown regions of L1, L2, and L3 after MLP data prediction; (**e**) Result of successive profiles after MLP data prediction; (**f**) Result of 3D resistivity model; (**g**) MLP predicted tunnel location and path.

#### **5. Conclusions**

Most of the extensively hidden underground tunnels in Kinmen Island have been abandoned and unoccupied to avoid the damage to human life and property caused by tunnel breakdown and prevent accidental destruction of the precious cultural heritage. This study used 2D-ERT for in situ measurement to find the location of tunnels hidden underground. To test the reliability of 2D-ERT in tunnel survey, a numerical model was built. The actual configuration of the onsite strata and the shape of the underground tunnel were simulated to analyze and validate the accuracy of in situ measurement. In addition, to enhance the measurement result interpretation ability and resolution to determine the actual position and path of the tunnel, this study used an artificial intelligence algorithm to estimate and build a 3D electrical resistivity model.

The results indicate that three survey lines were used for 2D ERT measurement in this study. As the tunnel is a reinforced concrete structure, it is characterized by unlikely electrical conduction and high resistivity, especially in the 2D ERT result. The three survey lines have high resistivity regions in the profiles. According to the numerical simulation result, this high resistivity region is highly likely to be the tunnel location, and the tunnel location, height, and width have been displayed clearly. The three survey lines are displayed in 2.5D mode. The high resistivity regions of various profiles are selected and connected in line, and the selected region and path are similar to in situ tunnel location and path. According to the results of this study, using 2D-ERT to survey underground tunnels is highly feasible.

With the 2D-ERT test result showing a 2D profile, the tunnel location and path cannot be interpreted directly. As the 3D electrical resistivity survey is limited by space, time, and economy, the measurement is not yet frequently used. In recent years, artificial intelligence has been extensively used in various domains. To establish a 3D resistivity map, this study used the MLP of artificial intelligence deep learning algorithm for calculation and used the computational ability of learning, data prediction. The error between the actual value and prediction value of MLP was reduced to complete the building of in situ 3D electrical resistivity model. According to a series of findings, the tunnel location and path can be displayed using an artificial intelligence algorithm. The measurement result interpretation ability and resolution are enhanced effectively. This study can further discuss different types of tunnels in the future to enhance the feasibility of 2D-ERT for underground tunnel surveys.

**Author Contributions:** Conceptualization, Y.-C.H.; methodology, Y.-C.H. and Y.-X.Z.; modelling, W.-C.H.; software, Y.-X.Z. and W.-C.H.; data analysis, Y.-X.Z., W.-C.H. and Y.-C.H.; conclusions, Y.-C.H. and Y.-X.Z.; Field test, W.-C.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the Ministry of Science and Technology of the Republic of China, Taiwan, for financially supporting this research under Project number MOST 109-2221-E-507-004.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

