Method for rapid survey of natural radioactivity environment and target area delineation based on satellite remote sensing

Abstract

The application discloses a satellite remote sensing based natural radioactivity environment rapid survey and target area delineation method, and relates to the technical field of geophysical prospecting. The method comprises the following steps: acquiring multi-source satellite remote sensing data and historical measured label data of a target area; extracting a radioactivity correlation feature set; training the radioactivity correlation feature set to a subarea mapping model of a natural radioactivity level according to a geological unit; performing global inversion to generate a radioactivity level prediction surface and a prediction uncertainty surface; determining a region with a prediction value exceeding a preset risk threshold as a high-value region, and delineating a part of the high-value region with an uncertainty lower than a preset confidence threshold as a determined target area, and the remaining part as a to-be-checked region; incorporating detailed survey data of the to-be-checked region into the label and incrementally updating the subarea mapping model, and re-inverting and delineating based on the subarea mapping model. The application has the advantages that a globally classifiable radioactivity level distribution can be given without adding large-scale measurement, and detailed survey forces can be directed to an untrustworthy region.

Description

Technical Field

[0001] This application relates to the field of geophysical exploration technology, and in particular to a method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing. Background Technology

[0002] Natural radioactivity originates from nuclides such as uranium, thorium, and potassium in rocks and soil. Its content is controlled by the lithology and tectonic development of bedrock, and is higher in areas with granite and alkaline rocks than in surrounding areas. People living in areas with high levels of natural radioactivity for extended periods will receive additional radiation through indoor radon inhalation and external radiation from building materials. Therefore, it is essential to first understand the spatial distribution of this level in national radiation environmental baseline surveys, residential area radon risk screening, and building stone source area investigations.

[0003] These types of surveys cover continuous land areas at the county level or even larger. However, ground-based methods that can provide traceable quantitative values ​​are limited by manpower and time constraints, and a single operation can only cover a small portion. Therefore, in practice, a low-cost, rapid census of the entire area is conducted first to delineate key areas, and then limited detailed survey resources are allocated to these key areas. The viability of this approach depends on whether the rapid census can provide a horizontal distribution that covers the entire area and can be used for hierarchical ranking, and whether it can indicate which conclusions need to be verified on the ground.

[0004] Existing rapid survey methods fall into two categories. Ground-based point-by-point surveys, vehicle-mounted mobile surveys, and airborne gamma spectroscopy measurements rely on detectors to receive surface gamma photons and invert nuclide content according to energy windows, providing traceable values. However, the detectors must be transported to the vicinity of the measured location to obtain readings, thus binding coverage to transportation costs. Ground-based surveys have a single-point field of view on the order of ten meters, vehicle-mounted mobile surveys are restricted by road networks, and airborne surveys have a much larger spacing between survey lines than the sampling interval along the survey lines. The results all depend on spatial interpolation between measurement points, and since radioactivity levels exhibit a step change in the lithological contact zone, gradual interpolation would inevitably smooth out this zone. Remote sensing alteration anomaly extraction does not require transporting detectors, reducing coverage costs by several orders of magnitude. However, the output is a dimensionless relative intensity, and the responses to topographic shadows and surface humidity are on the same order of magnitude as actual alteration. As vegetation cover increases, the surface signal is suppressed by the canopy, and the relationship between alteration and nuclide enrichment is statistically correlated rather than causally equivalent, so the output can only serve as a clue. The two categories of coverage and quantification are at opposite ends of the spectrum. To obtain traceable values, detectors must be sent to the field and accept the resulting reduced coverage. To achieve full coverage, dimensions and traceability must be abandoned. Therefore, it is impossible to provide a horizontal distribution that covers the entire area, can be graded and ranked, and indicates where the conclusions are credible. Ground-based detailed survey forces have no way of determining where to direct their efforts. Summary of the Invention

[0005] In order to provide a comprehensive, hierarchical, and reliable distribution of natural radioactivity levels without the need to move detectors point by point, this application provides a method for rapid survey and target delineation of natural radioactivity environments based on satellite remote sensing.

[0006] This application provides a method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing, employing the following technical solution: A rapid survey and target delineation method for natural radioactive environments based on satellite remote sensing, comprising the following steps: S1. Acquire multi-source satellite remote sensing data covering the target area, as well as historical measured tag data carrying the measured values ​​of natural radioactivity levels in a portion of the target area; S2. Extract a radioactive correlation feature set from the multi-source satellite remote sensing data, wherein the radioactive correlation feature set includes multiple proxy features characterizing the geological background of natural radionuclide enrichment; S3. Divide the target area into multiple geological units, and train a partitioning mapping model from the radioactivity association feature set to the natural radioactivity level in each geological unit using the historical measured label data; S4. Perform a global inversion on the target region using the partition mapping models described above to generate a radioactivity level prediction surface and a prediction uncertainty surface; S5. The region where the value of the radioactivity level prediction surface exceeds the preset risk threshold is identified as a high-value region. The part of the high-value region where the value of the prediction uncertainty surface is lower than the preset confidence threshold is identified as the target area. The remaining part of the high-value region is identified as the area to be verified. S6. The detailed survey data obtained by performing a detailed ground survey on the area to be verified is added as a new label and incorporated into the historical measured label data. The zoning mapping model of the geological unit to which the area to be verified belongs is incrementally updated. Based on the zoning mapping model after the incremental update, S4 and S5 are re-executed, and the determined target area is output.

[0007] By adopting the above technical solution, gamma rays cannot penetrate the entire atmosphere, and satellites cannot directly read the radioactivity values ​​on the Earth's surface. Therefore, observable measurements such as alteration, tectonics, and topography are used to construct a set of radioactive correlation features. Historical measured values ​​from some areas are used to link proxy features with the true level. The values ​​originate from ground measurements and the coverage is expanded to all pixels. Regional training ensures that non-stationary relationships controlled by lithology are established within each unit. The prediction uncertainty surface and the radioactivity level prediction surface participate in the segmentation in parallel, and areas where the model is unreliable are listed separately as areas to be verified. The feedback of measured values ​​from detailed surveys drives incremental updates with new labels. Each round of detailed surveys narrows the unreliable range once, and the general survey and detailed surveys form a closed loop.

[0008] Optionally, in step S5, the predicted uncertainty surface is further decomposed into a cognitive uncertainty component and a random uncertainty component. The cognitive uncertainty component represents the uncertainty that can be converged through supplementary observations due to label scarcity, and the random uncertainty component represents the uncertainty that cannot be converged through supplementary observations due to label noise propagation. For the portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, supplementary observation suggestions are output. The supplementary observation suggestions include suggested observation locations, suggested observation types, and observation priorities. For the portion of the area to be verified where the random uncertainty component is greater than the cognitive uncertainty component, data quality limitation labels are output.

[0009] By adopting the above technical solution, uncertainty is decomposed into two components based on whether it can converge through supplementary observations. Supplementary observation suggestions are given at the point where the cognitive uncertainty component is dominant, while the data quality is limited at the point where the accidental uncertainty component is dominant. This distinguishes between those that are worth investigating in detail and those that are difficult to converge even with investigation.

[0010] Optionally, in S6, the detailed survey measured data is the data obtained by performing a detailed ground survey on the portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, according to the supplementary observation suggestion; after each round of re-executing S4 and S5, it is evaluated whether the preset convergence condition is met: if there is no portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component in this round, it is determined that the preset convergence condition is met; if there is a portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component in this round, this round is considered to be... The area of ​​the region to be verified is deemed to meet the preset convergence condition if the decrease in area of ​​the region to be verified relative to the previous round and the decrease in the total amount of cognitive uncertainty components in the current round relative to the total amount of cognitive uncertainty components in the previous round are both lower than the preset convergence limit. The first round is deemed not to meet the preset convergence condition. When the preset convergence condition is met, the next round of incremental update is terminated, and an undecidable region list is output. The undecidable region list includes the remaining regions to be verified at this time, and the proportion of the random uncertainty components in the regions to be verified in the value of the predicted uncertainty surface.

[0011] By adopting the above technical solution, the closed loop obtains the termination condition. When the part dominated by the cognitive uncertainty component disappears or the decrease of both components is lower than the preset convergence limit, the iteration stops. The remaining undetermined part is explicitly transferred in the list of undeterminable areas, avoiding endless detailed investigation.

[0012] Optionally, the historical measured tag data includes historical airborne gamma spectrum data, and S3 includes sub-steps S31-S34: S31. Perform a label timeliness check on the historical measured label data. The label timeliness check is based on the multi-temporal multi-source satellite remote sensing data to perform surface disturbance detection on each label location, and downweight or remove the labels that detect surface disturbance to obtain the checked labels. S32. The radioactive correlation feature set is weighted and aggregated to the detection field of view (PDV) scale according to the detection field of view (PDV) response function of each measurement point in the historical airborne gamma spectrum data to obtain the patch scale features. The effective width of the PPV response function is determined point by point according to the flight altitude record and the terrain undulation of the location of each measurement point. S33. Sampling is performed on the verified labels according to the characteristic distribution quota within each geological unit to correct the survey corridor bias of the historical airborne gamma spectral data and obtain a training sample set; S34. Based on the training sample set and the patch scale features, train the partition mapping model at the patch scale of the detection field of view; In step S4, the original pixel resolution of the multi-source satellite remote sensing data is used to apply each of the partition mapping models.

[0013] By adopting the above technical solution, historical airborne gamma spectrum data is converted into usable weak area labels after tag timeliness verification, detection field spot aggregation and offset correction. Training is completed at the detection field spot scale and applied back to the original pixel resolution. Scale misalignment does not enter the radioactivity level prediction surface.

[0014] Optionally, the historical measured tag data also includes ground measured point data and geochemical unmixed tags. The geochemical unmixed tags are obtained by unmixing and allocating the elemental content data of stream sediments through the catchment area. In step S34, the training sample set is also subjected to weighted training according to the noise variance of each tag. The noise variance corresponding to the historical airborne gamma spectrum data is calculated by the Poisson statistical characteristics of gamma counts, the flight altitude attenuation factor, and the action width of the detection field of view spot response function. The noise variance corresponding to the ground measured point data is determined by the accuracy index of the measuring instrument. The noise variance corresponding to the geochemical unmixed tags is determined by the residual of the unmixing allocation in the catchment area. The noise variance of each tag is unified to the same content dimension before participating in the weighted training.

[0015] By adopting the above technical solution, the noise variance of three types of labels with different sources and different accuracies is calculated according to their respective measurement physics and then weighted and trained after unifying the dimensions. Labels with low accuracy are no longer weighted equally with labels with high accuracy to skew the model.

[0016] Optionally, in step S2, feature-level multi-temporal robust synthesis is further performed on the radioactive correlation feature set. The feature-level multi-temporal robust synthesis includes: taking robust statistics as feature values ​​for the multi-temporal sequences of each proxy feature, outputting the temporal variation coefficient of each pixel as the feature quality weight, and performing periodic component decomposition on the multi-temporal sequences to distinguish between periodic phenological changes and abrupt surface disturbances; and providing the detection results of the abrupt surface disturbances to the surface disturbance detection in step S31.

[0017] By adopting the above technical solution, the proxy features complete multi-temporal robust synthesis at the feature layer, the temporal variation coefficient is converted into feature quality weights to continue to participate in modeling, and the periodic component decomposition distinguishes phenological fluctuations from surface abrupt changes and provides them back to the label timeliness verification.

[0018] Optionally, S2 includes sub-steps S21-S24: S21. Based on the comparison result between vegetation coverage and a preset coverage threshold, the pixels of the target area are divided into low-coverage pixels and high-coverage pixels; S22. Extract the radioactive correlation feature set from the low-coverage pixels through the main spectral channel, wherein the main spectral channel is based on the alteration spectral features of the short-wave infrared band; S23. Extract the radiometric correlation feature set from the high-coverage pixels through an alternative channel, wherein the alternative channel is based on at least one of synthetic aperture radar backscattering features, topographic structure features, and neighborhood lithological inference features; S24. Select pixels from the low-coverage pixels where both the main spectral channel and the alternative channel are available, and where the terrain shadow and surface humidity do not exceed the preset interference limit, as bridging samples. Extract the radiometric correlation feature set from the bridging samples through the alternative channel. Fit the mapping relationship between the output of the alternative channel and the output of the main spectral channel using the bridging samples to complete the bridging calibration, so that the output of the alternative channel and the output of the main spectral channel have the same dimensions.

[0019] By adopting the above technical solution, the high vegetation cover pixels are modified to extract features through alternative channels, and bridging calibration is completed with overlapping pixels that are available in both channels and have limited interference. The outputs of the two channels fall on the same dimension, and the vegetation area is no longer removed by the entire mask.

[0020] Optionally, in step S24, the fitting residuals of the bridging calibration on each bridging sample are statistically analyzed, and a channel difference uncertainty model that monotonically increases with vegetation coverage is established for the fitting residuals according to the vegetation coverage of each bridging sample; in step S4, the channel difference uncertainty components are obtained by taking values ​​of the channel difference uncertainty model according to the vegetation coverage of each pixel, and the channel difference uncertainty components are injected into the prediction uncertainty surface.

[0021] By adopting the above technical solution, the fitting residual of the bridging calibration itself is used to build a monotonically increasing model based on vegetation coverage and inject the prediction uncertainty surface. The unreliability of the location with a heavier share of the alternative channel inference is taken into account in the diversion process.

[0022] Optionally, after step S5, the following source type determination operation is further performed on each patch contained in the determined target region: The extrapolated proportion is obtained by statistically analyzing the proportion of pixels whose values ​​of the radioactive correlation feature set within each of the aforementioned patches fall outside the training feature distribution of the partitioning mapping model. Detect the occurrence time and range evolution information of each patch in the multi-temporal multi-source satellite remote sensing data; The patches whose extrapolation ratio exceeds the preset extrapolation limit and whose occurrence time is detected are identified as human sources, classified into human radioactive source target areas, and source term investigation suggestions are output. The source term investigation suggestions include the occurrence time and the range evolution information. The patches whose extrapolation ratio does not exceed the preset extrapolation limit and whose occurrence time is not detected are identified as geological sources, included in the geological target area, and ground detailed investigation suggestions are output.

[0023] By adopting the above technical solution, the source type is distinguished by two criteria: the appearance time of the map patch and the proportion of model extrapolation. The target area of ​​man-made stockpile is transferred to the source term investigation, and the target area of ​​geological origin is transferred to the detailed ground investigation. The two types of subsequent procedures are no longer mixed.

[0024] Optionally, between S3 and S4, a transition zone is further included on both sides of the contact zone between two adjacent geological units; in S4, for each pixel located in the transition zone, the prediction results of the respective partition mapping models of the two adjacent geological units are weighted and fused according to the membership degree of each pixel to the two adjacent geological units; the width of the transition zone is set differently according to the geological type of the contact zone, and the width of the transition zone corresponding to fault contact is smaller than the width of the transition zone corresponding to intrusive gradient contact.

[0025] By adopting the above technical solution, the model output of adjacent units on both sides of the contact zone is fused with membership degree weighting, the pseudo-jump caused by the hard boundary of the partition is erased, and the width of the transition zone varies with the contact type so that the real step at the fault is preserved.

[0026] Optionally, in step S5, delineation is performed on a per-patch basis, a minimum area constraint is applied to each patch, the consistency between the principal axis of the shape of each patch and the orientation of the local structure at its location is statistically analyzed, and the target area priority is reduced for patches whose consistency is lower than a preset consistency limit; when performing morphological operations, the orientation of the structural elements is adaptively determined according to the orientation of the local structure.

[0027] By adopting the above technical solution, the delineation results are extended along the structural direction, and isolated fragments generated by pixel-by-pixel threshold segmentation are suppressed.

[0028] Optionally, the unmixing assignment of the catchment area includes: extracting the catchment area of ​​each river sediment sampling point based on the digital elevation model, determining the assignment weight according to the area, lithology and erosion intensity of each pixel in the catchment area, and assigning the element content data of the river sediment to each pixel in the catchment area using regularized inversion to obtain the geochemical unmixing label.

[0029] By adopting the above technical solution, the mixed readings at the outlet of the catchment area are restored to the shared values ​​of each pixel upstream, and the existing geochemical results are converted into weak area labels.

[0030] Optionally, in S3, for the geological unit whose number of historical measured label data is lower than a preset sample size threshold, the global inversion is performed using the partition mapping model of another geological unit adjacent to the geological attribute of the geological unit, and the value of the geological unit in the prediction uncertainty surface is improved.

[0031] By adopting the above technical solution, the geological units without labels can still obtain predicted values, and these values ​​are simultaneously marked as low confidence and diverted to the area to be verified.

[0032] Optionally, each of the partition mapping models adopts a hierarchical structure and shares the same set of model parameter priors. The geological units with a smaller amount of historical measured label data are trained with the partition mapping model according to the constraints of the model parameter priors.

[0033] By adopting the above technical solution, small sample geological units can complete training by sharing prior parameters, thereby reducing overfitting while preserving the differences between regions.

[0034] Optionally, in step S4, feature space extrapolation detection is also performed on each pixel, and pixels whose values ​​of the radioactive association feature set fall outside the training feature distribution of the partition mapping model are included in the area to be checked.

[0035] By adopting the above technical solution, predictions made by the model outside the training distribution are excluded from the target area, and the failure of the extrapolation position becomes visible.

[0036] Optionally, the proxy features include alteration spectral features, tectonic linear volume density features, tectonic linear volume convergence features, topographic factor features, and thermal inertia features.

[0037] By adopting the above technical solutions, the proxy characteristics are respectively pointed to the nuclide enrichment background from mineral alteration, tectonic flow channels and surface material composition.

[0038] Optionally, in S3, the geological unit is obtained by fusing regional geological maps with remote sensing lithological classification results.

[0039] By adopting the above technical solution, the zoning boundary not only inherits the existing geological mapping results, but also supplements the insufficient accuracy of the remote sensing lithology classification.

[0040] Optionally, after step S6, a comprehensive score is performed on each patch contained in the determined target area based on the value of the predicted surface of the radioactivity level, the patch area, the population exposure weight, and the accessibility of the transportation, and a detailed investigation priority list is output according to the comprehensive score.

[0041] By adopting the above technical solutions, the target area was transformed from a parallel list into a sequential work list, and detailed investigations were carried out step by step according to risk and input-output ratio.

[0042] In summary, this application includes at least one of the following beneficial technical effects: 1. Replace the gamma readings that cannot be directly measured by satellite with satellite proxy features, establish a mapping within the geological unit and perform full-domain inversion by using sparse historical labels, and divert the target area by dual constraints of the radioactivity level prediction surface and the prediction uncertainty surface, so as to achieve both full coverage and scalable values.

[0043] 2. The detailed investigation and measurement data feedback is used to drive incremental updates of the partition mapping model with new labels. Combined with uncertainty decomposition, convergence criteria and undecidable area list, the general survey and detailed investigation form a closed loop of convergence in each round, and the unconverged part is explicitly transferred.

[0044] 3. Label timeliness verification, detection field spot response function aggregation, noise variance weighting and channel bridging calibration make existing historical data usable without adding new measurements, and the deviations caused by vegetation shading and human source mixing are included in the visible uncertainty. Attached Figure Description

[0045] Figure 1 This is a schematic diagram of the overall process of a rapid survey and target area delineation method for natural radioactive environments based on satellite remote sensing in one embodiment of this application.

[0046] Figure 2 As shown in one embodiment of this application Figure 1 A detailed flowchart of step S2.

[0047] Figure 3 As shown in one embodiment of this application Figure 1 A detailed flowchart of step S3.

[0048] Figure 4 As shown in one embodiment of this application Figure 1 A schematic diagram of the decision-making and diversion process in step S5. Detailed Implementation

[0049] The present application will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative of the application and are not intended to limit the scope of the application.

[0050] This application discloses a method for rapid surveying and target area delineation of natural radioactivity environments based on satellite remote sensing. This application uses observable satellite surrogate features as substitutes for readings of natural radioactivity levels. It links these surrogate features to the true levels using existing sparse measured values ​​within the target area, and outputs the reliability of the inference result along with the result itself. First, multi-source satellite remote sensing data and historical measured label data are acquired, and a set of radioactivity-related features is extracted. Then, a regional mapping model is trained according to geological units, and a global inversion is performed on the target area to obtain a radioactivity level prediction surface and a prediction uncertainty surface. Finally, the two surfaces are used together to constrain the delineation, separating high-value areas into target areas and areas to be verified. Detailed measured data from the areas to be verified are fed back as new labels to drive incremental updates of the model. Thus, a comprehensive and hierarchically sortable distribution of natural radioactivity levels can be obtained without transporting detectors point by point, allowing ground-based detailed surveys to be directed to locations where conclusions are not yet credible.

[0051] The following section provides a detailed explanation of the terminology, application environment, overall method flow, specific implementation of each step, and device mapping.

[0052] The terms "first," "second," etc., in the specification and claims of this application are used to distinguish similar objects and do not describe a specific order; "comprising" and "having" are intended to cover non-exclusive inclusion.

[0053] Natural radioactivity level refers to the content of naturally occurring radioactive nuclides such as uranium, thorium, and potassium per unit mass or volume of surface media, or the absorbed dose rate of gamma radiation produced by the decay of these nuclides in near-surface air. This quantity is determined by the material composition of rocks and soils, and is spatially distributed in a blocky manner depending on the bedrock lithology and structural development, with magnitude differences possible on both sides of the lithological contact zone. This quantity can be directly measured by ground-based gamma spectrometers, dose rate meters, or large-volume scintillation crystal detectors flying at low altitudes. However, gamma rays decay exponentially in air, and after passing through the entire atmosphere, they are insufficient to form a recognizable signal at orbital altitude, thus making direct reading by spaceborne detectors impossible. Therefore, this method does not include a direct spaceborne measurement step.

[0054] A surrogate feature is an observable measurement that can be directly calculated from multi-source satellite remote sensing data and has a stable correlation with the geological background of natural radionuclide enrichment. The surrogate feature itself is not a radioactivity value; its function is to act as a substitute for the reading of natural radioactivity levels, converting quantities that are not observable from space into observable quantities. Multiple surrogate features organized according to a unified spatial grid constitute a radioactivity-correlated feature set. Each pixel in the radioactivity-correlated feature set corresponds to a set of feature values, and these feature values ​​collectively characterize the geological background conditions of radionuclide enrichment in that pixel.

[0055] Historical measured label data refers to existing data within a target area that contains measured values ​​of natural radioactivity levels for a portion of the area. These values ​​originate from direct measurements taken at ground level or low altitude, thus allowing for traceability to the measuring instruments. This type of data is spatially sparse and unevenly distributed, with a coverage rate lower than that of full satellite imagery, but its quantity is sufficient to support the establishment of a statistical relationship between proxy characteristics and natural radioactivity levels. The radioactive geological background is quasi-static on a geological timescale; measured values ​​collected many years ago remain usable in undisturbed locations on the surface. This allows existing data to be reused as labels without the need for remeasurement.

[0056] A geological unit refers to a continuous block within a target area whose lithological composition and diagenetic background are consistent. The relationship between proxy characteristics and natural radioactivity levels within the same geological unit can be considered stable, while this relationship can differ significantly between different geological units. The boundaries of geological units are usually distributed along lithological contact zones, fault zones, or unconformities.

[0057] A zonal mapping model is a regression model trained on historically measured labeled data covered by a single geological unit, derived from a set of radioactivity-related features to the level of natural radioactivity. Each geological unit has its own zonal mapping model, with parameters independent of each other. The differences in lithologically controlled relationships mean that these models are not compressed into the same set of parameters. The model's input is a set of feature values ​​for a pixel, and its output is the inferred level of natural radioactivity for that pixel.

[0058] The radioactivity level prediction surface refers to a continuous surface of values ​​obtained by inferring values ​​from each pixel of the target area using the mapping model for each region. Its spatial resolution is consistent with the original pixel resolution of the multi-source satellite remote sensing data used. The prediction uncertainty surface is another surface of values ​​with the same grid and resolution as the radioactivity level prediction surface. The value of each pixel in the uncertain surface represents the degree of unreliability of the inferred value at that pixel. The two surfaces are output side by side, so that the inferred values ​​and their reliability levels are readable simultaneously on the same grid.

[0059] The target area is defined as the region where the predicted radioactivity level exceeds a preset risk threshold, but the predicted uncertainty level is below a preset confidence threshold. It indicates a location that both indicates a high level of natural radioactivity and whose inferences are credible. The area to be verified is defined as the region where the predicted radioactivity level exceeds a preset risk threshold, but the predicted uncertainty level is not below a preset confidence threshold. It indicates a high level of natural radioactivity, but whose inferences are not yet directly credible. These two types of areas together constitute a diversion of the high-value area.

[0060] Detailed ground survey refers to the process of dispatching personnel with gamma spectrometers or dose rate meters to the site to conduct point-by-point measurements according to the aforementioned area list. It has a small spatial coverage area but reliable single-point measurements. The measured data from the detailed ground survey, i.e., the measured values ​​of natural radioactivity levels produced by the ground survey, plays two roles in this method: it is both the final determination of the natural radioactivity level at that location and a new label that can be fed back into the model.

[0061] The application environment of this method is described below. This method can be applied to scenarios that require first understanding the distribution of natural radioactivity levels over a large area before arranging ground operations, such as national radiation environment baseline surveys, radon risk screening in residential areas, investigation of building stone source areas, and health risk early warning in areas with high geological background. Taking a county with both granite and Quaternary overburden as the target area, this county contains both biotite granite outcrops and thick Quaternary overburden areas and metamorphic sandstone distribution areas, with a total area larger than the area that a detailed ground survey can cover in one operation cycle. This method is at the forefront of this type of survey process chain. Its output area list serves as the basis for subsequent detailed ground surveys, while the output of the detailed ground surveys flows back to this method, forming a sequential process pair.

[0062] The following is combined Figure 1 The overall process of this method is described in detail. This method includes steps S1 to S6, which are executed sequentially.

[0063] S1. Acquire multi-source satellite remote sensing data covering the target area, as well as historical measured label data carrying the measured values ​​of natural radioactivity levels in some areas within the target area.

[0064] In practice, multi-source satellite remote sensing data includes multispectral and hyperspectral optical imagery, synthetic aperture radar imagery, thermal infrared imagery, and digital elevation models. These data types provide observations on surface mineral composition, surface roughness and dielectric properties, diurnal temperature variation, and topographic features. Historical measured data consists of three sources, each differing in spatial representation, dimensions, and precision.

[0065] Tag point 1 is a measurement point for historical airborne gamma-ray spectroscopy measurements. These measurement points are used by helicopters or fixed-wing aircraft carrying large-volume scintillation crystal detectors, collecting data along the measurement line. A single reading covers a circular patch on the ground with a diameter on the same order of magnitude as the flight altitude, and the value is given as nuclide content or equivalent content. This type of data is dense along the measurement line and sparse perpendicular to the line, exhibiting a striped distribution within the target area.

[0066] Tag point 2 is a ground-based measurement point. These types of measurement points are obtained point by point by operators using instruments. Each reading covers a range of about ten meters around the measurement point, and the measurement value can be directly traced back to the instrument's accuracy specifications. This type of data point is few in number and scattered in location, mostly concentrated near existing engineering sites and known anomalies.

[0067] Point 3 is a geochemical sampling point for river sediments. These sampling points are located within river valleys, and their elemental content readings reflect the results of erosion, transportation, and mixing of materials throughout the upstream catchment area, and therefore do not directly correspond to the surface content at the sampling point's location. This type of data has accumulated significantly in regional geological surveys and covers a wide area.

[0068] After acquiring the above data, before proceeding with feature extraction, the spatiotemporal reference of the multi-source satellite remote sensing data is unified. This involves resampling the images from each source to the same spatial grid and registering them to the same coordinate system, while also assigning the observation times from each source to a unified time axis. For optical images, they are synthesized pixel-by-pixel based on cloud cover, observation angle, and image quality within a certain time window to obtain a composite image with fewer clouds. For areas with perennial clouds and rain, locations with insufficient effective optical observations are supplemented by concurrent synthetic aperture radar (SAR) images to ensure that spatial coverage is not disrupted by cloud distribution. After these preparations, the multi-source satellite remote sensing data constitutes a grid-aligned, temporally comparable set of observations across the entire domain, which is used by S2 to extract surrogate features.

[0069] S2. Extract a radioactive correlation feature set from multi-source satellite remote sensing data. The radioactive correlation feature set includes multiple proxy features that characterize the geological background of natural radionuclide enrichment.

[0070] In practical implementation, feature extraction focuses not on finding the spectral response of radioactivity itself, but on finding observable traces left on the surface by the geological conditions underlying nuclide enrichment. Uranium and thorium migrate with fluids and accumulate in specific tectonic locations during diagenesis and subsequent hydrothermal activity. This process simultaneously alters the mineral composition of the surrounding rocks and changes the degree of surface fragmentation and thermophysical properties. Therefore, the surface mineral alteration, the density of the tectonic framework, and differences in thermal inertia can all serve as indirect indicators of the nuclide enrichment background. The radioactive correlation feature set is thus organized into a unified grid based on these indirect indicators.

[0071] In some embodiments, the proxy features include alteration spectral features, tectonic linear volume density features, tectonic linear volume convergence features, topographic factor features, and thermal inertia features.

[0072] Alteration spectral characteristics, calculated from multispectral or hyperspectral images in the visible-near-infrared to short-wave infrared bands, are used to characterize the relative enrichment of aluminum-containing hydroxyl minerals and iron-stained minerals. These minerals are products of hydrothermal fluid activity, and their spatial distribution coincides with the fluid migration path, thus being correlated with the precipitation location of nuclides carried by the fluid. For example, the diagnostic absorption depth near 2.2 micrometers in the short-wave infrared is taken as the hydroxyl alteration intensity, and the ratio of the red to blue bands in the visible light is taken as the iron staining intensity; these two together constitute the two components of the alteration spectral characteristics.

[0073] The linear volume density feature is constructed by statistically analyzing linear structural traces extracted from the image's edge response and topographic shadows using a digital elevation model. The value is the cumulative length of the linear volume traces per unit area. Linear structures correspond to fracture and joint development zones, serving as channels for the ascent of deep fluids; the denser the channels, the larger the contact area between the fluid and the surrounding rock. This feature is statistically output on a grid using a moving window, with the window size matching the minimum discernible area of ​​the target region.

[0074] The intersection degree characteristic of linear bodies is constructed by statistically analyzing the number of intersecting nodes of linear bodies with different orientations per unit area. Unlike density characteristics, which only reflect the number of channels, intersection degree reflects the connectivity and superposition relationship between channels. At intersection nodes, the fluid residence time is longer, and the conditions for nuclide precipitation are more sufficient. Locations with high density and low intersection degree are not equivalent to locations with both high density and high intersection degree in terms of mineralization indication; therefore, the two characteristics are listed separately rather than combined.

[0075] Topographic features, derived from the digital elevation model, include slope, aspect, topographic relief, and curvature. These parameters characterize the degree of erosion or overburdening of surface materials, affecting whether bedrock signals can be observed by satellites. They also participate in subsequent discrimination of topographic shadows and illumination differences, identifying the portion of proxy feature values ​​caused by topography rather than geology.

[0076] Thermal inertia characteristics, calculated by combining diurnal temperature variation and albedo from thermal infrared images, reflect the ability of the surface medium to resist temperature changes. The thermal inertia of dense crystalline rock masses differs significantly from that of loose, weathered overburden, thus providing a supplementary criterion for lithological tendency in locations where spectral information is limited.

[0077] The aforementioned proxy features are calculated pixel by pixel on the global grid. The multiple feature values ​​of each pixel constitute the feature vector of that pixel. The set of feature vectors of all pixels is the radial correlation feature set, which is used for both S3 training and S4 inversion.

[0078] Vegetation canopy exhibits strong moisture absorption in the shortwave infrared band. In areas with high canopy coverage, surface signals are suppressed by the canopy, and the alteration spectral characteristics no longer reflect the surface mineral composition. To address this situation, refer to... Figure 2 In some embodiments, S2 includes sub-steps S21 to S24.

[0079] S21. Based on the comparison between vegetation cover and a preset cover threshold, the pixels in the target area are divided into low-coverage pixels and high-coverage pixels. Vegetation cover is calculated from the normalized vegetation index using a pixel-based binary model. Each pixel is compared with the preset cover threshold, and the comparison result serves as the switching switch for the subsequent two sub-steps. This division does not discard high-coverage pixels, but rather changes their feature source, thus ensuring that no pixels in the entire area are removed. For example, if the preset cover threshold is 0.4, pixel 1 has a vegetation cover of 0.2 and is classified as a low-coverage pixel, while pixel 2 has a vegetation cover of 0.7 and is classified as a high-coverage pixel.

[0080] S22. For low-coverage pixels, a radiometric correlation feature set is extracted using the main spectral channel, which is based on the alteration spectral characteristics of the shortwave infrared band. Low-coverage pixels have a high proportion of exposed surface, and the radiance received in the shortwave infrared band is dominated by rock and soil. The diagnostic absorption depth can be directly converted into alteration intensity. The output dimension of this channel is the alteration intensity itself, which is the dimensional benchmark for the global features. The output of S23 must ultimately be aligned to this benchmark.

[0081] S23. For high-coverage pixels, a radiometric correlation feature set is extracted using an alternative channel. This alternative channel is based on at least one of the following: synthetic aperture radar (SAR) backscattering features, topographic features, and neighborhood lithology inference features. SAR backscattering features utilize the penetrating power of microwaves through the canopy to extract surface roughness and dielectric properties from the backscattering coefficient and polarization decomposition results. These two properties are related to the degree of bedrock weathering. Topographic features are derived from a digital elevation model (DEM), where the canopy undulates with the terrain without altering the terrain itself; this source is unaffected by cover. Neighborhood lithology inference features are extrapolated from the lithology of surrounding low-coverage pixels using a distance-weighted approach, leveraging the continuity of lithology within the same geological unit to fill gaps in the pixel's data. Any one of these three sources can be used, or they can be used in combination; the number of sources used depends on the types of data actually available for that pixel.

[0082] S24. Select bridging samples from low-coverage pixels where both the main spectral channel and the alternative channel are available, and where the terrain shadow and surface humidity do not exceed preset interference limits. Extract radiometric correlation feature sets from the bridging samples through the alternative channel. Fit the mapping relationship between the output of the alternative channel and the output of the main spectral channel using the bridging samples to complete the bridging calibration, ensuring that the output of the alternative channel and the output of the main spectral channel have consistent dimensions. The bridging sample must simultaneously meet three conditions: the pixel is a low-coverage pixel and therefore the main spectral channel is available; the radar and terrain data for the pixel are complete and therefore the alternative channel is also available; and the terrain shadow intensity and surface humidity index of the pixel do not exceed preset interference limits. The third condition is used to exclude locations where both channels are disturbed by the same external factor. If sampling is performed in a shadow area or a high-humidity area, the fitted mapping relationship will solidify the disturbance into the calibration coefficients. The fitted mapping relationship is then applied to the output of the alternative channel of all high-coverage pixels, so that the feature values ​​of the two types of pixels fall on the same dimension. Therefore, during training in S3, it is not necessary to distinguish which channel the pixel originally belonged to.

[0083] In its implementation, S24 also statistically analyzes the fitting residuals of the bridging calibration on each bridging sample and establishes a channel difference uncertainty model that monotonically increases with vegetation cover based on the vegetation cover of each bridging sample. Although the vegetation cover of the bridging samples is all below the preset cover threshold, there is still a distribution range below the threshold. For samples with cover close to the threshold, the proportion of effective signal in the main spectral channel has decreased, and the inference share undertaken by the substitute channel has increased accordingly, resulting in a larger fitting residual. After statistically analyzing the residuals according to the cover level and fitting them into a monotonically increasing curve, this curve can be extrapolated to the cover range above the threshold, providing an error estimate for channel replacement in high-coverage pixels. This estimate is injected as a component into the prediction uncertainty surface in S4.

[0084] The above describes one implementation method of splitting the dual channels according to coverage and calibrating them using an overlap domain bridging. In other embodiments, S2 may not perform channel splitting, but instead stack the observations from the spectral and radar sources, along with the vegetation coverage itself, as components into the feature vector. The partition mapping model then learns the effect of coverage on the reliability of each component during training. This alternative method does not include a bridging calibration step; dimensional consistency is handled by the weights within the model, achieving the same goal of extracting radiometric correlation features within the vegetation-covered area.

[0085] In some embodiments, S2 further performs feature-level multi-temporal robust synthesis on the radioactive correlation feature set. Feature-level multi-temporal robust synthesis includes: taking robust statistics as feature values ​​for the multi-temporal sequences of each proxy feature, outputting the temporal variation coefficient of each pixel as the feature quality weight, and performing periodic component decomposition on the multi-temporal sequences to distinguish between periodic phenological changes and abrupt surface disturbances; and providing the detection results of abrupt surface disturbances to the surface disturbance detection in S31.

[0086] Synthesis is performed at the feature layer rather than the image layer. That is, a set of surrogate features is first calculated for each image period, and then statistical analysis is performed across periods along the feature dimension. For example, alteration spectral features are calculated once for each of the 12 observations within a year, and the median of the 12 values ​​is taken as the feature value of that pixel. Deviations caused by thin clouds, aerosols, or temporary water accumulation in individual periods are not included in the final value.

[0087] The temporal variation coefficient is obtained by the ratio of the standard deviation to the mean of the feature values ​​over multiple periods for the same pixel. A larger value indicates that the feature value of that pixel is more unstable. This coefficient is not used to filter out pixels, but is output as a feature quality weight along with the feature value. It participates in the weighting in S3 training and S4 inversion, so that the influence of locations with poor observation conditions on the model is reduced accordingly.

[0088] Periodic component decomposition breaks down a multi-period sequence into two parts: a periodic term and a residual term. The periodic term corresponds to seasonal phenological changes that recur and reappear in the same period of the following year; the residual term contains abrupt and non-reversible changes, which correspond to actual disturbances such as excavation, filling, or hardening of the land surface. The distinction between the two is passed downstream as input for S31 to determine whether the surface at the location of the tag point has changed.

[0089] S3. Divide the target area into multiple geological units, and train a partitioning mapping model from radioactivity correlation feature set to natural radioactivity level in each geological unit using historical measured label data.

[0090] In practical implementation, the relationship between proxy characteristics and natural radioactivity levels is controlled by lithology. The nuclide content corresponding to the same alteration intensity value is different in a granite background and in a sedimentary rock background. If a single model is fitted globally, the parameters can only take the compromise values ​​of each lithological zone, and the inference of any lithological zone is skewed by the other lithological zone. Partition training decomposes the global relationship fitting problem into several relationship fitting problems within a unit. The relationship within each unit can be regarded as stationary, and the model parameters do not need to take into account two incompatible material backgrounds.

[0091] In some embodiments, geological units are obtained by fusing regional geological maps with remote sensing lithological classification results. Regional geological maps provide lithological boundaries verified in the field; their classification is reliable, but the map scale is often coarse, and small rock masses and later alteration zones may be incorporated into adjacent units. Remote sensing lithological classification is based on pixel-by-pixel differentiation according to spectral and textural differences in the image, providing rich spatial detail but lacking prior knowledge of stratigraphic attribution. During fusion, the boundaries of the regional geological map serve as the framework. Within this framework, subdivided subregions with significantly different spectral characteristics are determined by the remote sensing lithological classification results. Locations where the two classifications are inconsistent are determined using the regional geological map and marked as boundaries to be verified. Taking the target area as an example, geological unit 1 is a biotite granite unit, geological unit 2 is a Quaternary overburden unit, and geological unit 3 is a metamorphic sandstone unit.

[0092] Historical airborne gamma-ray spectrum data in historical measured tag data do not correspond to satellite pixels in terms of spatial scale and acquisition time. Using them directly as tags would introduce both scale misalignment and outdated time phases into the model. To address this issue, refer to... Figure 3 In some embodiments, historical measured tag data includes historical airborne gamma spectrum data, and S3 includes sub-steps S31 to S34.

[0093] S31. Perform label timeliness verification on historical measured label data. This verification is based on multi-temporal, multi-source satellite remote sensing data, performing surface disturbance detection on each label location. Labels with detected surface disturbances are downweighted or removed to obtain verified labels. The radioactive geological background is quasi-static on a geological timescale. The reason for historical label failure is not due to changes in the radioactive field itself, but rather to alterations to the surface of the label location by subsequent human activities. Surface disturbance detection utilizes all temporal observations from the year the label location was acquired to the present, receiving the abrupt surface disturbance discrimination results output by the S2 periodic component decomposition. Locations where a step change is detected during this period and does not recover are considered to have changed surface conditions. For labels determined to have changed surface conditions, the medium described by their values ​​is no longer the same as the medium currently observed by the satellite. Therefore, their weight in training is reduced according to the disturbance intensity, and those with disturbance intensity exceeding a preset limit are directly removed. For example, tag point 2 is located on the edge of a quarry. After the year of collection, the albedo at this location showed a step increase and did not recover. It was determined that the surface had changed and the intensity of the disturbance exceeded the preset limit. Therefore, the tag was removed. The remaining points retained their original weights and were used together to form the tags after verification.

[0094] S32. Based on the detection field-of-view (PDG) response function of each measurement point in the historical airborne gamma spectrum data, the radioactive correlation feature set is weighted and aggregated to the PEG scale to obtain patch-scale features. The effective width of the PEG response function is determined point-by-point according to the flight altitude record and the terrain undulation of the location. A single reading in airborne measurement is not the value of a specific point, but rather the average of a ground area below that point weighted by contribution. The contribution decreases as the horizontal distance from the nadir point of the measurement point increases, and the attenuation pattern is the PEG response function. This sub-step does not break down the labels to the pixel level, but aggregates the pixel-level surrogate features upward according to the same PEG response function, so that the features and labels are paired at the same scale, i.e., the scale misalignment is resolved by the forward modeling direction rather than the inverse modeling direction. The effective width is determined according to the radar altitude record of each measurement point. Within the same sortie, the flight altitude varies with terrain undulation, and the effective width in the elevated section of mountainous areas is greater than that in the plains section. The terrain undulation itself is further corrected. In locations with large undulations, the actual distance from the detector to the ground is unevenly distributed within the PEG, and the effective width is correspondingly widened. Instead of taking the average of all flights, features are determined point-by-point, ensuring that each label's paired features correspond to its actual observed geometry. The aggregated features are called patch-scale features, and they correspond one-to-one with the verified labels.

[0095] S33. Within each geological unit, sampling is performed on the verified labels according to the characteristic distribution quota to correct the survey corridor bias of historical airborne gamma spectral data, resulting in a training sample set. Airborne surveys fly along the survey line, and the labels are spatially concentrated in the corridors on both sides of the survey line. The terrain features, vegetation conditions, and alteration intensity ranges traversed by the corridors do not represent the overall distribution of the geological unit. Directly sampling all of them would leave the model without constraints on feature regions not covered by the corridors. This sub-step first statistically analyzes the value distribution of all pixels within the geological unit on each proxy feature, divides it into several feature regions according to this distribution, and determines the sample quota for each region. Then, sampling is performed on the verified labels falling into each region according to the quota. Labels in excessively dense regions within the corridor are downsampled, while labels in sparse regions are either fully sampled or have their sampling weight increased. The distribution of the sampled samples in the feature space is consistent with the overall feature distribution of the geological unit, thus correcting the bias caused by the survey line corridor in the feature space. This sub-step is performed on each geological unit separately, and the resulting samples are merged into a training sample set.

[0096] S34. A partition mapping model is trained at the patch scale in the detection field of view based on the training sample set and patch scale features. Training is carried out at the patch scale, the input is the patch scale features, and the target quantity is the verified label value of the corresponding measurement point. The model can be a regressor that can output ensemble member predictions, such as gradient boosting regression tree, random forest, or Gaussian process regression.

[0097] In S4, the mapping models for each region are applied using the original pixel resolution of multi-source satellite remote sensing data. Training is completed at the patch scale, while application reverts to the original pixel resolution. This is because the model fits a functional relationship between surrogate feature values ​​and natural radioactivity levels. This relationship itself does not change with the observation scale; only the spatial averaging of the input features changes. Patch-scale pairing avoids the spurious details introduced by forcibly assigning large-scale labels to small-scale pixels, while pixel-scale application preserves the spatial resolution of the satellite observations themselves.

[0098] In some embodiments, the historical measured tag data also includes ground measured point data and geochemical unmixed tags. The geochemical unmixed tags are obtained by unmixing and allocating the elemental content data of stream sediments through the catchment area. In S34, the training sample set is also weighted according to the noise variance of each tag. The noise variance corresponding to the historical airborne gamma spectrum data is calculated by the Poisson statistical characteristics of gamma count, the flight altitude attenuation factor and the action width of the detection field of view spot response function. The noise variance corresponding to the ground measured point data is determined by the accuracy index of the measuring instrument. The noise variance corresponding to the geochemical unmixed tags is determined by the residual of the unmixing allocation in the catchment area. The noise variance of each tag is unified to the same content dimension before participating in the weighted training.

[0099] The reliability of the three types of labels differs significantly. If they are trained with equal weight, the class with a larger number of labels but lower accuracy will dominate parameter fitting, while the class with a smaller number of labels but higher accuracy will be overwhelmed. Weighted training uses the reciprocal of the noise variance of each label as its weight in the loss function, making reliability an explicit condition for participating in the fitting.

[0100] The noise variance of historical airborne gamma-ray spectrum data is estimated from the physical processes involved in its measurement. The estimation is based on a set of nominal observation conditions, namely, a pre-defined nominal flight altitude and a nominal field-of-view (FOR) width. Deviations of the actual observation conditions at each measurement point from this benchmark are accounted for by a correction factor. The noise variance is determined by the following formula: in, The noise variance of the measurement point label is denoted by the subscript 'a', which indicates that the label comes from historical airborne gamma spectrum data and is used to distinguish it from the noise variance of ground measurement point data and geochemical unmixed labels. The proportionality coefficient for converting count fluctuations to the square of the content dimension is determined by the content scale relationship of the instrument and the nominal observation conditions. The dimension is the product of the square of the content and the length. denoted as the linear attenuation coefficient of air for photons in the corresponding energy window. This represents the difference between the actual flight altitude and the nominal flight altitude at the measuring point; the product of the two is dimensionless. The atmospheric attenuation correction factor characterizes the effective count attenuation factor caused by the absorption of gamma photons by the air along the way from the ground to the detector. Its exponential form comes from the physical fact that the intensity of photons decreases exponentially along the way in a uniform absorbing medium. Here, the height difference is used as the increment of the attenuation path length, corresponding to the case of the detector observing vertically downward. When observing oblique path, the actual path is longer than this increment and must be calculated separately according to the observed zenith angle. The equivalent count received by the measuring point within the integration time is the original count after being reduced to the nominal observation conditions; The effective width of the detection field-of-view spot response function at this measuring point.

[0101] This formula gives three monotonic relations, each corresponding to an independent physical mechanism, and they are not counted repeatedly. Gamma counts follow Poisson statistics; the relative fluctuations of the counts increase as the count decreases. Therefore, the equivalent count after reduction... The smaller the value, the larger the noise variance; this term accounts for statistical fluctuations. When the actual flight altitude is higher than the nominal flight altitude, the path of photons being absorbed by the air becomes longer, and the effective count reaching the detector decreases exponentially, leading to an atmospheric attenuation correction factor. As it increases, therefore The larger the value, the greater the noise variance; this term accounts for the deviation of the actual observed height from the reference. The deviation has been factored into the nominal height and is not included. The pattern repeats itself. A wider effective width results in a larger average ground area for a single reading; spatial averaging itself suppresses random fluctuations, therefore... The larger the value, the smaller the noise variance; this term reflects the spatial averaging of the observation geometry. These three relationships correspond to statistical fluctuations, atmospheric attenuation, and spatial averaging mechanisms, respectively, and together determine the reliability of the label.

[0102] The noise variance of ground-based measured data is determined by the accuracy specifications of the measuring instrument. The instrument's nominal accuracy is given in the form of relative or absolute error, which can be directly converted into variance. This source is preferred over the statistical variance of repeated observations at the same measuring point because historical data typically only records single readings and does not retain records of repeated observations. The instrument's nominal accuracy, however, can always be traced back to the instrument model and is still available when data conditions are limited. If a batch of data does indeed retain records of repeated observations, then the statistical variance of the repeated observations is used instead of the nominal accuracy; whichever is more readily available and closer to the actual measurement conditions is chosen.

[0103] The noise variance of geochemical unmixing tags is determined by the residuals of the unmixing allocation in the catchment area. Unmixing allocation distributes the reading of a sampling point to multiple upstream pixels. The difference between the allocation result and the original reading after it is re-converged in the catchment area is the residual. The larger the residual, the greater the deviation between the allocation and the actual mixing process.

[0104] The variances of the three types of noise are not naturally consistent in terms of dimensions. Airborne data are given as count rates or equivalent content, ground points as dose rates or content, and geochemical tags as elemental content. Comparison is only possible after unifying them to the same content dimension. Dimensional conversion is completed based on the respective instrument calibration relationships and the conversion coefficient from nuclide content to dose rate. Only the converted variances are used in weighted training.

[0105] In some embodiments, the unmixing assignment of catchment areas includes: extracting the catchment areas of each stream sediment sampling point based on a digital elevation model; determining the assignment weights based on the area, lithology, and erosion intensity of each pixel within the catchment area; and using regularized inversion to assign the elemental content data of stream sediments to each pixel within the catchment area, thus obtaining geochemical unmixing labels. Since there are numerous pixels within the catchment area but only one observation, the assignment problem is undetermined. Therefore, a regularization term is used to constrain the inversion by using a gradual change in values ​​between adjacent pixels, making the solution reasonable in terms of material origin. This method restores a mixed reading at the catchment outlet to the shared values ​​of upstream pixels, transforming geochemical results that could only provide qualitative indications into weakly distributed area labels.

[0106] Geological unit 2 is a Quaternary overlying unit, and the number of historically measured labels falling within it is less than required for modeling. In some embodiments, in S3, for geological units where the number of historically measured label data is lower than a preset sample size threshold, a regional mapping model of another geological unit adjacent to the geological unit is used to perform a global inversion, and the value of the geological unit in the prediction uncertainty surface is increased. Adjacent geological attributes are determined based on diagenetic background, major rock-forming mineral assemblage, and provenance relationship, rather than solely on spatial proximity. Borrowing the model prevents the unit from being left blank on the result surface due to a lack of labels, and simultaneously increasing the uncertainty value makes this borrowing fact explicitly marked. Therefore, the high-value portion of the unit enters the verification area rather than the target area during the S5 diversion.

[0107] In other embodiments, each zonal mapping model adopts a hierarchical structure, sharing the same set of prior model parameters. Geological units with fewer historical measured label data are trained with zonal mapping models constrained by the prior model parameters. The model parameters of all units across the entire domain are considered to come from the same upper-level distribution. Geological units 1 and 3 with sufficient labels are mainly determined by their own samples, while geological unit 2 with few labels is more constrained by the upper-level distribution, and its parameters are pulled towards common values ​​across the entire domain. The difference between this approach and the aforementioned borrowing of neighboring unit models is that the small number of samples in this unit still participate in the fitting and affect the parameters, rather than being discarded. Both approaches achieve the goal of obtaining usable zonal mapping models within geological units with insufficient labels.

[0108] The zoning mapping models for each geological unit are independent of each other, and different models provide inferred values ​​on both sides of the unit boundary. Therefore, value jumps that are inconsistent with actual material changes may occur on the boundary line. To address this situation, in some embodiments, between S3 and S4, a transition zone is further included on both sides of the contact zone between two adjacent geological units. In S4, for each pixel located within the transition zone, the prediction results of the respective zoning mapping models of the two adjacent geological units are weighted and fused according to the pixel's membership degree to the two adjacent geological units. The width of the transition zone is set differently according to the geological type of the contact zone; the width of the transition zone corresponding to fault contact is smaller than the width of the transition zone corresponding to intrusive gradual contact.

[0109] The transition zone is set after the geological unit division is completed and before the global inversion begins, that is, at the point where the zoning results are determined but the model has not yet been applied pixel by pixel. This arrangement allows the spatial extent of the transition zone to be directly derived from the determined unit boundaries, while also enabling the fusion process to be completed on-site during the inversion process, without having to generate a result surface with jumps first and then perform post-processing.

[0110] Membership is determined by the vertical distance from the pixel to the center line of the contact zone. Pixels located on the center line receive half the membership for each of the units on both sides. The farther away from the center line, the closer the membership for the unit on the same side is to one, and the closer the membership for the other side is to zero. Pixels outside the transition zone have a membership of zero or one, meaning the result is given by the model of that unit. Fusion is a weighted sum of the outputs of the two models based on membership. The output surface continuously changes from one side to the other when crossing the contact zone.

[0111] The width of the transition zone is set according to the geological type of the contact zone, based on the fact that the spatial scale of material change differs among contact zones of different origins. Fault contact is formed by tectonic displacement, where the lithologies on both sides directly meet at the fault plane, and material change occurs within a few meters to tens of meters. The step change in natural radioactivity level at this point is a real geological phenomenon. If an excessively wide transition zone is given, fusion will erase this real step change as a gentle slope. Intrusive gradual contact is formed by magma intruding into the surrounding rock and contacting and replacing it. There is a contaminated zone within the contact zone where the composition gradually transitions, and its width can reach tens to hundreds of meters. The actual material change at this point is gradual, so the width of the transition zone is correspondingly larger. Taking the target area as an example, the boundary between geological unit 1 and geological unit 3 is a fault contact in the northern section, so the width of the transition zone is set to a smaller value; the boundary in the southern section is an intrusive gradual contact formed by the intrusion of the rock mass into the surrounding rock, so the width of the transition zone is set to a larger value. The differentiated setting ensures that fusion only erases the pseudo-jumps caused by artificial division of the zones, while retaining the real steps caused by the geological processes themselves.

[0112] S4. Perform a global inversion on the target region using the mapping models of each region to generate a radioactivity level prediction surface and a prediction uncertainty surface.

[0113] In practice, the global inversion is performed pixel by pixel. First, the corresponding zoning mapping model is selected according to the geological unit to which the pixel belongs. Pixels located in the transition zone use the models on both sides according to the aforementioned membership degree. Then, the feature vector obtained by the pixel in S2 is input into the model, and the model output is taken as the value of the pixel on the radiometric level prediction surface. After all pixels have been processed, a continuous measurement surface with the same grid and resolution as the satellite image is obtained.

[0114] The prediction uncertainty surface and the radioactivity level prediction surface are generated synchronously in the same inversion. The value of each pixel is synthesized from several uncertainty components. The synthesis method is to add the variances of each component and take the square root. Each component corresponds to the model's own inference confidence, the noise of the label itself, and the bias introduced in the feature extraction process. The three are statistically uncorrelated, so they are treated as additive variances.

[0115] The three sources mentioned above correspond to the dichotomy made later based on whether convergence can be achieved through supplementary observations. The model's own inference is based on the density of training samples within the feature interval, which can be tightened by adding measured labels and is classified as a cognitive uncertainty component. The noise of the labels themselves is propagated from measurement errors through the model and does not disappear with supplementary measurements, thus it is classified as a random uncertainty component. Among the biases introduced in the feature extraction stage, the channel difference uncertainty component introduced by channel replacement is caused by neither insufficient training samples nor label noise propagation, and therefore does not participate in the comparison between the cognitive uncertainty component and the random uncertainty component. Its effect is reflected in the total amount of the prediction uncertainty surface, that is, after this component increases the total amount, the corresponding pixels are more likely to fall into the verification area in the splitting stage. In the embodiment without channel replacement, the prediction uncertainty surface is composed of both the cognitive uncertainty component and the random uncertainty component.

[0116] In some embodiments, the quantification of model inference mastery is taken as the prediction dispersion among the members of the ensemble model. When the partition mapping model is trained in an ensemble manner, it contains multiple base learners, each providing its own inference value for the same pixel. If the feature value of the pixel falls within a densely covered region of the training samples, the inferences of each member tend to be consistent, with low dispersion; if it falls within a sparse or uncovered region of the training samples, each member extrapolates freely in the unconstrained direction, and the dispersion increases accordingly. Using the prediction variance among members as the value of this component allows the model's judgment of its own mastery to be given by the model structure itself, without relying on external priors.

[0117] The above describes one implementation method for quantifying model inference. In other embodiments, quantile regression can also be used, whereby the upper and lower quantiles of the target quantity are fitted simultaneously during training, and the difference between the two quantiles is used as the interval width of the pixel and converted into an uncertainty component; or, in the feature space of the training samples, the nearest neighbor samples are retrieved for each pixel, and the local variance of the fitting residuals of the nearest neighbor samples is used as the value of this component. Quantile regression directly gives the width of the conditional distribution, while the nearest neighbor residual rule is based on sample density and local fitting quality. Both can be summarized by the same higher-level step, thus achieving the goal of outputting the degree of uncertainty of the inferred value at the same time as outputting the inferred value.

[0118] The component of label noise propagated through the model is obtained by passing the noise variance of each label in the S34 weighted training to the output along the model's response relationship. For labels participating in fitting the feature interval of a given pixel, the larger their noise variance, the greater the jitter in the model output within that interval. This component is characterized by not decreasing with additional observations because its cause is not insufficient sample size but rather measurement errors inherent in the samples themselves.

[0119] In some embodiments, in S4, the channel difference uncertainty component is obtained by taking values ​​based on the vegetation cover of each pixel using the channel difference uncertainty model, and the channel difference uncertainty component is injected into the prediction uncertainty surface. The monotonically increasing relationship between the residual obtained from the bridging sample fitting and the coverage in S24 is used as the lookup curve in this step. The corresponding residual estimate is obtained for each pixel based on its vegetation cover, and this estimate is the channel difference uncertainty component of that pixel. The features of low-coverage pixels come directly from the main spectral channel without channel replacement, and their components are taken at the coverage of that pixel according to the channel difference uncertainty model. Due to the low coverage, the value is very small. The fitting residual statistically obtained from the bridging sample reflects the uncertainty of the calibration relationship itself, which is separate from whether the pixel has actually undergone channel replacement. The features of high-coverage pixels come from the substitute channel and are converted by bridging calibration. The higher the coverage, the heavier the inference share undertaken by the substitute channel, and the larger the value of this component. After injection, although the inferred values ​​for densely vegetated areas are still given, their unreliability is simultaneously displayed on the results surface. Therefore, the high values ​​at these locations will not be directly identified as the target area during S5 diversion.

[0120] In some embodiments, S4 further performs feature space extrapolation detection on each pixel, classifying pixels whose values ​​of the radioactive correlation feature set fall outside the training feature distribution of the zoning mapping model into the verification area. Extrapolation detection compares the pixel's feature vector with the feature distribution of the training samples for the geological unit, checking dimension by dimension whether the value falls within the training value range, and checking whether the nearest distance between the vector and the training sample set in the feature space exceeds a preset range. Falling outside the distribution means that the model is extrapolating rather than interpolating at this point, and its output lacks sample support; even if the integrated members happen to be consistent, it does not constitute a reliable basis. These pixels are directly classified into the verification area without threshold comparison, allowing locations outside the model's cognitive boundary to be explicitly transferred to ground work, rather than being mixed into a defined target area with a seemingly definite value.

[0121] S5. The area where the value of the radioactivity level prediction surface exceeds the preset risk threshold is identified as a high-value area. The part of the high-value area where the value of the prediction uncertainty surface is lower than the preset confidence threshold is identified as the target area. The remaining part of the high-value area is identified as the area to be verified.

[0122] In the specific implementation, refer to Figure 4 The delineation process consists of two comparisons. The first comparison is performed on the radioactivity level prediction surface, comparing each pixel with a preset risk threshold. Pixels with values ​​exceeding this threshold are identified as high-value areas, and this comparison determines the locations of concern regarding natural radioactivity levels. The second comparison is performed on the prediction uncertainty surface, only for pixels already falling into the high-value area. Pixels with values ​​below a preset confidence threshold are delineated as the confirmed target area, while those above the threshold are delineated as areas to be verified. This comparison determines locations within these areas with sufficient inference confidence. Therefore, the two surfaces are not primary and secondary, but rather participate in the segmentation in parallel. Any pixel must pass both comparisons simultaneously to enter the confirmed target area. The confirmed target area pixels and areas to be verified pixels obtained from the pixel-by-pixel comparison are grouped into several connected regions based on spatial connectivity. Each connected region is a patch, and these patches constitute the basic operational unit for subsequent source type identification and priority ranking based on area objects.

[0123] The preset risk threshold is set according to the purpose of the survey and the natural background level of the region. When used for radon risk screening in residential areas, the upper quantile of the background statistics under the local geological background can be taken, for example, the 95th percentile of the predicted values ​​of all pixels in the county, so that the threshold reflects the degree of elevation relative to the local normal. When used for the investigation of building stone source areas, the corresponding original rock content level can be calculated by back-calculating the limit standard of radionuclides in building materials. The preset confidence threshold is set according to the acceptable cost of false alarms. The lower the threshold, the more stringent the conditions for entering the target area, the larger the area moved out to the area to be investigated, and the workload of detailed ground investigation increases accordingly, while the false alarm rate decreases. The two thresholds control risks in different directions. The preset risk threshold mainly controls the risk of underreporting, while the preset confidence threshold mainly controls the risk of false alarms. Neither threshold can replace the role of the other.

[0124] For example, if the predicted value of a pixel exceeds a preset risk threshold, but its prediction uncertainty surface value is lower than a preset confidence threshold, that pixel enters the designated target area. Another pixel, whose predicted value also exceeds the preset risk threshold, is located in a densely vegetated area, resulting in a larger channel difference uncertainty component. The combined uncertainty value is not lower than the preset confidence threshold, and this pixel enters the area to be verified. The values ​​of two pixels on the radioactivity level prediction surface may be similar; their different paths are determined by the prediction uncertainty surface.

[0125] In some embodiments, S5 further decomposes the prediction uncertainty surface into a cognitive uncertainty component and a random uncertainty component. The cognitive uncertainty component represents the uncertainty that can be converged through supplementary observations due to label scarcity, while the random uncertainty component represents the uncertainty that cannot be converged through supplementary observations due to label noise propagation. Supplementary observation suggestions are made for the portion of the output in the verification area where the cognitive uncertainty component is greater than the random uncertainty component. The supplementary observation suggestions include suggested observation locations, suggested observation types, and observation priorities. Quality-restricted labels are made for the portion of the output data in the verification area where the random uncertainty component is greater than the cognitive uncertainty component.

[0126] The decomposition is based on the different responses of the two types of uncertainty to supplementary observations. Cognitive uncertainty stems from insufficient training samples within the characteristic interval, resulting in a lack of constraints on the model. Once new measured labels are added to this interval, the model's inference in that interval is pinned, and this component decreases accordingly. Random uncertainty, on the other hand, originates from measurement errors carried by the labels themselves, propagated to the output through the model. Even if one hundred more points are measured at the same location, as long as the instrument accuracy and observation geometry remain unchanged, this component will not decrease. Therefore, for the same area to be verified, the degree of improvement possible differs depending on where the two types of uncertainties are dominant.

[0127] For regions where cognitive uncertainty is dominant, supplementary observation suggestions are output. Suggested observation locations are determined by the pixels within that region whose feature values ​​deviate most from the densest training sample area, as these locations provide the greatest constraint increment from supplementary measurements. Suggested observation types are determined based on the surface conditions and required dimensions at that location; for exposed bedrock, ground gamma-ray spectroscopy measurements can be suggested to obtain nuclide content, while for areas with thicker overburden, combined sampling analysis can be recommended. Observation priority is determined by a comprehensive ranking based on the predicted value of the pixel, the magnitude of the cognitive uncertainty component, and the area of ​​the associated map patch, ensuring that limited detailed survey slots are prioritized for locations with the greatest impact.

[0128] For portions where random uncertainty dominates, the output data quality is marked as limited. This marking does not indicate any remeasurement action, but rather states that the unreliability of the inference at this point stems from the measurement conditions of the existing data itself, and convergence cannot be achieved through remeasurement without changing the data source. Based on this, the operator can temporarily place this portion in a pending verification state or arrange for alternative observation methods.

[0129] The above describes one implementation method for delineation based on pixel-by-pixel threshold comparison. In other embodiments, in S5, delineation is performed on a per-patch basis, with a minimum area constraint applied to each patch. The consistency between the principal axis of each patch's shape and the orientation of its local structure is statistically analyzed, and patches with consistency below a preset consistency limit are given lower target area priority. During morphological operations, the orientation of the structural elements is adaptively determined according to the local structural orientation. Nuclides are enriched along fractures and joints under structural control, and the morphology of patches in the real target area often extends along the structural orientation, while isolated small patches generated by pixel-by-pixel threshold segmentation often originate from noise. The minimum area constraint filters out fragmented patches with excessively small areas, while the principal axis consistency criterion downgrades patches whose morphology does not conform to the structural framework. The orientation of the structural elements ensures that the connectivity of morphological operations along the orientation direction is stronger than that in the vertical direction, avoiding the truncation of a narrow, long target area extending along a fracture during the closing operation. This alternative method also achieves the goal of generating a list of workable areas from the radioactivity level prediction surface and the prediction uncertainty surface.

[0130] Following S5, a source type discrimination operation is performed on each patch within the identified target area. This involves calculating the extrapolated proportion of pixels within each patch that fall outside the training feature distribution, detecting the occurrence time and range evolution information of each patch during multi-temporal observations, and classifying the patches as either anthropogenic or geological sources based on two criteria. The patches are then categorized into anthropogenic radioactive source target areas and geological target areas, respectively, and corresponding recommendations are output. The following provides a detailed explanation of these operations.

[0131] The extrapolated proportion is obtained by statistically analyzing the percentage of pixels whose values ​​of the radioactive correlation feature set within each patch fall outside the training feature distribution of the partitioning mapping model.

[0132] Detect the occurrence time and range evolution information of each patch in multi-temporal, multi-source satellite remote sensing data.

[0133] Image patches whose extrapolation ratio exceeds the preset extrapolation limit and whose occurrence time is detected are identified as human sources, included in the human radioactive source target area, and source term investigation suggestions are output. The source term investigation suggestions include information on the occurrence time and range evolution.

[0134] Map patches whose extrapolation ratio does not exceed the preset extrapolation limit and whose occurrence time is not detected are identified as geological sources, included in the geological target area, and ground detailed investigation suggestions are output.

[0135] The locations of high-level natural radioactivity within the demarcated target area can be categorized into two types based on their origin. One type is of primary geological origin, meaning the bedrock itself is rich in radionuclides; the other is of anthropogenic origin, meaning quarry waste, mine tailings, or backfill soil transported to this location and forming a dump. The subsequent procedures for these two types of target areas differ. Geologically originating target areas require detailed ground surveys to determine the distribution range and concentration levels, while anthropogenically originating target areas require source tracing investigations to determine the material's origin and the responsible party. Therefore, they must be distinguished before handover.

[0136] The two criteria approach the issue from the perspectives of model state and time, respectively. The extrapolation ratio reflects whether the feature values ​​of a map patch fall within the range seen by the model. Materials from man-made stockpiles originate elsewhere, and their mineral assemblage, grain size, and water content differ from the original surface of the local geological unit. Therefore, the surrogate feature values ​​often fall outside the feature distribution of the training samples for that unit, resulting in a high extrapolation ratio. In contrast, the surface materials of target areas with primary geological origins belong to that geological unit, and their feature values ​​fall within the training distribution, resulting in a low extrapolation ratio. This criterion uses the model's own training feature distribution as the cognitive boundary, and its judgment is based on the model state itself, without relying on a prior sample library or morphological template of the stockpile site.

[0137] The occurrence time is obtained through backtracking from multi-temporal observations. This involves examining the surface within the map patch's area period by period in a time series to determine if a step-like, irreversible change occurred. The occurrence time is determined by pinpointing the period from non-existence to existence. Primary geological bodies exist continuously on the remote sensing-observable timescale, so their occurrence time cannot be detected. However, man-made dumping sites have a specific initiation date. This criterion assumes that the starting phase of the multi-temporal observation sequence is earlier than the formation date of the map patch. For historical dumping sites whose formation date is earlier than the starting phase of the available sequence, their occurrence time also cannot be detected. The extent evolution information records the expansion, contraction, or displacement of the map patch's boundary since its appearance. Its morphology can reflect whether the dumping site is continuously being filled, has ceased filling, or is being gradually cleared.

[0138] Both criteria must be met simultaneously for a determination to be made. If the two criteria do not correspond, the source type of the patch will not be determined; it will remain within the designated target area and be transferred according to the original procedures. Its source type will be determined after on-site verification of the material morphology and occurrence. For example, patch 1 is located in a granite outcrop area of ​​geological unit 1. Its extrapolation ratio is low, and the surface showed no abrupt changes in all periods during multi-temporal retrospective analysis, with no detected occurrence time. Therefore, it is determined to be a geological source, included in the geological target area, and a detailed surface investigation recommendation is output. For example, patch 2 is located at the mouth of a gully downstream of a quarry. Its extrapolation ratio exceeds the preset extrapolation limit. During multi-temporal retrospective analysis, the area changed from bare soil to gravel accumulation in a certain year and expanded year by year thereafter. Information on the time and extent of its occurrence was detected, leading to its determination to be an anthropogenic source. It is included in the anthropogenic radioactive source target area, and the output source term investigation recommendation includes records of the occurrence time and the extent of the annual expansion, allowing investigators to trace the material source.

[0139] S6. The detailed survey data obtained from the ground survey of the area to be verified is added as a new label and incorporated into the historical measured label data. The zoning mapping model of the geological unit to which the area to be verified belongs is incrementally updated. Based on the zoning mapping model after the incremental update, S4 and S5 are re-executed to determine the target area.

[0140] In practice, the reason for the emergence of areas awaiting verification is that the model's inferences at these locations are insufficient to be trusted. The only way to eliminate this situation is to obtain traceable measured values ​​at these locations. Detailed ground surveys are an inherent step in the environmental investigation process chain, following the general survey. Their results were originally only used to answer the question of the natural radioactivity level at these locations. This step gives the same batch of results a second role: as new labels fed back into the model. This feedback does not incur additional field costs because the detailed survey work must be carried out in the areas awaiting verification regardless.

[0141] After new labels are incorporated into historical measured label data, they form the updated label set together with existing labels. During incorporation, a noise variance is assigned to the new labels according to the weighting rule of S34. The detailed measured data are directly obtained from ground instruments, and their noise variance is determined according to the accuracy specifications of the measuring instruments. Its magnitude is lower than that of historical airborne gamma-ray spectroscopy data and geochemical unmixed labels, thus receiving a higher weight during training. Incremental updates are only performed on the geological units to which the area to be verified belongs; the zoning mapping models of other geological units are unaffected. This avoids the computational overhead of full-domain retraining and also prevents new samples from one unit from perturbing the parameters of unrelated units.

[0142] After the update is complete, S4 and S5 are re-executed. During the re-inversion, the sample density of the feature interval where the new label is located increases, and the cognitive uncertainty component of that interval decreases. Pixels that were originally classified as pending verification areas because their uncertainty was not lower than the preset confidence threshold will have the following fates: Pixels that are confirmed to be high-value after detailed investigation will be transferred to the target area after their uncertainty decreases; Pixels that are confirmed to be below the preset risk threshold after detailed investigation will have their inferred value in the feature interval lowered, and these pixels will be removed from the high-value area; Pixels whose unreliability mainly comes from the random uncertainty component will remain in the pending verification area; Pixels not covered by detailed investigation in this round will be re-evaluated according to the updated partition mapping model, and those whose uncertainty is still not lower than the preset confidence threshold will continue to remain in the pending verification area for processing in subsequent rounds. Each cycle thus reduces the portion of the pending verification area that can be converged through supplementary testing. The first two fates will ensure that the corresponding pixels are no longer in the pending verification state, and the general survey and detailed survey processes will form a closed loop that connects them.

[0143] In some embodiments, in S6, the detailed measured data is the data obtained by performing a detailed ground survey on the portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, according to the supplementary observation recommendations. After each round of re-execution of S4 and S5, the preset convergence condition is evaluated: if there is no portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, it is determined that the preset convergence condition is met; if there is a portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, it is determined that the preset convergence condition is met if the decrease in the area of ​​the area to be verified in this round relative to the area of ​​the area to be verified in the previous round and the decrease in the total amount of the cognitive uncertainty component in this round relative to the total amount of the cognitive uncertainty component in the previous round are both lower than the preset convergence limit, and the first round is determined to be not to meet the preset convergence condition; when the preset convergence condition is met, the next round of incremental update is terminated, and an undecidable area list is output. The undecidable area list includes the remaining areas to be verified at this time, and the proportion of the random uncertainty component in the predicted uncertainty surface value in the areas to be verified.

[0144] The detailed investigation will no longer cover the entire area to be verified, but will instead focus only on the portion where the cognitive uncertainty component is dominant, as suggested by the supplementary observation recommendations. Even if measured, the portion where the random uncertainty component is dominant will not make the model's inferences in that interval more reliable; allocating detailed investigation slots there is tantamount to wasting resources on a direction that will not converge.

[0145] The convergence determination provides two termination paths. The first is the disappearance of the region dominated by the cognitive uncertainty component. This means that in the current verification area, there is no longer a portion where the cognitive uncertainty component is greater than the random uncertainty component. At this point, the unreliability of all remaining verification areas stems from the measurement conditions of existing data, and further measurements would be pointless. The closed loop has no further investment direction, and convergence is determined. The second is the decrease in magnitude path. This means that while the verification area still has a portion dominated by the cognitive uncertainty component, the decrease in the area of ​​the verification area relative to the previous round and the decrease in the total amount of cognitive uncertainty relative to the previous round are both lower than the preset convergence limit. This indicates that the improvement gained from continued investment is minimal, and convergence is determined. Both decrease magnitudes refer to relative decreases, i.e., the ratio of the difference between the previous and current round values ​​to the previous round value. The preset convergence limit is a dimensionless ratio. These two paths correspond to the two scenarios of no investment opportunities and inefficient investment, respectively. Termination occurs upon the success of either path. Since there is no data from the previous round in the first round, the two decrease rates cannot be calculated. Therefore, the first round is judged to be unsatisfied in terms of the decrease rate path. The disappearance path of the dominant region of cognitive uncertainty component does not depend on the data from the previous round, and its judgment is equally valid and takes precedence in the first round.

[0146] Upon termination, a list of undetermined areas is output. This list is not a discarded area, but rather an explicitly transferred area. It includes the remaining areas to be verified, and the proportion of random uncertainty components in the predicted uncertainty surface values ​​for each area. The closer this proportion is to one, the more likely the unreliability stems entirely from the data itself, allowing the team to determine whether to change the observation method rather than continue with additional measurements. Transferring this list ensures that unanswered questions during the survey phase become a transferable work item, rather than simply presenting a value without credibility on the results map.

[0147] In the extreme case where the dominant component of cognitive uncertainty is empty, the closed loop is executed in a degenerate manner for one round. Supplementary observations suggest that this is an empty set, and the detailed experimental data obtained accordingly is also an empty set. The empty set is merged into the label set without changing the content of the label set. Incremental updates degenerate into identity updates under the condition of zero new labels, and the model parameters remain unchanged. The results of re-executing S4 and S5 are the same as the previous round. Subsequently, the loop terminates at the evaluation point by hitting the empty set criterion. Therefore, S6 has a definite execution result under any input conditions and does not reach a state where it cannot proceed.

[0148] On the other hand, the number of pixels in the target area is limited, and the labels added in each round of detailed investigation accumulate monotonically. The cognitive uncertainty component decreases as the sample density of its feature interval increases. As a result, the part dominated by the cognitive uncertainty component disappears within a limited number of rounds, and the closed loop will not continue indefinitely.

[0149] In some embodiments, after S6, a comprehensive score is performed on each patch within the identified target area based on the predicted radioactivity level, patch area, population exposure weight, and accessibility. The comprehensive score is then used to sort and output a priority list for detailed investigation. The predicted value reflects the risk intensity, the patch area reflects the impact range, the population exposure weight is determined by the resident population density and residential land ratio within the patch and its surrounding area, and accessibility is determined by the distance to the nearest accessible road and the difficulty of traversing the terrain. The first two factors determine the priority of the location, and the latter two determine the ease of implementation. The four factors are weighted and ranked, transforming the target area identification from a list of parallel items into a ordered list of tasks.

[0150] The aforementioned closed loop terminates based on a convergence criterion. In other embodiments, it may also terminate at a preset maximum number of rounds, i.e., stopping once the set number of rounds has been completed; or it may terminate based on a detailed survey budget, i.e., stopping once the cumulative detailed survey workload reaches the budgeted amount. Both methods do not rely on the statistics of uncertainty components, making implementation simpler, and similarly achieving the goal of ending the cycle of census and detailed survey within a limited investment.

[0151] The above describes the specific implementation of each step of this method. The following section further explains the mapping relationship between this method and the executing entity. This method can be executed by a natural radioactive environment survey system, which includes a data acquisition module, a feature extraction module, a zoning modeling module, an inversion output module, a target area delineation module, and a closed-loop update module.

[0152] The data acquisition module is configured to acquire multi-source satellite remote sensing data covering the target area, as well as historical measured label data carrying the measured values ​​of natural radioactivity levels in some areas within the target area, and to perform spatiotemporal benchmark unification and multi-temporal synthesis on the multi-source satellite remote sensing data.

[0153] The feature extraction module is configured to extract a set of radioactive correlation features from multi-source satellite remote sensing data and split the pixels into the main spectral channel and the alternative channel according to vegetation cover, so as to complete the bridging calibration between the two channels with bridging samples.

[0154] The zoning modeling module is configured to divide the target area into multiple geological units, perform label timeliness verification, detection field patch-scale aggregation and feature distribution quota sampling on historical measured label data, and train the zoning mapping model in each geological unit according to the noise variance of each label.

[0155] The inversion output module is configured to perform a global inversion on the target area using the mapping models of each partition, generate a radiometric level prediction surface and a prediction uncertainty surface, and then weight and fuse the prediction results of the two models in the transition zone of adjacent geological units according to their membership degree.

[0156] The target area delineation module is configured to identify areas where the predicted radioactivity level exceeds a preset risk threshold as high-value areas. These high-value areas are then divided into identified target areas and areas to be verified based on the predicted uncertainty surface and the preset confidence threshold. For areas to be verified, supplementary observation suggestions and data quality limitation labels are output according to the dominance of the cognitive uncertainty component and the random uncertainty component. Finally, source type discrimination is performed on each patch included in the identified target area.

[0157] The closed-loop update module is configured to incorporate the detailed survey measured data as new labels into the historical measured label data, perform incremental updates on the zoning mapping model of the geological unit to which the area to be verified belongs, drive the inversion output module and the target area delineation module to be re-executed, output a list of undeterminable areas when the preset convergence conditions are met, and perform comprehensive scoring and priority ranking on each map patch included in the determined target area, outputting a detailed survey priority list.

[0158] The same or similar parts in the various embodiments of this application may be referred to each other; the technical features of each embodiment may be combined to form new embodiments when there is no conflict.

[0159] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is used as an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device / system can be divided into different functional units or modules to complete all or part of the functions described above.

[0160] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing, characterized in that, The steps are as follows: S1. Acquire multi-source satellite remote sensing data covering the target area, as well as historical measured tag data carrying the measured values ​​of natural radioactivity levels in a portion of the target area; S2. Extract a radioactive correlation feature set from the multi-source satellite remote sensing data, wherein the radioactive correlation feature set includes multiple proxy features characterizing the geological background of natural radionuclide enrichment; S3. Divide the target area into multiple geological units, and train a partitioning mapping model from the radioactivity association feature set to the natural radioactivity level in each geological unit using the historical measured label data; S4. Perform a global inversion on the target region using the partition mapping models described above to generate a radioactivity level prediction surface and a prediction uncertainty surface; S5. The region where the value of the radioactivity level prediction surface exceeds the preset risk threshold is identified as a high-value region. The part of the high-value region where the value of the prediction uncertainty surface is lower than the preset confidence threshold is identified as the target area. The remaining part of the high-value region is identified as the area to be verified. S6. The detailed survey data obtained by performing a detailed ground survey on the area to be verified is added as a new label and incorporated into the historical measured label data. The zoning mapping model of the geological unit to which the area to be verified belongs is incrementally updated. Based on the zoning mapping model after the incremental update, S4 and S5 are re-executed, and the determined target area is output.

2. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 1, characterized in that, In step S5, the prediction uncertainty surface is further decomposed into a cognitive uncertainty component and a random uncertainty component. The cognitive uncertainty component represents the uncertainty that can be converged by supplementary observation due to the scarcity of labels, and the random uncertainty component represents the uncertainty that cannot be converged by supplementary observation due to the propagation of label noise. For the portion of the area to be verified where the cognitive uncertainty component is greater than the random uncertainty component, supplementary observation suggestions are output. The supplementary observation suggestions include suggested observation locations, suggested observation types, and observation priorities. The output data in the area to be verified where the random uncertainty component is greater than the cognitive uncertainty component is marked as having quality limitations.

3. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 2, characterized in that, In step S6, the detailed survey measured data is the data obtained by performing a detailed ground survey on the portion of the area to be verified where the cognitive uncertainty component is greater than the accidental uncertainty component, according to the supplementary observation suggestion; After each round of re-execution of S4 and S5, it is evaluated whether the preset convergence condition is met: if there is no part in the area to be checked in this round where the cognitive uncertainty component is greater than the accidental uncertainty component, it is determined that the preset convergence condition is met; When there is a portion in the area to be verified in this round where the cognitive uncertainty component is greater than the random uncertainty component, the decrease in the area of ​​the area to be verified in this round relative to the area of ​​the area to be verified in the previous round and the decrease in the total amount of the cognitive uncertainty component in this round relative to the total amount of the cognitive uncertainty component in the previous round are both lower than the preset convergence limit, which is considered to satisfy the preset convergence condition, and the first round is considered not to satisfy the preset convergence condition; when the preset convergence condition is satisfied, the next round of incremental update is terminated, and an undecidable area list is output. The undecidable area list includes the remaining areas to be verified at this time, and the proportion of the random uncertainty component in the area to be verified in the value of the predicted uncertainty surface.

4. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 1, characterized in that, The historical measured tag data includes historical airborne gamma spectrum data, and step S3 includes the following sub-steps: S31. Perform a label timeliness check on the historical measured label data. The label timeliness check is based on the multi-temporal multi-source satellite remote sensing data to perform surface disturbance detection on each label location, and downweight or remove the labels that detect surface disturbance to obtain the checked labels. S32. The radioactive correlation feature set is weighted and aggregated to the detection field of view (PDV) scale according to the detection field of view (PDV) response function of each measurement point in the historical airborne gamma spectrum data to obtain the patch scale features. The effective width of the PPV response function is determined point by point according to the flight altitude record and the terrain undulation of the location of each measurement point. S33. Sampling is performed on the verified labels according to the characteristic distribution quota within each geological unit to correct the survey corridor bias of the historical airborne gamma spectral data and obtain a training sample set; S34. Based on the training sample set and the patch scale features, train the partition mapping model at the patch scale of the detection field of view; In step S4, the original pixel resolution of the multi-source satellite remote sensing data is used to apply each of the partition mapping models.

5. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 4, characterized in that, The historical measured tag data also includes ground measured point data and geochemical unmixed tags. The geochemical unmixed tags are obtained by unmixing and allocating the elemental content data of stream sediments through the catchment area. In step S34, the training sample set is also subjected to weighted training according to the noise variance of each tag. The noise variance corresponding to the historical airborne gamma spectrum data is calculated by the Poisson statistical characteristics of gamma counts, the flight altitude attenuation factor, and the action width of the detection field of view spot response function. The noise variance corresponding to the ground measured point data is determined by the accuracy index of the measuring instrument. The noise variance corresponding to the geochemical unmixed tags is determined by the residual of the unmixing allocation in the catchment area. The noise variance of each tag is unified to the same content dimension before participating in the weighted training.

6. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 4, characterized in that, In step S2, feature-level multi-temporal robust synthesis is also performed on the radioactive correlation feature set. The feature-level multi-temporal robust synthesis includes: taking robust statistics as feature values ​​for the multi-temporal sequences of each proxy feature, outputting the temporal variation coefficient of each pixel as the feature quality weight, and performing periodic component decomposition on the multi-temporal sequences to distinguish between periodic phenological changes and abrupt surface disturbances; and providing the detection results of the abrupt surface disturbances to the surface disturbance detection in step S31.

7. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 1, characterized in that, S2 includes the following sub-steps: S21. Based on the comparison result between vegetation coverage and a preset coverage threshold, the pixels of the target area are divided into low-coverage pixels and high-coverage pixels; S22. Extract the radioactive correlation feature set from the low-coverage pixels through the main spectral channel, wherein the main spectral channel is based on the alteration spectral features of the short-wave infrared band; S23. Extract the radiometric correlation feature set from the high-coverage pixels through an alternative channel, wherein the alternative channel is based on at least one of synthetic aperture radar backscattering features, topographic structure features, and neighborhood lithological inference features; S24. Select pixels from the low-coverage pixels where both the main spectral channel and the alternative channel are available, and where the terrain shadow and surface humidity do not exceed the preset interference limit, as bridging samples. Extract the radiometric correlation feature set from the bridging samples through the alternative channel. Fit the mapping relationship between the output of the alternative channel and the output of the main spectral channel using the bridging samples to complete the bridging calibration, so that the output of the alternative channel and the output of the main spectral channel have the same dimensions.

8. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 7, characterized in that, In step S24, the fitting residuals of the bridging calibration on each bridging sample are also statistically analyzed, and a channel difference uncertainty model that monotonically increases with vegetation coverage is established for the fitting residuals according to the vegetation coverage of each bridging sample; in step S4, the channel difference uncertainty components are obtained by taking values ​​of the channel difference uncertainty model according to the vegetation coverage of each pixel, and the channel difference uncertainty components are injected into the prediction uncertainty surface.

9. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 1, characterized in that, Following step S5, the following source type determination operation is performed on each patch contained in the determined target region: The extrapolated proportion is obtained by statistically analyzing the proportion of pixels whose values ​​of the radioactive correlation feature set within each of the aforementioned patches fall outside the training feature distribution of the partitioning mapping model. Detect the occurrence time and range evolution information of each patch in the multi-temporal multi-source satellite remote sensing data; The patches whose extrapolation ratio exceeds the preset extrapolation limit and whose occurrence time is detected are identified as human sources, classified into human radioactive source target areas, and source term investigation suggestions are output. The source term investigation suggestions include the occurrence time and the range evolution information. The patches whose extrapolation ratio does not exceed the preset extrapolation limit and whose occurrence time is not detected are identified as geological sources, included in the geological target area, and ground detailed investigation suggestions are output.

10. The method for rapid survey and target delineation of natural radioactive environments based on satellite remote sensing according to claim 1, characterized in that, Between S3 and S4, the method further includes: setting transition zones on both sides of the contact zone between two adjacent geological units; in S4, for each pixel located in the transition zone, the prediction results of the respective partition mapping models of the two adjacent geological units are weighted and fused according to the membership degree of each pixel to the two adjacent geological units; the width of the transition zone is set differently according to the geological type of the contact zone, and the width of the transition zone corresponding to fault contact is smaller than the width of the transition zone corresponding to intrusive gradient contact.

Images