logo

Research on dynamic three-dimensional terrain correction methods of quantitative inversion for airborne gamma-ray spectrometer

ACCELERATOR, RAY AND APPLICATIONS

Research on dynamic three-dimensional terrain correction methods of quantitative inversion for airborne gamma-ray spectrometer

He-Xi Wu
Wei-Cheng Li
Rui Qiu
Chao Xiong
Yi-Ming Lyu
Yi-Qiang Xing
De-Hao Zhang
Zong-Shuo Tao
Yang Wang
Nuclear Science and TechniquesVol.37, No.5Article number 84Published in print May 2026Available online 11 Feb 2026
11500

Aerial surveys are dynamic and continuous processes, and there are different height distributions of the ground in the measurement area, which leads to problems such as overlapping measurement areas and inaccurate altitude correction during the survey process. Commonly used terrain correction methods are based on the concept of finite elementization of ground surface radioactive sources, using GPS coordinates, radar altitude, and ground elevation distribution information from aerial surveys, combined with the sourceless efficiency calibration method to construct a response matrix, which is then inverted for surface nuclide content. However, most of the sourceless efficiency calibration methods used are numerical calculations that consider the body detector as a point detector and do not consider the changes in intrinsic detection efficiency under different incident directions of gamma rays. Therefore, when the altitude of the measurement area varies significantly or the flight altitude of the aerial survey is relatively low, such sourceless efficiency calibration method calculations tend to have a large bias, which affects the accuracy of the terrain correction. To address the above problems, this study employs a novel sourceless efficiency calibration method based on the Boolean operation of the ray deposition process and simplifies the traditional body source measurement model to a surface source measurement model to achieve fast and accurate efficiency calibration. Then, through the discretization of the measurement process, the static measurement process is superposed as equivalent to the dynamic measurement process, and the dynamic measurement response matrix is built and optimized based on the calibration method. Finally, the PSO-MLEM algorithm was used to solve the dynamic measurement response matrix to achieve dynamic terrain correction of aerial survey data. Analysis of the Baiyun’ebo test area revealed that, after applying dynamic terrain correction, the inverted anomalies in uranium (eU), thorium (eTh), and potassium (K) concentrations were closer to ground measurements (within 5.72%–30.79%) and exhibited clearer anomaly boundaries compared to traditional height-based corrections. However, owing to the inherent statistical fluctuations and characteristics of matrix inversion, higher measurement values tend to absorb lower ones, potentially enlarging the anomalous regions. Nevertheless, the high-anomaly regions after inversion largely coincided with the ground truth validation, demonstrating that the proposed method can effectively correct airborne gamma spectrometry data.

Airborne gamma-ray spectrumDynamic Three-dimensionalTerrain correction
1

Introduction

The airborne gamma spectrum is a method for achieving the characterization or quantitative analysis of radionuclides by processing environmental gamma-ray measurement data [1]. As a key step in processing airborne gamma spectrum data, the inversion method of surface nuclide content has a crucial impact on the accuracy of airborne measurement results [2, 3].

In China’s nuclear industry standard, the height correction of the energy window counts of the energy spectrum is often performed by combining the ground standard source calibration data with the ray attenuation principle and elevation data during the flight [4]. The conventional height correction is only parameterized by the radar height, without considering the influence of the actual ground undulation and the uneven distribution of nuclide content, resulting in an eU content error exceeding 30% in rough terrain [5] and a maximum eTh content error of 21.5% across nine typical terrain types [6]. Xiong [7] established a 2D terrain influence coefficient method correction model by dividing the undulating surface into multiple oblique sources and then accumulating them, while Wan et al. [8] obtained the optimal results of the terrain influence coefficient method computation based on the segmentation of the body model after analyzing the model in depth. Bai [9] carried out terrain correction under the 3D spatial information subdivision of the aerial survey process through Google Earth. Brian et al. [5] proposed the use of a response matrix for three-dimensional quantitative inversion and the use of finite element concepts to correct the interactions between the survey area and the construction of an airborne gamma spectrum response matrix through the radar height and the height of the survey area terrain, so that the inversion results have been significantly improved. This method was also validated for the inversion of 137Cs activity in the environment by Moudud et al. [10]. However, the actual airborne gamma-ray survey is a continuous process, and the energy spectrum of a recording point is a synthesis of the contributions of a continuous section of the survey area under continuous flight during the measurement time. Therefore, Sun and Liu et al. [11-13] proposed a method of combining aerial flight paths with the sourceless efficiency calibration method to construct the dynamic terrain corrected response matrix, and the validity of this method was verified in aerial surveying experiments in the LSS area of Gansu and the EGRY area of Inner Mongolia.

With the research of many scholars, the theory of dynamic terrain correction has been gradually improved, so the main difficulty of the terrain correction method depends on the detection efficiency calculation of radioactive finite elements [14]. The sourceless efficiency calibration method can be divided into three categories: numerical calculation, Monte Carlo simulation, and a combination of the two [15]. Wu [16] applied the numerically calculated sourceless efficiency calibration method on the airborne gamma spectrum, and found that the numerical calculations were in good agreement with the measured values after comparison. Xiong et al. [17] derived a sourceless efficiency calibration method based on numerical integration through the axial symmetry of the geometry of the detector and the object to be measured and the principle of superposition of the gamma field, which proved to be accurate in the measurement experiments of NaI and HPGe detectors. Zhang et al. [18] constructed a numerical model of ray deposition by studying the deposition process of rays in the sensitive body of the detector combined with the probability of various types of reactions, which was verified on a point source experiment of CeBr3 detector. Numerical calculations have the advantage of being fast and easy when performing sourceless efficiency calibration; however, it is difficult to deal with the effect of multiple scattering within the detector on the photoelectric peaks, and in the face of complex measurement systems, numerical calculation methods are often difficult to provide the analytical equations or need to optimize the model to a certain extent. Therefore, the value of numerical calculations is less credible, and the calculation process is complicated.

The second type of sourceless efficiency calibration method, based on Monte Carlo simulation, is often applied to simulate difficult measurement conditions owing to its modeling flexibility [19]. Noack [20] explored how to reduce the error caused by deep penetration in Monte Carlo simulations and found that the bias method can effectively reduce the simulation error. Ghassoun et al. [21] explored the relationship between point-source and surface-source detection efficiency, measurement distance, and detector size using a Monte Carlo simulation, providing a theoretical basis for optimizing the simulation of large-volume body sources. Zhao et al. [22] proposed a geometrical transformation method of converting a body source into a line source for the simulation of detection efficiency of large body sources in Monte Carlo simulation, which is useful for the reduction of variance in the simulation results of large body sources. Yang et al. [23] tested the accuracy of the sourceless efficiency scale performance of LabSOCS and ANGLE under the 241Am source and proved its HPGe gamma-ray measurement validity. Frosio et al. [24] optimized the model in ISOCS/LabSOCS using the effect of each geometry in the measurement model on the uncertainty of the measurement results of the radioactive source and showed that the method can reduce the uncertainty by a maximum of eight times. Although Monte Carlo simulations can eliminate the effect of multiple scattering on the photoelectric peaks, their computational efficiency does not meet the practical needs when simulating complex measurement systems.

The third type of combined method utilizes the advantages of numerical computation in solving the incident ray penetration distances at different angles and Monte Carlo simulation in handling the complex particle transport process in the detector to achieve the complementary advantages of the first two methods. Zhang et al. [25] obtains the angular injection distribution of rays from an infinite body-like source by numerical computation, and combines it with a Monte Carlo simulation to simulate the detection efficiency under different angular distributions, energies, and heights. The results proved the high accuracy of the method in the measurement of natural nuclides such as U-series, Th-series, and K after ground and flight calibration experiments in Shijiazhuang. Based on the above method, we improved the Monte Carlo simulation part of the above process by simulating the deposition process of gamma rays in the detector using Geant4 and combining it with Boolean operations to quickly obtain the detection efficiency of the rays in different incidence modes, which improved the efficiency of this method by at least 162 times compared with that of the traditional Monte Carlo simulation. In the point and body source validation experiments of the CeBr3 detector, the error was kept within 5% [26].

In summary, this study introduces a novel sourceless efficiency calibration method based on Boolean operations in the ray deposition process to calibrate the static airborne gamma-ray survey model. Based on this, multiple static measurement models were superimposed to approximate the dynamic measurement model, upon which the dynamic terrain response matrix was constructed. This approach provides a solid theoretical foundation for the accurate inversion of surface radionuclide concentrations in aerial survey data.

2

Method

2.1
Sourceless efficiency calibration method

In this study, a sourceless efficiency calibration method based on the Boolean operation of the ray deposition process was used to calculate the detection efficiency of the airborne gamma spectrometer [26]. The interaction between rays and matter for analysis, rays from the radioactive source to be recorded by the detector energy spectrum photoelectric peaks of the process can be split into two stages: 1) before the ray enters the detector sensitive body, the attenuation of the ray produced by each layer of shielding (including the detector envelope), called the AE part; 2) when rays enter the detector sensitive body, the probability of depositing the energy completely (intrinsic detection efficiency) due to Compton scattering, the electron pair effect, and the photoelectric effect, is called the part. Therefore, the formula for the detection efficiency of a radioactive source can be expressed as Eq. (1).pic (1)where ηE,V is the photoelectric peak detection efficiency of the radiation source, Ω is the angular range of the ray emission, and V is the spatial range of the radiation source.

The parameters involved in the shielding process include the point source location distribution V, shield geometry and material information (including air, airplane floor, detector mounting box, detector housing, vibration damping, and reflective layers), ray energy E, and ray angular distribution Ω (elevation and azimuthal angles). To solve the AE part quickly and flexibly, the continuous variables in the measurement environment can be discretized (Fig. 1 shows a schematic diagram of the discretization of the detector sensitivities of the airborne gamma spectrometer GR820, which consists of 15 NaI(Tl) detectors with dimensions of 10.16 cm × 10.16 cm × 40.64 cm.)

Fig. 1
(Color online) Scatter diagram of the sensitive body surface of the GR820
pic

After the same discretization of each shielding body, the angle and distance relationship between the coordinates of each scattering point and the equation of the ray in a certain emission direction were used to obtain the ray penetration distance dl in each shielding body. Thus, Eq. (1) can be converted as follows:pic (2)where μE,l is the line absorption coefficient of the shielding body l for rays of energy E, obtained using the winXcom software.

For the part, Geant4 can be used to obtain and record the deposition process of NE particles with energy E in the infinitely large sensitive body at a specific incidence angle and position as the original coordinates Dn, and then spatially transform the original deposition coordinates according to the direction of ray incidence (ϕin, θin) and the coordinates of the point of incidence Sin, as expressed in Eq. (3).pic (3)where are the deposited coordinates after the transformation of the coordinate system, and Rz and Rx are the rotational multipliers of the coordinate system around the z-axis and x-axis, respectively.

Then, is intersected with the sensitizer in a graphical Boolean operation, as shown in Fig. 2 (using the rectangular detector of the GR820 as an example).

Fig. 2
(Color online) Schematic of Boolean operation of with a single sensitive body
pic

Points located inside the sensitive body are set to 1, and those outside the sensitive body are set to 0, and are stored in subscript order in the column vector Pn. Thus, the intrinsic detection efficiency , for which all the deposition processes occur inside the sensitive body, can be computed using Eq. (4).pic (4)where represents the number of rays that are fully deposited within the detector, Pn is the number of particles with full energy deposition during the n-th reaction, ones is a unit column vector with the same length as Pn, and NE is the total number of incident particles. Finally, the detection efficiency can be solved by combining Eq. (4) with Eq. (2). During efficiency calibration for narrow-beam gamma rays in the energy of 59.54–2620 keV, the computational efficiency of the method is approximately 900–7200 times higher than that of conventional Monte Carlo simulations, and it further increases with longer ray path length (CPU @ R5800H).

2.2
Simplification of the calibration model

The ground gamma-ray injection rate model usually simplifies the ground gamma-ray source as an infinite circular table body source, and the detector as a point detector for injection rate calculation [27]. However, as a 3D model, the body source model increases the calculation time for the detection efficiency. Therefore, based on the concept of the radioactivity spectrum balance, the body source model can be simplified to a surface source model to improve the calculation efficiency. In the quantitative inversion of the airborne gamma spectrum, the main focus is on the counts in the range of U-series, Th-series, and 40K characteristic peaks; therefore, the simplification process only analyzes the balance process of the characteristic gamma rays of each nuclide. Gamma rays are emitted from underground rocks with a uniform distribution of radionuclides and density ρ. The characteristic gamma-ray injection at point O can be considered to be formed by the superposition of gamma rays emitted by the spherical shell rock with r=0~R. Therefore, the injection at point O can be expressed as:pic (5)where the rock is considered a homogeneous medium with density ρ (g/cm3), radionuclide content q (g/g), absorption coefficient of gamma rays μ (cm-1), gamma ray emission probability P, and decay constant λ.

From Eq. (5), when R tends to infinity, the injection rate at point O is equal to the injection rate of the infinite body source and is isotropic in its spatial distribution. Therefore, the surface injection under the infinite radioactive body source is numerically equal to that at point O in this case, and theoretically, each point conforms to the stereoscopic angular distribution of the 2π point source. And in Eq. (5), all parameters are constants except R. Taking the surface injection rate at an infinite body source as a benchmark, the natural gamma-ray surface injection rate saturation versus R for 1.46, 1.76, and 2.62 MeV is shown in Fig. 3 (with a medium density of 2.2 g/cm3 and a substrate of SiO2, for example).

Fig. 3
(Color online) Variation pattern of gamma-ray surface injection rate with body source radius
pic

The natural gamma ray surface injection rate reaches 95% of the gamma ray surface injection rate under an infinite body source before R=50 cm. In reality, the range of lithological changes and the range of aerial surveys are generally much larger than 50 cm; therefore, this study simplified the aerial gamma body source model to the surface source model and carried out dynamic terrain correction based on the surface sourceless efficiency calibration, while the relationship between the content of natural radionuclides on the surface and the surface injection rate can be expressed by Eq. (6).pic (6)

2.3
Dynamic terrain correction method
2.3.1
Dynamic terrain aerial survey forward model

In the aerial survey process, the distribution of surface radioactivity can be divided into finite elements according to influencing factors such as measurement point location, time, and measurement distance (Fig. 4, with nj parts of latitude and nk parts of longitude). When the spectrometer flights from one recording point ti-1 to the next recording point ti, its measurement range moves dynamically with the spectrometer; therefore, all radioactive finite elements inside the total measurement range of a single flight will contribute to the spectrometer counts, and the detection efficiency of each finite element changes with the spectrometer coordinates.

Fig. 4
(Color online) Aerial survey schematic
pic

Therefore, in combination with the simplified model in Sect. 2.2, the counts of a single spectrum can be considered as a superposition of the contributions of all radiation field injection rates within the measurement range of the spectrometer during this flight segment. Moreover, because there is a difference in the radioactivity content of each geologic body, the spectrum count rate Ci of the i-th measurement can be expressed by Eq. (7):pic (7)where t is the measurement time, [x, y] is the measurement range of the ith measurement, A is the ray injection of a single radioactive finite element, and η is the detection efficiency in the corresponding measurement coordinates at time t.

To ensure the accuracy of the inversion results in aerial surveys, the range of radioactive finite elements is often larger than the accuracy of the map elevation, resulting in the terrain relief of each radioactive finite element. Moreover, the intersection region of the measurement range and the corresponding radioactive finite element changes dynamically over time. Therefore, after discretizing the single measurement time t and converting the inversion object into a radioactive finite element A(j, k) containing a height distribution, the contribution of the radioactive finite element to the counts of the spectrometer at the ith measurement can be expressed by Eq. (8).pic (8)where is the intersection area of the (j, k) radioactive finite element with the measurement range during the i-th measurement relative to time, and η(ti, j, k) is the average detection efficiency of the intersection area relative to time. Furthermore, the spectrometer count at the i-th measurement can be converted from Eq. (7) to Eq. (9).pic (9)Therefore, the entire measurement process from the first measurement to the n-th measurement can be expressed in the matrix form of Eq. (10).pic (10)where B is the response of the radioactive finite element to the count rate of the spectrometer.

In radioactive measurements, radionuclide decay is considered a random process. Consequently, even under constant source and measurement conditions, the count rate of the spectrometer is not fixed but exhibits statistical fluctuations. If a nucleus decays with a probability rate τ, then the probabilities of decay and non-decay within a time interval t are (1 - e-τ t) and e-τ t, respectively. When the spectrometer has a detection efficiency η, the probability p that a decay is both emitted and detected, and the probability q that it does not produce a count can be defined as shown in Eq. (11).pic (11)Thus, the radionuclide counts obey a binomial distribution. Given measurements on N0 radioactive nuclei, the spectrometer count ξ = n at time t can be expressed as Eq. (12).pic (12)During aerial surveys, N0 is much greater than p. By applying a limiting approximation to Equation (12), it can be shown that P(ξ = n) follows a Poisson distribution with a mean λ = N0 p, as expressed in Eq. (13).pic (13)In summary, after incorporating the background count rate δ and Poisson statistical fluctuations, the dynamic terrain aerial survey forward model can be expressed by equations such as Eq. (14).pic (14)Therefore, this study employs the PSO-MLEM algorithm [28] based on Poisson statistical fluctuations to solve the matrix.

2.3.2
Construction of the Dynamic Terrain Response Matrix

As derived in Sect. 2.3.1, the construction of the dynamic terrain response matrix B is achieved by superimposing the contributions of each radioactive finite element under static measurement conditions to obtain their total contributions under dynamic measurement conditions. As shown in Fig. 5, during static measurements, the intersection area between each radioactive finite element and the measurement range varies with the detector position. Moreover, terrain variations within each finite element also significantly affect the detection efficiency. Therefore, the main challenges in constructing the response matrix are as follows: (1) determining the static measurement range and obtaining the spatial information of all radioactive finite elements within it; (2) rapidly computing the detection efficiency of surface sources after acquiring the spatial distribution of the finite elements; and (3) determining the minimum number of convergent segments required to approximate dynamic measurements through the superposition of static measurements.

Fig. 5
(Color online) Static measurement schematic
pic

Among naturally occurring radionuclides, the 2.62 MeV gamma ray from the Th series nuclide 208Tl undergoes less attenuation in matter compared to the 1.46 MeV gamma ray from 40K and the 1.76 MeV gamma ray from the uranium series nuclide 214Bi. Consequently, its detection efficiency is significantly affected by the measurement range and incident angle within the sensitive body of the detector. Because the response matrix construction process is largely consistent across different gamma-ray energies, this section focuses on constructing the response matrix using the 2.62 MeV gamma-ray as an example. The detector considered is the GR820 airborne gamma spectrometer, with a flight altitude of 100 m.

1. Spatial information of radioactive finite elements within the static measurement range. First, the static measurement range of the detector was determined. Based on the previously described sourceless efficiency calibration method, a model of the GR820 detector was constructed. Assuming a measurement altitude of 100 m, the surface detection intensity for 2.62 MeV gamma rays was calculated over a series of radial distances ranging from 100 m to 1000 m in steps of 100 m. The resulting trend of the detection efficiency as a function of the measurement radius is shown in Fig. 6, with the 98% saturation point indicated by a red asterisk.

Fig. 6
Detection intensity vs. measurement radius at 100 m altitude
pic

After determining the measurement range at a given flight altitude, it is necessary to identify the area and spatial distribution of the radioactive finite elements within this range that contribute to the spectrometer count. As shown in Fig. 8, each radioactive finite element exhibits an altitude distribution. To account for altitude variation within each finite element, a vertical correction is applied by further subdividing the element into altitude subelements—regions with different elevations but the same radioactive activity. The detection efficiency of each altitude subelement was calculated individually based on its intersection with the measurement range. The contributions from all altitude sub-elements belonging to the same radioactive finite element are then aggregated to determine the total contribution of that element to the spectrometer, as shown in Fig. 7.

Fig. 7
(Color online) Schematic diagram of height correction for radioactive finite elements
pic

To efficiently record the intersection information between each altitude subelement and the measurement range (including the intersection area and its geographic coordinates), the measurement range is simplified as a square with side length 2R and center coordinates (Ox, Oy). Using Eq. (15), the vertex coordinates of each altitude subelement (xm,n,l, ym,n,l) are processed to compute their intersection with the square measurement range. The resulting intersection regions are all rectangles (where l=1 or 2 indicates the minimum and maximum values of x or y in the rectangle). Consequently, for each altitude subelement within the measurement range, the intersection area Sm,n and the coordinates of the four rectangular vertices were recorded. These parameters are essential for the subsequent calculation of the detection efficiency.pic (15)In this equation, sort refers to arranging the elements in ascending order. The two middle values from the sorted set are assigned as the minimum and maximum x and y coordinates of the intersection rectangle vertices.

In addition to obtaining the intersection information of the altitude subelements within the measurement range, an element with a smaller angle can still be hidden, although the terrain must be very rugged for this, as illustrated in Fig. 8. Because gamma-ray attenuation in rock and soil is significantly greater than that in air, the contributions from obstructed altitude sub-elements during each measurement were set to zero. To determine whether an altitude sub-element is obstructed, the angle θ between the sub-element and the spectrometer is first calculated. Then, all other altitude sub-elements along the line of sight between the sub-element and the spectrometer were evaluated. If any of these intermediate sub-elements have an angle to the spectrometer greater than θ, the original sub-element is considered obstructed, and its contribution is set to zero. If all the intermediate angles are smaller than θ, the subelement is considered unobstructed.

Fig. 8
(Color online) Schematic diagram of measurement obstruction
pic

2. Fast computation of surface source detection efficiency during aerial surveys. After obtaining the spatial information of each altitude sub-element that contributes to the spectrometer count, the surface source detection efficiency for each valid altitude sub-element must be calculated. However, thousands of measurement records are typically generated during practical aerial surveys. Each record must be decomposed into multiple static measurement superpositions, and each static measurement involves dozens of altitude subelements. This results in an extremely large computational workload when constructing the dynamic terrain response matrix.

To address this issue, the characteristics of the detection efficiency distribution within the measurement range during aerial surveys were analyzed to optimize the efficiency-computation process. Owing to spatial symmetry, the distribution within the first quadrant of a 500 × 500 m area was sufficient to represent the detection efficiency pattern within a 500 m radius at a flight altitude of 100 m. A total of 250,000 points were calculated with a grid resolution of 1 m. The resulting detection efficiency distribution is shown in Fig. 9.

Fig. 9
(Color online) Detection efficiency distribution of the GR820 spectrometer at 100 m altitude (2.62 MeV)
pic

From Fig. 9, it can be observed that although the overall distribution of detection efficiency resembles a two-dimensional Gaussian distribution, a closer inspection of local regions (e.g., the area within [100–160, 60–140] in Fig. 9) reveals that the distribution is approximately planar. Within this region, the detection efficiency at the geometric center can be considered an effective approximation of the average detection efficiency for the entire planar area. To identify sub-regions that can be approximated as planar segments, a profile line was drawn outward from the center, and piecewise linear fitting was performed, as shown in Fig. 10(a). Segments with a coefficient of determination R2 > 0.99 were retained. Using this method, the 500 m-radius measurement range was divided into 12 linearly fitted segments, each corresponding to a region where the detection efficiency exhibited approximately linear variation. The geometric center of each segment was then used to approximate the average detection efficiency within that region. Consequently, the entire measurement range was partitioned into 12 quasi-linear regions, as illustrated in Fig. 10(b).

Fig. 10
(Color online) Linear segmentation of detection efficiency. (a) Profile-based linear segmentation; (b) Measurement range division
pic

3. Minimum number of discrete segments for dynamic measurements. According to Eq. (14), the count rate within a characteristic energy window statistically follows a Gaussian distribution determined by the product of the radioactive intensity and detection efficiency. Even when flying over the same area at different speeds (i.e., with different measurement times but the same flight distance), the detection efficiency for a given ground region remains constant because it depends only on the geometry. Therefore, discretizing a dynamic measurement into multiple static measurements is primarily related to flight altitude and flight distance, rather than measurement duration. For example, if a flight distance of 10 m requires n discrete segments, then a flight distance of 20 m would require 2n segments by linear superposition. As shown in Figs. 6 and 9, the count rate of the spectrometer is mainly contributed by a small region directly beneath the detector, and this region is most sensitive to variations in flight distance. Thus, as shown in Fig. 11, this section analyzes the influence of the number of discrete segments n on detection efficiency within the left and right 50×50 m zones, under the condition of a GR820 spectrometer flying at 100 m altitude.

Fig. 11
(Color online) Schematic diagram of discrete segment number testing
pic

Based on the static measurement parameters derived in the previous two sections, the number of discrete segments was set from 1 to 20. The corresponding variations in the detection efficiency within the left- and right-side regions under different discretization numbers are shown in Fig. 12,.

Fig. 12
(Color online) Influence of discretization number on detection efficiency at different flight altitudes
pic

As shown in Fig. 12, when the discretization number is 1, the detection efficiency corresponds to the static measurement value at the detector position (25, 25, 100), resulting in a lower initial detection efficiency on the left side compared to the right side. As the discretization number increases, the detection efficiencies of the left and right regions gradually converge. Using a ±2% deviation as the convergence criterion, the minimum discretization number required for convergence at a flight altitude of 100 m was determined to be 4.

In summary, after obtaining the spatial distribution of radioactive finite elements within the measurement range as derived in Step 1, and developing a fast surface source detection efficiency calibration method based on linear segment midpoint approximation in Step 2, the construction method for the static measurement response matrix was established. Then, using the minimum discretization number for the dynamic measurement obtained in Step 3, the static response matrices are superimposed to construct the dynamic terrain response matrix. The complete workflow is illustrated in Fig. 13, and the relative error of each step remains within 2%, according to the principle of error propagation, the overall error is expected to remain within 3.5%.

Fig. 13
(Color online) Workflow diagram for constructing the dynamic terrain response matrix
pic
3

Experiment

3.1
Airborne and ground gamma spectrometry measurement

The experimental area was located in the Baiyun’ebo iron ore mining district, Darhan Muminggan United Banner, Baotou City, Inner Mongolia Autonomous Region. The site features abundant mineral samples and low vegetation cover, making it suitable for subsequent ground gamma-ray measurements and geological mapping to validate airborne survey data. The measurement instrument used was the GR820 airborne gamma spectrometer, and all operations strictly followed the “Specifications for Airborne Gamma-Ray Spectrometry” (EJ/T 1032-2005). The selected survey area spans from 110°00′11″E to 110°04′11″E and 41°37′29″N to 41°53′31″N, covering an area of approximately 165 km2. The layout of the survey lines is illustrated in Fig. 14. A total of 20 main survey lines were arranged longitudinally (east-west) with a spacing of approximately 350 m, and five transverse tie lines were arranged latitudinally (north-south) with an average spacing of 7.4 km, yielding a total of 10,085 measurement points. To ensure flight safety, the flight altitude was increased slightly. Approximately 81.40% of the survey was conducted at altitudes between 80 and 120 meters, with an average flight altitude of about 100 meters. The test zone included 10 main survey lines and 1 tie line, with a total of 987 measurement points.

Fig. 14
(Color online) Parameters of the aerial survey. (a) Aerial survey; (b) Terrain map of the survey area; (c) Flight altitude statistics chart
pic

To verify the effectiveness of the dynamic terrain correction algorithm, a test area was selected within the airborne survey region, as shown in Fig. 14(a) within the red box. This area is located near the Baiyun’ebo iron ore deposit, covering the coordinates 110°01′05″E to 110°03′04″E and 41°47′45″N to 41°50′50″N, with an approximate length of 5.7 km and a width of 2.7 km. Satellite imagery indicates that the terrain in this region is relatively flat, with clearly visible strip-shaped ore belts in the central and upper parts. The lower left corner of the area is adjacent to the Baiyun’ebo mining zone. During the survey, the ground was dry and the vegetation coverage was low, making the area well-suited for ground-based gamma spectrometry (the sampling locations are marked by red lines in the figure). The ground investigation section included 10 airborne gamma survey lines, comprising 2,083 measurement points. Measurements were conducted using the GR-320 ground gamma spectrometer, strictly following the “Technical Specifications for Ground Gamma-Ray Spectrometry” (DZ/T 0205–1999).

3.2
Evaluation and Analysis of the Method

Based on the spatial information provided by the airborne survey system, the dynamic terrain response matrix (with detection efficiency corrected to sensitivity) was constructed using the method described in the previous section. The PSO-MLEM inversion algorithm was then applied to the background-corrected airborne gamma data for eU, eTh, and K to estimate their surface concentrations. The matrix construction parameters are summarized in Table 1, and the inversion results are presented as contour maps in Figs. 15, 16, and 17.

Table 1
Inversion parameters for dynamic terrain correction
No. Parameter Indicator
1 Latitude × Longitude 500 × 20
2 Altitude Sub-Element 30 × 30 m
3 Measurement Range 5 × Flight Altitude
4 Linear Segmentation 12
5 Number of Discretizations 4
Show more
Fig. 15
(Color online) Comparison of eU content inversion results (a)Ground measurement results; (b) Altitude correction results (c) Dynamic terrain correction results
pic
Fig. 16
(Color online) Comparison of eTh content inversion results. (a) Ground measurement results; (b) Altitude correction results; (c) Dynamic terrain correction results
pic
Fig. 17
(Color online) Comparison of K content inversion results: (a) Ground measurement results; (b) Altitude correction results; (c) Dynamic terrain correction results
pic

As shown in Fig. 15, four anomalous zones were delineated based on the eU content from ground measurements. After applying conventional altitude correction to the airborne survey results, anomalous contours were observed in all zones except Zone 2. Although the eU content values were relatively low, the spatial variation pattern indicated that the airborne survey exhibited good accuracy. After the dynamic terrain correction, all anomalous zones became directly observable through contour maps. A comparison of the average eU content in each anomalous zone (Table 2) shows that the values obtained after dynamic terrain correction were generally closer to the ground measurement. Compared with altitude correction, the dynamic terrain-corrected inversion results showed an improvement in accuracy ranging from 3.12% to 25.79%.

Table 2
Comparison of eU content and relative error in anomaly zones
Anomalous Ground Altitude correction Terrain correction
zones (ppm) (ppm) Err.(%) (ppm) Err.(%)
1 11.27 5.56 -50.66 7.20 -36.10
2 8.87 6.35 -28.41 10.31 16.17
3 9.92 5.85 -40.97 8.41 -15.18
4 11.66 6.99 -40.06 7.36 -36.94
Show more

As shown in Fig. 16, five anomalous zones were delineated based on eTh content from ground measurements. After applying the altitude correction to the airborne survey data, anomalous contours were observed in all zones except Zones 1 and 4; however, the contour boundaries were relatively indistinct. Following the dynamic terrain correction, more distinct contours emerged in all zones, except Zone 4. Nevertheless, owing to the limited density of flight lines, Zone 1 exhibited a certain degree of spatial deviation. A comparison of the average Th content in each anomalous zone, as summarized in Table 3, indicates that the dynamic terrain correction improved the inversion accuracy by 7.50% to 48.82% compared with the altitude correction. Moreover, the dynamically corrected results show that the inversion accuracy for eTh content was significantly higher than that for eU. This can be attributed to the generally higher Th concentrations in the survey area and the stronger penetration capability of the characteristic gamma rays from the Th decay series, which effectively reduces the impact of statistical fluctuations on the inversion accuracy.

Table 3
Comparison of eTh content and relative error in anomaly zones
Anomalous Ground Altitude correction Terrain correction
zones (ppm) (ppm) Err.(%) (ppm) Err.(%)
1 44.86 20.79 -53.64 42.69 -4.82
2 42.68 40.69 -4.67 47.04 10.20
3 51.73 22.54 -56.43 40.91 -20.91
4 48.28 25.38 -47.44 37.04 -23.28
5 48.24 42.64 -11.63 46.28 -4.07
Show more

As shown in Fig. 17, five anomalous zones were delineated based on the K content from ground measurements. After applying altitude correction to the aerial survey results, only Zones 3 and 5 exhibited clear anomalies, whereas no obvious anomalies were observed in the other locations. After dynamic terrain correction, corresponding anomalies appeared in Zones 3, 4, and 5; however, positional deviations were observed in Zones 1 and 2, which are likely related to the sparser survey line layout of the aerial survey compared to ground measurement. A comparison of the average K content at each anomalous zone (Table 4) indicates that, compared with altitude correction, dynamic terrain correction improved the inversion accuracy by 9.10% to 36.52%.

Table 4
Comparison of K content and relative error in anomaly zones
Anomalous Ground Altitude correction Terrain correction
zones (%) (%) Err.(%) (%) Err.(%)
1 6.28 3.53 -43.78 5.80 -7.65
2 5.96 4.55 -23.62 6.41 7.68
3 5.63 5.06 -10.08 5.91 5.05
4 6.05 3.78 -37.49 4.34 -28.38
5 6.88 4.94 -28.16 5.81 -15.59
Show more

Based on the comprehensive analyses shown in Figs. 15, 16, 17 and Tables 2, 3, 4, compared with altitude correction, dynamic terrain correction produced U, Th, and K anomaly zone values that were closer to ground measurements, with accuracy improvements ranging from 3.12% to 48.82%, and also yielded clearer anomaly contours. However, owing to statistical fluctuations and the inherent nature of matrix-based inversion, each survey point is influenced by neighboring points within the measurement range. Elevated measurement values can absorb lower values, resulting in an apparent increase in the number of anomalous zones, as shown by the yellow blocks in Fig. 17(c). Nevertheless, the locations of the high-anomaly zones (represented by red and white blocks in the figure) generally aligned well with those from ground surveys.

It should also be noted that, apart from unavoidable measurement errors such as radar altimeter noise, GPS positioning inaccuracies, and variations in flight attitude, the aerial survey used fewer and non-corresponding survey lines compared to ground measurements. Slight positional deviations between the anomalous points on the dynamically terrain-corrected contour maps and those derived from ground measurements were also observed.

4

Summary

Based on a novel sourceless efficiency calibration method, this study derived the static measurement model of airborne gamma spectrometry into a dynamic measurement model and further improved the inversion accuracy of surface radionuclide concentrations by constructing and solving a dynamic terrain response matrix. It addresses the influence of the numerical calibration accuracy of the detection efficiency on the inversion results in traditional terrain correction methods, as well as the low simulation efficiency of Monte Carlo programs for airborne survey modeling. Airborne and ground gamma spectrometry experiments conducted in the Baiyun’ebo iron ore district revealed that: (1) after applying dynamic terrain correction, the anomalous zones of eU, eTh, and K concentrations increased significantly; using ground measurements as a reference, the corrected radionuclide concentrations showed higher accuracy than those obtained from conventional airborne inversion, with a maximum relative error reduction of 30.58%, and the anomaly boundaries became clearer and more consistent with the distribution observed in ground surveys; (2) due to statistical fluctuations in measurement and the inherent characteristics of matrix inversion, elevated values tend to absorb nearby lower values, resulting in the expansion of anomaly zones, although the locations of high-concentration anomalies remain largely consistent; (3) due to slight deviations between the airborne and ground survey line layouts, the positions of the inversion anomalies show minor discrepancies compared to the ground-measured anomaly zones.

References
1.H.S. Wu, Nuclear Technology Exploration. (Atomic Energy Press, Beijing, 1998), p.87. (in Chinese)
2.Y.X. Cheng, N.P. Wang, S.L. Hou, Nuclear Radiation Field and Radioactive Exploration. (Geology Press, Beijing, 2005), p. 84. (in Chinese)
3.M.H. Safaa, A.S.Y. Mohamed, S.G. Safwat,

Radioactive mineralization detection using remote sensing and airborne gamma-ray spectrometry at Wadi Al-Miyah area, Central Eastern Desert, Egypt

. Egypt. J. Remote Sens. Space Sci. 25 37-53 (2022). https://doi.org/10.1016/j.ejrs.2021.12.004
Baidu ScholarGoogle Scholar
4.EJ/T 1032–2005, Specification for airborne gamma-ray spectrometry, (China Commission of Science Technology and Industry for National Defense, Beijing, 2005). (in Chinese)
5.B. Minty, R. Brodie,

The 3D inversion of airborne gamma-ray spectrometric data

. Exploration Geophysics 47, 150157 (2016). https://doi.org/10.1071/EG14110
Baidu ScholarGoogle Scholar
6.Q.X. Zhang, Dissertation, Chengdu University of Technology, 2010. (in Chinese).
7.S.Q. Xiong,

Study on terrain effects and correction methods in the interpretation of airborne gamma spectrometric anomalies

. Earth Science Frontiers 5, 12 (1998). (in Chinese)
Baidu ScholarGoogle Scholar
8.J. H Wan, S.Q. Xiong, Z.G. Fan,

The status and prospects of airborne gamma-ray spectrometry technology and its application

. Geophysical and Geochemical Exploration 36, 386391 (2012). (in Chinese) https://doi.org/10.11720/wtyht.2012.3.12
Baidu ScholarGoogle Scholar
9.Y.F. Bai, Dissertation, Chengdu University of Technology, 2015. (in Chinese)
10.M.M. Hasan, J. Rutten, J. Camps et al.,

Minimum detectable activity concentration of radio-cesium by a LaBr3(Ce) detector for in situ measurements on the ground-surface and in boreholes

. Appl. Radiat. Isotop. 185, 110247 (2022). https://doi.org/10.1016/j.apradiso.2022.110247
Baidu ScholarGoogle Scholar
11.K. Sun, Dissertation, Chengdu University of Technology, 2022. (in Chinese)
12.Q.S. Liu, Dissertation, Chengdu University of Technology, 2022. (in Chinese)
13.K. Sun, C. Xiong, L.Q. Ge et al.,

Line unit correction method for airborne gamma-ray spectrum anomaly information based on wavelet transform

. Nucl. Tech. 41, 100501 (2018) (in Chinese). https://doi.org/10.11889/j.0253-3219.2018.hjs.41.100501
Baidu ScholarGoogle Scholar
14.S. Xu, Dissertation, Chengdu University of Technology, 2019. (in Chinese)
15.Z.N. Tian, Dissertation, Tsinghua University, 2015. (in Chinese)
16.H.X. Wu, Dissertation, Chengdu University of Technology, 2016. (in Chinese)
17.C. Xiong, X.Y. Yuan, X. Wan et al.,

Study on sourceless efficiency calibration scheme based on numerical integration method

. Nuclear Electronics and Detection Technology 44, 146152 (2024). (in Chinese). https://doi.org/10.3969/j.issn.0258-0934.2024.01.020
Baidu ScholarGoogle Scholar
18.J. Zhang, Q.X. Zhang, W.C. Lai et al.,

The full-energy peak efficiency of gamma-ray spectrometry by numerical calculation method

. Nucl. Instrum. Meth. Phys. Res. Sect. A 1061, 169145 (2024). https://doi.org/10.1016/j.nima.2024.169145
Baidu ScholarGoogle Scholar
19.H.F. Liu, L.Q. Ge, X.C. Xie et al.,

Calculation of desired X-ray collection angle on XRF analyzer designed by Monte Carlo method

. Nucl. Tech. 38 060502 (2015). (in Chinese) https://doi.org/10.11889/j.0253-3219.2015.hjs.38.060502
Baidu ScholarGoogle Scholar
20.K. Noack,

Efficiency and reliability in deep-penetration Monte-Carlo calculations

. Ann. Nucl. Energy 16, 309316 (1991). https://doi.org/10.1016/0306-4549(91)90016-Q
Baidu ScholarGoogle Scholar
21.J. Ghassoun, A. Jehouani,

Russian roulette efficiency in Monte Carlo resonant absorption calculations

. Appl. Radiat. Isotop. 53, 881885 (2000). https://doi.org/10.1016/S0969-8043(00)00241-4
Baidu ScholarGoogle Scholar
22.J. Zhao, T.C. Feng,

Monte Carlo simulation of single-energy gamma radionuclide depth distribution for in-situ measurement

. Paper presented at Proceedings of the National Academic Symposium on Nuclear and Radiation Facility Decommissioning, Chinese Nuclear Society Radiation Protection Branch, Mianyang, 1 Nov 2007. (in Chinese)
Baidu ScholarGoogle Scholar
23.B.L. Yang, Q. Zhou, J. Zhang et al.,

Performances of different efficiency calibration methods of high-purity-germanium gamma-ray spectrometry in an inter-comparison exercise

. Nucl. Sci. Tech. 30, 37 (2019) https://doi.org/10.1007/s41365-019-0562-x
Baidu ScholarGoogle Scholar
24.T. Frosio, N. Menaa, P. Bertreix et al.,

A novel technique for the optimization and reduction of gamma spectroscopy geometry uncertainties

. Appl. Radiat. Isotop. 156, 108953 (2020). https://doi.org/10.1016/j.apradiso.2019.108953
Baidu ScholarGoogle Scholar
25.Q.X. Zhang, Y.L. Guo, S. Xu et al.,

A hybrid method on sourceless sensitivity calculation for airborne gamma-ray spectrometer

. Appl. Radiat. Isotop. 137, 6872 (2018). https://doi.org/10.1016/j.apradiso.2018.03.009
Baidu ScholarGoogle Scholar
26.W.C. Li, Q. X Zhang, H.X. Wu et al.,

A gamma ray sourceless efficiency calibration method based on the Boolean operation of the ray deposition process

. Appl. Radiat. Isotop. 193, 110620 (2023). https://doi.org/10.1016/j.apradiso.2022.110620
Baidu ScholarGoogle Scholar
27.Y. Gu, Dissertation (Chengdu University of Technology, 2010). (in Chinese)
28.Y.H. Li, R. Yang, Q.X. Zhang et al.,

Study on PSO and MLEM hybrid algorithm for NDP spectrum unfolding

. Atomic Energy Science and Technology 58, 11521159 (2024). (in Chinese) https://doi.org/10.7538/yzk.2023.youxian.0642
Baidu ScholarGoogle Scholar
Footnote

The authors declare that they have no competing interests.