Application of Airborne Magnetic Survey in Deep Iron Ore Prospecting—A Case Study of Jinling Area in Shandong Province, China

: With the increasing demand for mineral resources, there is an inevitable trend to carry out deep prospecting in existing old mines to ﬁnd a second or even third mining space. Deep prospecting is also an affordable and practical way to prolong the lives of mines and provide a sustainable supply of mineral resources. The magnetic survey is arguably the most effective method for iron ore prospecting. In this paper, a high resolution airborne magnetic (HRAM) survey for deep iron prospecting in the Jinling iron ore cluster (JIOC) was carried out in 2018, which renewed the ﬁeld magnetic data of the JIOC obtained in the 1980s. From previous studies, almost all iron deposits in the JIOC are spatially distributed in the contact zone between the intrusive rocks and the surrounding rocks. The key prospecting areas were inferred by delineating intrusive rock boundaries via boundary enhancement and edge detection methods, and one of the areas was veriﬁed by drilling.


Introduction
Iron is the earliest, most widely used and most abundant metallic element found in the world, and, in some ways, probably the most important mineral resource in humankind's history. The magnetic method is, almost certainly, the oldest geophysical exploration technique for iron ore prospecting [1]. "Swedes successfully located a new iron ore deposit using a simple declination compass in the 1640s" [2]. "Thalen-Tiberg magnetometer, invented by Thalen and modified by Tiberg in the 1880s, was a significant step forward of the magnetic survey" [2]. As an important branch of magnetic exploration methods, the aeromagnetic survey, taking the aeromagnetic induction magnetometer developed by A.A. Logachev in 1936 as an important milestone, has made outstanding contributions to the world's iron ore exploration [1][2][3]. Up to now, aeromagnetic exploration is still an effective method for iron ore exploration which is widely used throughout the world [4][5][6][7][8][9][10][11][12][13][14]. In China, it is easy to accept that the magnetic method is also the main means of iron ore exploration and that aeromagnetic exploration plays the most important role [15][16][17][18][19][20][21]. According to statistics, about 80% of China's iron deposits, especially concealed iron deposits, are related to aeromagnetic exploration [15,22,23].
With the growing demand for steel as reform and opening-up progresses, China is now the world's largest steel-producing and consuming country. But China's rich-iron ore

Regional Geology
The study area is located at the Luxi uplift (also called western Shandong uplift) of the northeastern North China Craton (NCC). The NCC is a main tectonic unit of the Chinese mainland with a history of 3.8 billion years. It is mainly bounded by three orogenic belts. The north and southwest boundaries are, respectively, the Inner Mongolia-Daxinganling orogen and the Qinling-Dabie orogen, and the southeast boundary is the Su-Lu orogen [37][38][39]. (Figure 1a) The intrusive rocks of JIOC are mainly early Cretaceous fine-grained hornblende diorite (K). There are four main intrusive rock bodies, monzonite, biotite diorite, hornblende diorite and quartz diorite, that form a heterogeneous body of medium-basicneutral-medium-basic rocks, which is known as the "Jinling Complex"(JC). The JC is distributed in a NE direction, with an intrusion time of 128 Ma and an exposed area of about 72 km 2 , and is the source rock of mineralization of iron ores in the area [35,42] ( Figure 2).  [37][38][39]). (b) Sketch map of Shandong Province (Modified after GEISZ and [40,41]).
The fold structures of the JIOC are mainly the Jinling anticline and Hutian syncline. The Jinling anticline is N45E in strike, 20 km in length and 10 km in width. The core of the Jinling anticline is JC, and the two wings are Ordovician, Carboniferous and Permian strata. The occurrence of the Jinling anticline closer to the JC is steeper, with a dip angle of 30° to 50° in dip angle, compared to 20° to 30° to the outside. The Hutian syncline is NE-NEE in strike and located to the southeast of JC. Its northwest wing is steep with a  [37][38][39]). (b) Sketch map of Shandong Province (Modified after GEISZ and [40,41]).
The Luxi uplift, located in the western part of Shandong province, is an important extensional structure in the southeastern NCC. It is spatially bounded to the north and west by the Jiyang depression and Linqing depression of the Bohai Basin, and to the east by the Tanlu Fault Zone as a contact zone connecting the Jiao Bei uplift, the Jiao Lai depression and the Jiao Nan-Weihai uplift from north to south, respectively. The frequent activities of folds, fractures and magma contributed to making the Luxi region an important concentration of skarn-type rich iron deposits in China [36,40,41] (Figure 1b).
The study area, JIOC, as mentioned above, is several kilometres away from the northeast of Zibo city. The strata of the JIOC from old to new are Ordovician (O), Carboniferous (C), Permian (P) and Quaternary (Q). The Ordovician strata are mainly carbonate of the Majiagou Group, with an overall thickness of 700-800 m, of which the most closely related to mineralization is relatively pure limestone. The Carboniferous strata mainly consist of iron-alumina mudstone and fine sandstone, with a thickness of about 400 m. The Permian strata mainly consist of sand shale and thin layer limestone, with a total thickness of about 500 m. The flood and slope deposits of the Quaternary, mainly yellow silt clay, silt clay and silt clay, with a thickness of 60-200 m, are widely distributed ( Figure 2). displacement of 250 m-500 m, is a normal fault located on the northeast part of JC ( Figure  2).
The statistical results of magnetic properties of the JIOC area provided by GEISZ indicate that the known iron-rich ores are mainly magnetite with very strong susceptibility and remanent magnetization, which are on average 1-2 orders of magnitude higher than diorites. The diorites have medium magnetism, and other strata are nonmagnetic. The results provide the basis for this HRAM (Table 1).
The intrusive rocks of JIOC are mainly early Cretaceous fine-grained hornblende diorite (K). There are four main intrusive rock bodies, monzonite, biotite diorite, hornblende diorite and quartz diorite, that form a heterogeneous body of medium-basic-neutralmedium-basic rocks, which is known as the "Jinling Complex"(JC). The JC is distributed in a NE direction, with an intrusion time of 128 Ma and an exposed area of about 72 km 2 , and is the source rock of mineralization of iron ores in the area [35,42] (Figure 2).
The fold structures of the JIOC are mainly the Jinling anticline and Hutian syncline. The Jinling anticline is N45E in strike, 20 km in length and 10 km in width. The core of the Jinling anticline is JC, and the two wings are Ordovician, Carboniferous and Permian strata. The occurrence of the Jinling anticline closer to the JC is steeper, with a dip angle of 30 • to 50 • in dip angle, compared to 20 • to 30 • to the outside. The Hutian syncline is NE-NEE in strike and located to the southeast of JC. Its northwest wing is steep with a dip angle of 20 • to 30 • and the southeast wing is gentle with a dip angle of 5 • to 10 • (Figure 2 The statistical results of magnetic properties of the JIOC area provided by GEISZ indicate that the known iron-rich ores are mainly magnetite with very strong susceptibility and remanent magnetization, which are on average 1-2 orders of magnitude higher than diorites. The diorites have medium magnetism, and other strata are non-magnetic. The results provide the basis for this HRAM (Table 1).

Airborne Magnetic Survey System
An Airbus AS350-B3e(H125) helicopter, carrying the aeromagnetic processing & storage system, airborne magnetometer sensor, radar altimeter, pilot navigation unit (PGU) and fluxgate magnetometer, was used as the aeromagnetic platform in the survey. The aeromagnetic input & processing system, radar altimeter and PGU were installed in the cabin. The magnetometer sensor and fluxgate magnetometer were installed in the stinger which was mounted under the cabin. Meanwhile, the diurnal magnetic monitoring (DMM, also called BaseMag) station was settled in a wide-open place nearby ( Figure 3).

Survey Lines
The terrain of the survey area is low in the north and high in the south. The highest place is Tushan Hill in the south-central part of the survey area, with an altitude of 169 m ( Figure 4a). The other elevation differences are generally less than 50 m. According to the distribution of geological structures, the direction of survey lines (SL) was designed as NW-SE (135°-315°), and the distance between SLs was 100m. The direction of tie lines (TL) was perpendicular to the main survey lines with a spacing of 1.5 km (Figure 4b). In Sept. 2020, the 194 SLs and 6 TLs, with a total length of 1200 km, had been flown in 2 days.

Test Flight (1) Compensation and Heading Error Test
To improve the quality of the aeromagnetic data, it is necessary to carry out an aeromagnetic compensation flight before the survey flight and after maintenance. The compensation flight is essential and effective for high-quality aeromagnetic survey, as it reduces magnetic interference and noise caused by the maneuvers of the aircraft flying within the Earth's magnetic field, and the heading error test flight ("Clover-leaf" test flight) is an easy way to determine and reduce the error when the aircraft changes heading during flight [43][44][45][46][47]. In this survey, the compensation and clover-leaf test flight were The integrated multi-parameter airborne console (IMPAC), integrated with data acquisition, navigation and aeromagnetic compensator, produced by PICO Envirotec, was used as an aeromagnetic processing and storage system. A CS-3 high sensitivity magnetometer sensor manufactured by Scintrex Ltd. (Concord, ON, Canada). was chosen to record magnetic field readings. A CG-T3 digital radar altimeter was selected to provide the height of flight. A Hemisphere Crescent R120 was used as the GPS instrument. A BILLINGSLEY AEROSPACE & DEFENSE TFM100G2 Triaxial fluxgate magnetometer was fixed in the middle of the stinger to obtain the flight manoeuvre information of the helicopter. A PICO PGIS-HSMAG station was placed on farmland at a certain position to monitor the diurnal magnetic change. The Geosoft Oasis Montaj (Version 8.0, manufactured by Seequent Limited, Toronto, ON, Canada) was used in processing and mapping the magnetic data.

Survey Lines
The terrain of the survey area is low in the north and high in the south. The highest place is Tushan Hill in the south-central part of the survey area, with an altitude of 169 m (Figure 4a). The other elevation differences are generally less than 50 m. According to the distribution of geological structures, the direction of survey lines (SL) was designed as NW-SE (135 • -315 • ), and the distance between SLs was 100m. The direction of tie lines (TL) was perpendicular to the main survey lines with a spacing of 1.5 km (Figure 4b). In September 2020, the 194 SLs and 6 TLs, with a total length of 1200 km, had been flown in 2 days.

Test Flight (1) Compensation and Heading Error Test
To improve the quality of the aeromagnetic data, it is necessary to carry out an aeromagnetic compensation flight before the survey flight and after maintenance. The compensation flight is essential and effective for high-quality aeromagnetic survey, as it reduces magnetic interference and noise caused by the maneuvers of the aircraft flying within the Earth's magnetic field, and the heading error test flight ("Clover-leaf" test flight) is an easy way to determine and reduce the error when the aircraft changes heading during flight [43][44][45][46][47]. In this survey, the compensation and clover-leaf test flight were carried out at an altitude of 3000 m to reduce the artificial magnetic influence on the ground ( Figure 5).
(2) LAG Test A LAG test flight is used to determine the magnetic field variations in opposite flight directions and to reduce the offset effect caused by the installation distance between the magnetometer sensor and GPS antenna ( Figure 6). The terrain of the survey area is low in the north and high in the south. The highest place is Tushan Hill in the south-central part of the survey area, with an altitude of 169 m (Figure 4a). The other elevation differences are generally less than 50 m. According to the distribution of geological structures, the direction of survey lines (SL) was designed as NW-SE (135°-315°), and the distance between SLs was 100m. The direction of tie lines (TL) was perpendicular to the main survey lines with a spacing of 1.5 km (Figure 4b). In Sept. 2020, the 194 SLs and 6 TLs, with a total length of 1200 km, had been flown in 2 days.

Test Flight (1) Compensation and Heading Error Test
To improve the quality of the aeromagnetic data, it is necessary to carry out an aeromagnetic compensation flight before the survey flight and after maintenance. The compensation flight is essential and effective for high-quality aeromagnetic survey, as it reduces magnetic interference and noise caused by the maneuvers of the aircraft flying within the Earth's magnetic field, and the heading error test flight ("Clover-leaf" test flight) is an easy way to determine and reduce the error when the aircraft changes heading during flight [43][44][45][46][47]. In this survey, the compensation and clover-leaf test flight were  (2) LAG Test A LAG test flight is used to determine the magnetic field variations in opposite flight directions and to reduce the offset effect caused by the installation distance between the magnetometer sensor and GPS antenna ( Figure 6). Figure 6. LAG test flight profiles were plotted on the same X and Y axes. The anomaly centres are aligned after a 0.5 s LAG correction was applied.

Corrections
(1) Diurnal Magnetic Correction Diurnal magnetic Correction (DMC) or BaseMag Correction (BMC) is used to eliminate the influence of periodic daily variation and short-term disturbance of the geomagnetic field [48,49]. The BaseMag Station (BMS) is generally placed nearby the survey area where the magnetic field is stable (Figure 7). In this survey, the BMS was  (2) LAG Test A LAG test flight is used to determine the magnetic field variations in opposite flight directions and to reduce the offset effect caused by the installation distance between the magnetometer sensor and GPS antenna ( Figure 6).

Corrections (1) Diurnal Magnetic Correction
Diurnal magnetic Correction (DMC) or BaseMag Correction (BMC) is used to eliminate the influence of periodic daily variation and short-term disturbance of the geomagnetic field [48,49]. The BaseMag Station (BMS) is generally placed nearby the survey area where the magnetic field is stable (Figure 7). In this survey, the BMS was located on a farm between the survey area and the landing point. The DMM and aeromagnetic were both set sampling rates at 10Hz. The DMC was computed point by point according to GPS time.

Corrections (1) Diurnal Magnetic Correction
Diurnal magnetic Correction (DMC) or BaseMag Correction (BMC) is used to eliminate the influence of periodic daily variation and short-term disturbance of the geomagnetic field [48,49]. The BaseMag Station (BMS) is generally placed nearby the survey area where the magnetic field is stable (Figure 7). In this survey, the BMS was located on a farm between the survey area and the landing point. The DMM and aeromagnetic were both set sampling rates at 10Hz. The DMC was computed point by point according to GPS time.
(1) Diurnal Magnetic Correction Diurnal magnetic Correction (DMC) or BaseMag Correction (BMC) is used to eliminate the influence of periodic daily variation and short-term disturbance of the geomagnetic field [48,49]. The BaseMag Station (BMS) is generally placed nearby the survey area where the magnetic field is stable (Figure 7). In this survey, the BMS was located on a farm between the survey area and the landing point. The DMM and aeromagnetic were both set sampling rates at 10Hz. The DMC was computed point by point according to GPS time.  (2) Lag Correction According to the flight results of the LAG test (Figure 6), the magnetic field of different directions is corrected to make the magnetic field value of each measuring point as close to the actual space position as possible.

(3) Background Removal
In aeromagnetic measurement, the measured magnetic field value can simply be understood as the vector sum of the normal geomagnetic background field and the magnetic anomaly field, and removing the normal geomagnetic background field can obtain magnetic anomaly information [46]. The International Geomagnetic Reference Field (IGRF) was established in 1965 by the International Association of Geomagnetism and Aeronomy (IAGA) and updated every 5 years [50,51]. IGRF correction is now the most common and widely used method of normal geomagnetic background field removal.
After the above processing and correction of the measured aeromagnetic data, the commonly-known anomaly magnetic map was obtained (Figure 8a).

Reduction to the Magnetic Pole
To eliminate the influence of geomagnetic inclination and declination on magnetic anomaly data, reduction to the magnetic pole (RTP) should be applied to convert inclination to the vertical and declination to the north. RTP can effectively eliminate the repeated superposition of positive and negative anomalies in the original data, thus reducing the complexity of the anomalies and making them better correlate with the geological units. Most of the anomalies are narrowed down after RTP processing, and the morphology becomes simpler and easier for interpretation. RTP processing was according to the geomagnetic inclination of 55.35 • and geomagnetic declination of −6.69 • in this region.
The colour-shaded grid maps of magnetic anomaly ( Figure 8a) and RTP (Figure 8b) show the excellent correspondence between the aeromagnetic data and the geological units. As clearly seen, Figure 8 shows that the strong positive magnetic field generated by the JC and iron ores present an expected result of this survey, according to the magnetic parameters (Table 1).

Low Pass Filter
Due to the superposition of different wavelength fields of geological bodies of different depths and scales, some anomalies become very complicated and are difficult to identify, which creates difficulties for geological interpretation. In the Jinling area, artificial objects such as factories generate many high-frequency magnetic signals, which would, more or less, generate interference in the boundary enhancement and other work. Therefore, a 3 × 3 convolution filter was applied to the RTP to properly suppress the high-frequency magnetic signals.
to the geomagnetic inclination of 55.35° and geomagnetic declination of −6.69° in this region.
The colour-shaded grid maps of magnetic anomaly (Figure 8a) and RTP (Figure 8b) show the excellent correspondence between the aeromagnetic data and the geological units. As clearly seen, Figure 8 shows that the strong positive magnetic field generated by the JC and iron ores present an expected result of this survey, according to the magnetic parameters (Table 1).

Boundary Enhancement and Edge Detection
Boundary enhancement often refers to the conversion of measured magnetic anomaly data using numerical conversion techniques, to highlight numerical features near the source boundary in the converted data for further use in the identification of the source boundary (edge detection).
Derivative analysis is the major research direction of boundary feature enhancement, which focuses on the various characteristics of the derivatives of anomalous data and has a long history and variety of methods. Methods such as the vertical derivative (VDR) [52] and the second-order vertical derivative (SVDR) [53][54][55] of the total field or components data can tighten the anomaly, and the zero value of VDRs refers to the edges of the body. The VDR and SVDR are formulated as follows: where T is the magnetic anomaly field and z is the vertical derivative direction. Note that as the depth of the field source increases, the magnetic anomaly becomes wider and slower, leading to an outward expansion of the zero value points of VDR and SVDR, thus deviating from the true geological boundaries. Conversely, the maximum value of the total horizontal gradient (THD) is located at the anomaly boundary [56,57]. In addition, the analytic signal (AS) method is also a basic method of boundary enhancement, and its maximum values can be used as markers to identify the boundary of a field source [58,59]. The equations of THD and AS are as follows: where T is the magnetic anomaly field, ∂T ∂x , ∂T ∂y and ∂T ∂z are the derivatives in x, y and z directions, respectively. THD and AS both have the disadvantages of being insensitive to weak anomalies adjacent to strong anomalies and deeply buried geological bodies.
Since the magnetic anomalies from deep field sources are often covered by shallow ones, the tilt angle method (Tilt) [60] has been proposed to identify the boundaries of sources at different depths: which is less affected by the burial depth of the field source and can simultaneously enhance the boundaries of sources with different anomaly amplitudes, but the position of the localized boundaries is seriously affected by the magnetization direction.
Based on the above methods, many improved methods have been developed, such as the following: The total horizontal derivative of the tilt angle (THD_Tilt) [61]: where Tilt is calculated in Equation (5). The maximum value of THD_Tilt is located at the edge of the body. THD_Tilt can approximately recognize the edges of deep anomalies but is not capable of separating anomalies very close together. The Theta method (θ) [62]: where VDR and AS are calculated in Equations (1) and (4), respectively. The value of Theta is between 0 and 1, and it uses the maximum value to locate the boundary position. The method is less affected by the burial depth of the field source and can enhance both deep and shallow anomaly boundaries, but it has the disadvantages of sensitivity to noise, and estimated edges that are more dispersed and wider than they actually are. The analytic signal of Tilt angle (AS_Tilt) [63]: which uses its maximum value position to locate the boundary. However, the method is extremely sensitive to noise and does not enhance the boundary well enough. The Tilt angle of analytic signal (Tilt_AS) [64]: The vertical derivative of analytic signal (VDR_AS) [65]: where both use their maximum value to delineate the boundary. Since they use AS as the input source, the two methods also suffer from insufficient enhancement of deep field sources.
Researchers are continuously comparing the advantages and disadvantages of boundary enhancement and edge detection methods and proposing newer methods [66][67][68][69]. From these studies above, we selected several commonly used methods and the STDR filter presented by [68] to detect the boundaries of HRAM in the JC area ( Figure 9).
Minerals 2021, 11, x 11 of 16 STDR for identifying the submarine geological boundaries in the northern Scandinavia coastal zone is shown in paper [68], where the detected boundaries are highly accurate and well recognized. The total horizontal derivative of STDR (THD_STDR) was also introduced in the paper [68] to get a better presentation of the detected boundaries.   The STDR filter was introduced as follows: where SVDR and THD are calculated in Equations (2) and (3), respectively. M is a constant given by the average total magnetic field intensity of the region. The maximum and minimum values of STDR delineate the range of positive and negative anomalies, respectively. STDR normalises the intensity of deep and shallow anomalies as well as different anomalies distributed very close to each other. A really excellent application of STDR for identifying the submarine geological boundaries in the northern Scandinavia coastal zone is shown in paper [68], where the detected boundaries are highly accurate and well recognized. The total horizontal derivative of STDR (THD_STDR) was also introduced in the paper [68] to get a better presentation of the detected boundaries. From Figure 9, it iss easy to see that VDR, THD, AS and Tilt_AS have poor effects on the two weak banded anomalies in the northern part of the study area and AS completely suppressed them. All methods reflect that the north part of the JC (NJC) is composed of multiple small rock masses. But VDR, THD, AS and Tilt_AS do not reflect the small rock body boundaries clearly. This is because the northern part of the study area is deeply covered and the anomalies generated by deep field sources are weak, and the NJC is composed of magnetically variable rock and ore bodies, producing complex magnetic anomalies of varying intensity and proximity. As described above, these four methods have insufficient resolution for the deep weak anomalies and adjacent anomalies. Comparatively, Tilt, Theta, STDR and THD_STDR have better descriptions of the boundaries. Consequently, by comparison, the performance of STDR is the most balanced, which not only enhances weak anomalies but also shows most internal details of the JC. Further, THD_STDR shows the edges with a higher resolution. Therefore, STDR is used to infer the JC boundaries, and mineralization areas are proposed accordingly (Figure 10a).

Results and Discussion
According to previous studies, the JC was formed by magma intruding from the north-east to south-west along with the fracture structures, with the intrusion becoming progressively shallower in depth [34]. Figure 10a shows the inferred distribution (blue dots area) of the rocks based on the newly measured HRAM data. The contact zones between the JC and the surrounding strata control the morphology of the iron ore bodies and indicate favourable locations for mineralisation. Therefore, we recommend several deep prospective areas in the outer edge of the inferred NJC internal minor rock bodies (zones 1-3 in Figure 10a) and the inferred deep intrusive rock zone in the northern part of the study area (zones 4-5 in Figure 10a). Deep exploration is suggested for implementation in zones 1-3 as a priority, as they are located within the JIOC. Zones 4-5 are the inferred concealed rocks areas that have not been previously investigated, therefore it is recommended that prospecting should first be carried out in Zone 4 and then enlarged to Zone 5 after some progress has been made. deep prospective areas in the outer edge of the inferred NJC internal minor rock bodies (zones 1-3 in Figure 10a) and the inferred deep intrusive rock zone in the northern part of the study area (zones 4-5 in Figure 10a). Deep exploration is suggested for implementation in zones 1-3 as a priority, as they are located within the JIOC. Zones 4-5 are the inferred concealed rocks areas that have not been previously investigated, therefore it is recommended that prospecting should first be carried out in Zone 4 and then enlarged to Zone 5 after some progress has been made.  Engineers from GEISZ carried out a drill verification firstly in area 1 in 2020 and have successfully discovered two new iron ore layers below a depth of 550 m (Figure 10b). The upper layer is 6.3 m thick and the lower one is 2.9 m thick. Both ore bodies are located in the skarn belt which developed along the intrusive rocks and carbonate contact zones. Based on the intrusive rock boundaries inferred from HRAM, the location of the ZK drill hole is on the eastern edge of the intrusive rock body and the main part of the newly discovered ore body should be located on the eastern side of the drill hole.

Conclusions
Based on the statistical mineralisation rule that all of the iron ore deposits of the JIOC are distributed in the contact zones between intrusive and carbonate rocks, the HRAM data was used to recognize the intrusive rock edges using boundary enhancement techniques, and new ore bodies were found at the outer border of the inferred rock masses. This has demonstrated the effectiveness of aeromagnetic surveys in the search for deep magnetic deposits. Due to various reasons, this deep iron ore prospecting has not been carried out with fine interpretation such as inversion, nor has it been possible to verify the inferred deep rock bodies in the northern part of the study area. It is recommended that subsequent researchers in this area could make better use of the newly measured HRAM data for fine processing interpretation such as 3D inversion to make greater breakthroughs.