Residual agricultural film three-dimensional distribution inversion method based on satellite radar collaborative calibration

Abstract

The invention discloses a residual agricultural film three-dimensional distribution inversion method based on satellite radar collaborative calibration, and belongs to the technical field of agricultural environment remote sensing monitoring, and the method comprises the steps: extracting surface environment characteristics and underground dielectric characteristics from a multispectral satellite image and a ground penetrating radar signal; inversion of surface agricultural film content hot area distribution is carried out by using a random forest and a Kriging interpolation model, and soil internal dielectric constant distribution is acquired from radar data based on a hybrid inversion method; the method comprises the following steps: constructing an encoder-decoder model based on Transform, and establishing a cross-modal mapping relation between satellite features and radar parameters through a multi-head attention mechanism; and fusing earth surface and underground inversion results to generate a three-dimensional space distribution diagram of the residual agricultural film. The method effectively overcomes the technical limitations that traditional optical remote sensing can only identify earth surface coverage and the on-site sampling space is discontinuous, and provides a reliable technical means for farmland residual film pollution assessment and treatment.

Description

TECHNICAL FIELD

[0001] The application relates to a residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration, and belongs to the technical field of agricultural environment remote sensing monitoring. BACKGROUND

[0002] Residual agricultural films are widely distributed on the soil surface and are broken and migrated to the soil interior under the action of plowing, freezing and thawing, forming a three-dimensional distribution pattern extending from the ground to the underground. Under the joint action of soil medium differences, plowing disturbance and environmental changes, residual agricultural films present significant spatial distribution unevenness and discontinuity at different depths and regions. How to accurately depict the three-dimensional spatial pattern of residual agricultural films and realize high-precision inversion and quantitative characterization has become a key bottleneck for fine identification and risk control of residual agricultural films.

[0003] At present, the distribution identification of residual agricultural films mainly adopts optical remote sensing interpretation and field sampling, and there are two limitations: first, optical remote sensing can only obtain the ground coverage area of residual agricultural films and cannot identify the distribution amount at different depths underground; second, although field sampling can obtain the distribution amount of residual agricultural films at different depths, the sampling points are sparse and dispersed, and it is difficult to form continuous high-precision distribution data. SUMMARY

[0004] The application aims to solve the technical problems in the prior art that optical remote sensing cannot identify the distribution amount at different depths underground and field sampling can obtain depth information but is not spatially continuous, and provides a residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration.

[0005] A residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration comprises the following steps: S1, satellite remote sensing and ground penetrating radar data extraction: feature factors are extracted from multispectral satellite images, the feature factors include terrain features, climate features, spectral principal component variables and plastic indexes, which together constitute a satellite remote sensing environmental variable set; A-scan and B-scan signals of the ground penetrating radar are obtained along the measuring line, the electromagnetic wave phase velocity, target depth and relative dielectric constant are calculated through hyperbolic fitting, the return envelope signal is subjected to anti-gain and geometric compensation, and the slope is calculated through sliding window logarithmic linear fitting to calculate the attenuation coefficient, the depth-sensitive dielectric features are constructed based on the depth weighting principle, and the radar feature subset is formed; S2, satellite remote sensing inversion of ground residual agricultural film hot area: based on the satellite remote sensing environmental variable set, a random forest model and a Kriging interpolation model are used for spatial prediction and mapping of the ground residual agricultural film content to generate a ground agricultural film content prediction grid; S3. Dielectric constant of residual agricultural film in soil inversion by ground penetrating radar: Construct and train a deep learning model to invert ground penetrating radar amplitude data into electromagnetic wave velocity, and then calculate the relative dielectric constant distribution of residual agricultural film in soil as the result of underground dielectric constant inversion. S4. Construction of a joint inversion model for spectral and dielectric constant based on multi-head attention mechanism: Construct a Transformer encoder-decoder model, take satellite remote sensing environmental variables as input, and learn the nonlinear mapping relationship between them and ground-penetrating radar amplitude and two-way time through multi-head attention mechanism to realize cross-modal joint inversion from satellite features to radar parameters; S5. 3D pixel-level visualization output of residual agricultural film: Integrate the predicted raster of surface agricultural film content and the inversion results of underground dielectric constant to generate and output a 3D spatial distribution map of residual agricultural film.

[0006] Specifically, in step S1, the terrain features include elevation, slope, and aspect; the climate features include precipitation, surface temperature, and temperature-vegetation drought index; the spectral principal component variables include absorption peak depth, inter-band ratio, and derivative spectral characteristics; and the plastic indices include advanced plastic greenhouse index and improved plastic greenhouse index.

[0007] The formula for calculating the Advanced Plastics Greenhouse Index (APGI) is as follows:

[0008]

[0009] in, For coastal band surface reflectance, The red band represents the surface reflectance. Near-infrared surface reflectance, The surface reflectance is the second band of shortwave infrared.

[0010] The formula for calculating the improved plastic greenhouse index (MPGI) is as follows:

[0011]

[0012] in, The surface reflectance in the super blue band. The surface reflectance is the first band of shortwave infrared.

[0013] The formula for the Temperature-Vegetation-Drought Index (TVDI) is:

[0014]

[0015] in, The surface temperature of the target pixel; and These represent the minimum and maximum surface temperatures of the wetted edge within the characteristic space of the large triangle, respectively.

[0016] Specifically, step S1 further includes a preprocessing step for the ground-penetrating radar amplitude data: using a graph-based window function to filter the raw radar data to remove direct wave interference, wherein the graph-based window function formula is:

[0017]

[0018] in, Let M be the index of the discrete sample points of the window function, where M = (N-1) / 2, and N is the length of the window function. It is a constant between 0 and 1.

[0019] Specifically, step S2 includes: S21. Read sample data containing geographic coordinates, target values ​​of agricultural film content and multiple environmental variables from the geographic database, and read the corresponding environmental variable raster data. Perform invalid value filtering, spatial reference and resolution unification and data type conversion on the data to form the model input dataset. S22. Divide the sample data into a training set and a validation set. Train the random forest model on the training set through multiple iterations and select the model with the highest determination coefficient on the internal validation set as the optimal model. S23. Simultaneously construct a vectorized optimized Kriging interpolation model, calculate the semi-variogram based on sample points and fit a spherical model, and calculate the prediction weights in batches through matrix operations. S24. The optimal random forest model and the vectorized kriging interpolation model are used to predict the environmental variable raster, respectively, to obtain two types of prediction results; based on the comparison of the accuracy of the two types of prediction results on the validation set, the best raster for predicting the agricultural film content on the ground surface is generated and stored in the geographic database.

[0020] Specifically, step S3 includes:

[0021] S31. Construct a training dataset covering the common electromagnetic wave velocity range in ground penetrating radar applications, select the electromagnetic wave velocity range of 0.054 to 0.212 m / ns, the corresponding relative permittivity range is 2 to 30, and generate a velocity model containing 3 to 10 random layers.

[0022] S32. Convert the velocity model in the spatial depth domain to a time depth domain model, using the following formula:

[0023]

[0024] Where t is the two-way propagation time of the ground-penetrating radar electromagnetic wave from transmission to target reflection and return to the receiver, d is the one-way propagation distance of the electromagnetic wave in the medium from the antenna to the reflecting interface, and v is the propagation speed of the electromagnetic wave in the corresponding medium.

[0025] S33. The two-dimensional time-domain finite difference method is used to perform forward modeling on the converted time-depth domain velocity model to generate the corresponding ground-penetrating radar amplitude data.

[0026] S34. Construct a deep learning-based inversion objective function. Using the trained nonlinear mapping function, convert the observed ground-penetrating radar amplitude data into corresponding radar wave velocity properties, thereby achieving the indirect inversion of the dielectric constant. The formula is as follows:

[0027]

[0028] in, Let c be the phase velocity of the electromagnetic wave in the target medium, and c be the speed of light. , is the relative permittivity of the residual agricultural film.

[0029]

[0030] Where m is the data tag, i.e., the speed of electromagnetic waves. To convert the observation data into a mapping function for radar wave velocity, where d is the measured ground-penetrating radar data, The set of parameters to be optimized.

[0031] Specifically, step S4 includes:

[0032] S41. Collaborative Sample Construction: Multispectral satellite raster data and ground-penetrating radar (GPR) survey line data are mapped from point to pixel. The satellite pixel location of each GPR sampling point is determined using bilinear interpolation. For each GPR sampling point, its corresponding 11-dimensional satellite feature vector is extracted to construct a training sample, where the input features are... The output labels are ;

[0033] S42. Preprocessing, normalization, and embedding representation of satellite variables:

[0034] Univariate normalization was performed on each of the 11 input satellite feature variables. Standardize using the following formula:

[0035]

[0036] in, , These are the mean and standard deviation, respectively, obtained from the training set.

[0037] The 11-dimensional normalized input vector of each pixel is mapped to an embedding vector of dimension d_model through linear projection, as shown in the formula:

[0038]

[0039] In the formula, Projection matrix , For the original input ;

[0040] Introducing spatial location coding and temporal coding, respectively for the pixel center coordinates and observation time Constructing learnable location embeddings:

[0041]

[0042] The final encoder input is:

[0043]

[0044] S43. Transformer Model Construction:

[0045] A Transformer model employing an encoder-decoder architecture is used. The encoder consists of multiple stacked layers, each containing a multi-head self-attention mechanism, a feedforward network, and a residual connection structure, used to encode input features and output satellite conditional features. ;

[0046] The decoder uses two learnable query vectors. As input, updates are performed through a self-attention mechanism:

[0047]

[0048] in, For the learnable linear projection parameters in the attention mechanism, the input representation is... These are mapped to query, key, and value vectors, respectively. The feature dimension of the key vector serves as the normalization scale in the scaled dot product attention. To control the magnitude of the dot product and stabilize the softmax distribution;

[0049] Get the updated query representation , As a query, the encoder output Perform cross-attention computation:

[0050]

[0051] Conditional representation is obtained

[0052] Entering the feedforward network:

[0053]

[0054] Where x represents each row vector of C, denoted as This is the learnable weight matrix for the first-level linear transformation; σ is the learnable weight matrix for the second-level linear transformation; σ is the nonlinear activation function. These are the learnable bias vectors for the first and second layer linear transformations, respectively.

[0055] The final output of the decoder is obtained through residual connections and layer normalization. ;

[0056] Two consecutive predicted values ​​are obtained through independent linear mapping heads:

[0057]

[0058] in , Each is a matrix The two row vectors correspond to two queries respectively. Updated task-specific representation; , These represent the learnable weights of the two regression heads, respectively. , These represent the learnable bias terms of the two regression heads, respectively; , These represent two consecutive predicted values ​​output by two independent linear regression heads;

[0059] S44. Loss Function and Training Strategy:

[0060] The model training uses a weighted mean squared error loss function and introduces an L2 weight regularization term. For the There are n samples, and the sample weights are 1. The real label is The model predicts that The weighted mean square error of the two objectives is defined as follows:

[0061]

[0062] in, Let A be the weighted mean square error of the target variable A. Let T be the weighted mean square error of the target variable T, and N be the total number of samples involved in the loss calculation.

[0063] The overall loss function is:

[0064]

[0065] in, , For loss weighting coefficients, The regularization coefficient is . This represents the set of model parameters.

[0066] Specifically, step S5 includes:

[0067] S51. Generate visualization results including two-dimensional and three-dimensional views; the two-dimensional view supports the overlay display of remote sensing base map, agricultural film distribution map, ground-penetrating radar survey line location and sampling point information, and supports the simultaneous display of the original ground-penetrating radar profile, processed profile and interpretation result profile; the three-dimensional view generates a three-dimensional spatial distribution map of residual agricultural film content in farmland by fusing satellite remote sensing classification results and depth information obtained by ground-penetrating radar inversion, and supports viewing the profiles of main survey lines and connecting survey lines in any area as well as interactive browsing of three-dimensional rendering results; S52. Export the agricultural film distribution map, burial depth map, and risk assessment map in raster or vector format; statistically analyze and output the spatial distribution characteristics of residual agricultural film, including the residual area of ​​different grades, average burial depth, and estimated total mass; automatically generate a standardized test report, which includes project overview, technical methods, result charts, and comprehensive conclusions.

[0068] The beneficial effects of this invention are:

[0069] This invention is based on multi-source information acquired collaboratively by satellite remote sensing and ground radar. It constructs a cross-modal feature association and fusion inversion framework, and performs joint analysis on the surface satellite spectral response and the underground radar dielectric constant. This enables high-precision inversion and quantitative characterization of the three-dimensional spatial distribution of residual agricultural film from the surface to the ground. It effectively overcomes the technical limitations of traditional optical remote sensing, which can only acquire two-dimensional information of the surface and has discontinuous sampling space. This provides reliable quantitative technical support for the assessment and precise treatment of agricultural film pollution. Attached Figure Description

[0070] Figure 1 This is a flowchart of the method of the present invention;

[0071] Figure 2 This is a map showing the distribution of residual agricultural film sampling points in a county-level city (A) in northern China.

[0072] Figure 3 This is a schematic diagram of satellite data preprocessing.

[0073] Figure 4 This is a multi-profile combined detection map from ground-penetrating radar;

[0074] Figure 5 A comparison chart showing the fitting of satellite remote sensing inversion results of agricultural film residue with sampled values;

[0075] Figure 6 The dielectric constant and simulated residual amount of soil film at different depths;

[0076] Figure 7 This is a three-dimensional spatial distribution map of the volume fraction of residual agricultural film. Detailed Implementation

[0077] Example 1: As Figure 1 As shown, a method for inverting the three-dimensional distribution of residual agricultural film based on satellite-radar collaborative calibration includes the following steps: Step 1, extraction of satellite remote sensing-ground penetrating radar data.

[0078] First, feature factors are extracted from multispectral satellite images. Principal component analysis (PCA) is used to eliminate redundant information between bands, compressing the multispectral image information into at least a few representative principal component variations to obtain multispectral principal component variables, including absorption peak depth, inter-band ratio, and derivative spectral characteristic values. Second, plastic indices are calculated based on reflectance information from different bands, constructing advanced plastic greenhouse indices and improved plastic greenhouse indices. Finally, topographic feature parameters of the target area are obtained, including elevation, slope, and aspect, as well as relevant climate feature parameters such as precipitation, surface temperature, and temperature-vegetation drought index. These features collectively constitute the satellite remote sensing environmental variable set {X}, the specific contents of which are shown in Table 1.

[0079] Table 1: Characteristic Factors for Retrieval of Residual Agricultural Film

[0080] Class Feature Factor Topographic Feature Elevation, Slope, Aspect Climatic Feature Precipitation, Land Surface Temperature, Temperature Vegetation Dryness Index Spectral Principal Component Variable Absorption Peak Depth, Band Ratio, Derivative Spectral Feature Plastic Index Advanced Plastic Greenhouse Index, Modified Plastic Greenhouse Index

[0081] The formula for calculating the Advanced Plastics Greenhouse Index (APGI) is as follows:

[0082]

[0083] In the formula, These represent the surface reflectance of the coastal and aerosol, red, near-infrared, and shortwave infrared bands, respectively.

[0084] The formula for calculating the Modified Plastics Greenhouse Index (MPGI):

[0085]

[0086] The formula for the Temperature Vegetation Drought Index (TVDI) is:

[0087]

[0088] In the formula, The surface temperature of the target pixel; and These represent the minimum and maximum surface temperatures of the wetted edge within the "large triangle" characteristic space, respectively.

[0089] Based on ground radar data, A-scan and B-scan signals are acquired along the survey line. Electromagnetic wave phase velocity v and target depth are calculated through hyperbolic fitting. The relative permittivity ε = (c / v)² is used. Inverse gain and geometric compensation are applied to the echo envelope signal, and the slope is obtained through sliding window logarithmic linear fitting. The attenuation coefficient α = -b / v is then calculated. A depth-sensitive dielectric feature Q(z) is constructed based on the depth-weighted principle, forming a radar feature subset related to agricultural film content. .

[0090] Preprocessing is performed on the amplitude data of the ground-penetrating radar (GPR). Direct waves, as the main interference signal in GPR, have amplitudes much larger than those of the target's reflected echo. Therefore, removing direct waves helps in target reflection identification and neural network training in step three. A graphical window function is chosen to filter the GPR data to remove direct wave interference. The formula for the graphical window function is as follows: 1.1.

[0091]

[0092] In the formula, M = (N-1) / 2, where N is the length of the window function. It is a constant between 0 and 1.

[0093] Step 2: Satellite remote sensing inversion of hot zones of residual agricultural film on the ground surface

[0094] First, sample point data containing longitude, latitude, target values ​​for agricultural film content, and various environmental variables obtained from satellite remote sensing are read from the Geographic Database (GDB) file, and invalid values ​​are filtered out. At the same time, the environmental variable raster data matching the spatial location of the sample points are read, and their spatial reference coordinate system and pixel resolution are unified. The data type is converted to floating-point format to ensure the consistency and comparability of subsequent model inputs.

[0095] Secondly, the sample data was proportionally divided into a training set (80%) and a validation set (20%), and the training set was further divided into an internal validation set (8:2) for model selection and performance evaluation. A random forest model was used for 50 iterations of training, with different random seeds and fixed model hyperparameters in each iteration to enhance the diversity and stability of the model training. After each training iteration, the RMSE and R2 scores on the internal validation set were calculated and recorded. Finally, the model with the highest internal validation R2 score was selected as the optimal random forest model and saved as a model file (DMHL_RF_Best.model).

[0096] A vectorized optimized Kriging interpolation model is constructed simultaneously. First, the distance matrix between sample points and the semi-variogram are calculated based on the spatial coordinates of the sample points and the target value. A spherical model is then used to fit the semi-variogram curve to obtain the variogram parameters. Second, a Kriging matrix A is constructed using the variogram and the distance between sample points, and its inversion is performed once to support subsequent multi-batch weight calculations. In the prediction stage, the pixels to be predicted are vectorized according to batches (batch_size=1000), and the distances and variogram values ​​between each pixel and the sample points are calculated in batches to construct the corresponding batch matrix b. The weights are then solved using matrix multiplication, and the prediction results are calculated.

[0097] For the random forest branch, the stacked raster data of environmental variables is reshaped into an input matrix of (n_pixels × n_features). Pixels containing NoData (any pixel with a feature value of -9999 is considered invalid) are filtered out, and the effective pixels are predicted using the optimal random forest model. The prediction results are then reshaped into a raster array. For the kriging branch, a vectorized kriging model is used to predict pixels within the template range in batches, generating a complete prediction array. This array is then converted to ArcGIS rasters using NumPyArrayToRaster, and the projection information is defined before saving it to the result FileGDB. Finally, based on the validation set, the prediction accuracy of the two models is compared, and the optimal raster for predicting agricultural film content is generated and stored in the result database. This method, by integrating machine learning optimization training with geostatistical spatial interpolation, effectively improves the accuracy and computational efficiency of characterizing the spatial heterogeneity of agricultural film residue, achieving high-precision and high-efficiency mapping of large-scale agricultural film residue hotspots.

[0098] Step 3: GPR inversion to determine the dielectric constant of residual agricultural film inside the soil

[0099] The inversion of the dielectric constant of ground-penetrating radar (GPR) can be transformed into the inversion of radar wave velocity. By obtaining the radar wave velocity and combining it with the two-way propagation time, the subsurface depth of the target medium can be further calculated. Based on this, a deep learning-based inversion objective function is constructed. Its goal is to use a trained nonlinear mapping function to convert the observed GPR amplitude data into corresponding radar wave velocity properties, thereby achieving the indirect inversion of the dielectric constant.

[0100]

[0101] In the formula, c is the speed of light. , is the relative permittivity of the residual agricultural film.

[0102]

[0103] In the formula, m represents the data tag (electromagnetic wave velocity). This represents the mapping function that converts observation data into radar wave velocity, where d represents the measured ground-penetrating radar data. This represents the set of parameters to be optimized (i.e., weights and biases).

[0104] First, a training dataset covering the common electromagnetic wave velocity range in ground-penetrating radar applications is constructed, selecting the electromagnetic wave velocity range of 0.054–0.212 m / ns, corresponding to a relative permittivity range of 2–30, and generating a velocity model containing 3–10 random layers. However, the generated electromagnetic wave velocity model is located in the spatial depth domain, while the amplitude data obtained from forward modeling is located in the temporal depth domain. To ensure that the neural network input and output remain consistent within the same depth domain, the spatial depth domain to temporal depth domain conversion needs to be achieved using the following formula, where the converted one-dimensional electromagnetic wave velocity serves as the output label for neural network training.

[0105]

[0106] Secondly, a two-dimensional time-domain finite-difference method was used for forward modeling to generate 9000 sets of radar amplitude data corresponding to one-dimensional electromagnetic wave velocities; the electromagnetic wave velocity was then converted into the relative permittivity. Finally, the relative permittivity distribution map of residual agricultural film in the soil was obtained.

[0107] Step 4: Construction of a spectral-dielectric joint inversion model based on multi-head attention mechanism

[0108] To achieve cross-modal nonlinear mapping of multispectral satellite characteristic variables to key parameters of ground-penetrating radar, a Transformer-based correlation model is constructed. The model construction process consists of four key stages: Collaborative sample construction; Feature preprocessing and embedding; Transformer encoding-decoding modeling; Loss function and training strategy.

[0109] Collaborative Sample Construction: Remote Sensing Pixel-Radar Dielectric Spatiotemporal Alignment

[0110] Point-to-pixel mapping is performed on multispectral satellite raster data and GPR survey line data. The satellite pixel location of each GPR sampling point is determined by bilinear interpolation. For each GPR sample point, its corresponding 11-dimensional satellite feature vector is extracted to construct a training sample, where the input features are... The output labels are

[0111] Preprocessing, normalization, and embedding representation of satellite variables

[0112] Univariate normalization was performed on each of the 11 input satellite feature variables. Standardize using the following formula:

[0113]

[0114] In the formula, , These are the mean and standard deviation, respectively, obtained from the training set. Then, the 11-dimensional normalized input vector of each pixel is mapped to an embedding vector of dimension d_model through linear projection, as shown in the following equation.

[0115]

[0116] In the formula, Projection matrix , For the original input Based on this, spatial location coding and temporal coding are introduced to separately address the pixel center coordinates. and observation time Constructing learnable location embeddings:

[0117]

[0118] The final Encoder input is:

[0119]

[0120] Transformer Model Construction: Mapping from 11-dimensional satellite variables to 2 GPR parameters

[0121] The model employs the standard Encoder–DecoderTransformer architecture. The Encoder consists of six stacked layers, each containing a multi-head self-attention mechanism, a feedforward network (FFN), and an Add&Norm residual connection structure. Its input tensor is... The formula for multi-head attention calculation is: The final encoded output is .

[0122] Satellite condition characteristics for obtaining Encoder output Next, let Query be two "query vectors" representing the prediction targets (A and T), used to drive the Decoder to extract from the satellite features of the Encoder. The Decoder uses these two learnable query vectors... As input, the two are first stacked into a query sequence. Inside the Decoder layer, First, it goes through a multi-head self-attention mechanism, using scaled dot product attention.

[0123]

[0124] Get the updated query representation This allows the two queries representing amplitude A and two-way time t to be correlated in the feature space. Subsequently, As a query, output to the Encoder Perform cross-attention by

[0125]

[0126] Conditional representation is obtained This representation explicitly integrates spectral, vegetation, surface structure, and temperature features from satellite imagery. Subsequently, Entering the feedforward network

[0127]

[0128] The final output of the Decoder is obtained through the residual connection layer of Add&Norm. In this output, the first row vector Used to generate radar amplitude, second row vector Used to generate two-way time. Finally, two consecutive predictions are obtained through independent linear mapping heads:

[0129]

[0130] in , At this point, Decoder has completed the process of generating conditions from satellite characteristic conditions to key radar parameters (amplitude and two-way travel time), and provides accurate upstream radar attribute inputs for the subsequent dielectric constant inversion model.

[0131] Loss function and training strategy:

[0132] The model training employs a weighted mean squared error loss function and introduces an L2 weight regularization term. For the There are n samples, and the sample weights are 1. The real label is The model predicts that The weighted mean square error of the two objectives is defined as follows:

[0133]

[0134] The overall loss function is:

[0135]

[0136] in, , For loss weighting coefficients, The regularization coefficient is . This represents the set of model parameters. After the model is trained, this Transformer correlation model can be applied to generate radar amplitude and two-way time prediction results at the same pixel scale.

[0137] Step 5: 3D pixel-level visualization output of residual agricultural film

[0138] Multidimensional visualization of results includes both two-dimensional and three-dimensional views. The two-dimensional view supports the overlay display of remote sensing base maps, agricultural film distribution maps, GPR survey line locations, and sampling point information, and can simultaneously display the original GPR profile, processed profile, and interpreted profile. Based on this, the remote sensing classification results are fused with the depth information obtained from GPR inversion to generate a three-dimensional spatial distribution map of residual agricultural film content in farmland. This supports viewing the profiles of main survey lines and connecting survey lines in any area, as well as interactive browsing of the three-dimensional rendering results.

[0139] Secondly, the system supports the export of results and the automatic generation of test reports. It can export results such as agricultural film distribution maps, burial depth maps, and risk assessment maps in commonly used raster and vector formats, and statistically output spatial distribution characteristics of residual agricultural film, including indicators such as residual area at different levels, average burial depth, and estimated total mass. Simultaneously, the system can automatically generate standardized test reports, covering project overview, technical methods, result charts, and comprehensive conclusions, providing intuitive support for farmland residual film surveys and management decisions.

[0140] Example 2:

[0141] This embodiment is based on the three-dimensional distribution inversion method of residual agricultural film based on satellite radar co-calibration described in Embodiment 1. Taking farmland in a county-level city A in a province in northern China as an example, the three-dimensional inversion of residual agricultural film is carried out. The specific process is as follows: (1) Overview of the study area and data preparation

[0142] A county-level city (A) in northern China was selected as a typical research area for agricultural film mulching. The total land area of ​​this area is approximately 7499 km², of which arable land accounts for 3363.28 km². The research scope is as follows: Figure 2As shown, the seasonal freezing period (October–November) after the autumn harvest in 2024 and before the spring planting of the following year (February) was selected as the key time period for monitoring and experimenting with residual agricultural film. After the autumn harvest, crops are harvested, farmland surfaces are exposed, and residual agricultural film mainly remains on the soil surface or is turned into shallow soil, making it an ideal window for large-scale remote sensing monitoring. Before spring planting, a large amount of residual agricultural film is covered by ice and snow, some frozen to the surface, and some embedded in the soil, providing favorable conditions for radar to conduct local three-dimensional detection of agricultural film distribution. By integrating multispectral remote sensing and radar detection data acquired from the same area at different times, complementary use of multi-source information is achieved, providing decision support for precision agriculture and the management of residual agricultural film.

[0143] like Figure 3 As shown, multispectral data from Sentinel-2 in October and November were selected. The Sentinel-2 satellite has a band range of 400-2400 nm, covering three characteristic bands of residual agricultural plastic film (1200 nm, 1730 nm, 2310 nm), with 13 bands (B1, B2, B3, B4, B5, B6, B7, B8, B8a, B9, B10, B11, B12). The revisit cycle for the two satellites working together in mid-latitude regions is 2-3 days, with a minimum spatial resolution of 20 m. The coverage area of ​​a single image is 290*290 km². A single image from one sampling point was selected and divided into... The block is a small area of ​​290*290m2.

[0144] The download and processing of Sentinel-2 remote sensing images were both implemented on the GEE cloud computing platform. The Sentinel-2 Level-2A surface reflectance (SR) data stored in GEE had already undergone radiometric calibration and atmospheric correction preprocessing. Therefore, in this embodiment, the preprocessing of Sentinel-2 remote sensing images mainly included cloud removal, mosaicking, and cropping. To characterize spectral and vegetation index features and fill in pixel missing values ​​caused by cloud pollution, the HANTS harmonic regression method was used in this case. Since the temporal variation of spectral reflectance is similar to a sine wave, trigonometric function fitting was performed on the time-series curves of mulched farmland on different reflectance bands and vegetation indices. The fitting time window was from the end of October to the end of November, and each band was considered as a time-related curve. The function, denoted as f(t):

[0145]

[0146] In the formula, the independent variable t is a date in a year, expressed as a decimal between 0 (January 1st) and 1 (December 31st); It is the intercept term; It is the order of the harmonic sequence; These are the coefficients of the cosine terms; These are the coefficients of the sine term; It is the angular frequency. In the above harmonic regression formula, the parameter... and Optimization is needed to balance fitting accuracy and prevent overfitting. This study evaluates the change in the root mean square error (RMSE) of the fit and selects... and As the optimal parameter combination, the final harmonic regression formula is obtained:

[0147]

[0148] Based on this, principal component analysis (PCA) was used to reduce the dimensionality of the harmonic regression curves fitted in each small region, extracting characteristic variables such as absorption peak depth, interband ratio, and derivative spectrum; subsequently, a single 290×290m... 2 The pixel features within a small area are averaged to obtain representative spectral feature values ​​for that area. Simultaneously, the advanced plastic greenhouse index and the improved plastic greenhouse index are calculated, and the presence of residual agricultural film in different cultivated land soils is determined by comparing the spectral features of mixed pixels with those of pure pixels containing residual agricultural film. Furthermore, topographic factors (elevation, slope, aspect) and climatic factors (precipitation, surface temperature, temperature-vegetation drought index) of the study area are obtained and combined with multispectral features to form a satellite remote sensing variable set. It is used to invert the content distribution of residual agricultural film on the ground surface.

[0149] The ground-penetrating radar used was the GS9000 multi-channel 3D ground-penetrating radar. The equipment was equipped with 35 vertically polarized (VV) antennas and 15 horizontally polarized (HH) antennas, forming a dual-polarized multi-channel array with a minimum antenna spacing of 2.5 cm. The GS9000 array module scanned the frequency band from 500 MHz to 3 GHz, with a scan frequency of 27,500 times / meter and a time window of 0 to 35 ns. The experiment preset the center frequency to 500 MHz, and the antennas moved along the polarization direction. The spacing between the survey lines was designed first using the antenna radiation range formula. The antenna radiation can be approximated by an elliptical model, as shown in the following formula:

[0150]

[0151]

[0152] In the formula: The wavelength of electromagnetic waves; The depth of the projection surface; It is the relative permittivity; The radius of the major axis of the ellipse; Let be the radius of the minor axis of the ellipse.

[0153] Under the estimated conditions, the pulse width of the 500MHz antenna is approximately 2ns, corresponding to a wavelength of 0.2m; taking the relative permittivity of dry farmland soil as 6 and the detection depth as 1m, the calculated values ​​are... Based on this, the survey line spacing was set to 0.8m. In each 290×290m... 2 Within a small area, radar amplitude and two-way time information for different media are collected along the survey line, ultimately obtaining waveform diagrams and envelope diagrams, such as... Figure 4 As shown.

[0154] For soil sampling of residual agricultural film, 100 typical small areas were selected as sampling points. A two-stage density separation method was used to extract residual agricultural film from the soil samples. The specific steps were as follows: 5g of dried soil sample was placed in a 500mL beaker or Erlenmeyer flask and digested with H2O2 in a constant temperature shaking incubator; the mixture was then filtered through an 11µm filter membrane. The particles on the filter membrane were then treated with 100mL of saturated NaCl solution (…). Rinse into a 250mL Erlenmeyer flask. After shaking for 1 hour and settling for 12 hours, the supernatant was filtered through an 11µm filter membrane via a vacuum filter. The remaining material was further added to a saturated NaBr solution, and the shaking and settling process was repeated, with the supernatant being filtered again. Finally, the filter membrane was placed in a clean glass petri dish and allowed to air dry at ambient temperature.

[0155] Quantitative characterization of residual agricultural film was performed using an optical stereomicroscope. The residual agricultural film particles were photographed, counted, and their volume estimated. Their abundance in soil samples was calculated using the following formula:

[0156]

[0157] In the formula, A represents the abundance of residual agricultural film in the soil, in mm. 3 / L; N is the volume of a single residual agricultural film, in mm. 3 V represents the total volume of filtered wastewater, in liters (L).

[0158] (2) Satellite remote sensing inversion of agricultural film content on the land surface

[0159] First, prepare the environment. On a host with ArcGIS Pro, set up a Python interpreter environment and install the necessary libraries (numpy, pandas, scipy, scikit-learn, joblib, arcpy). Read raster data corresponding to each environmental variable from the geodatabase. Use the first raster as a template to unify the spatial reference coordinate system and pixel resolution, and convert all rasters to 32-bit floating-point type before stacking them into a 3D data volume. Then, divide the sample data into training and validation sets in an 8:2 ratio. Based on the training set, start 50 iterations of optimization training of the random forest model. Each iteration uses a different random seed to build a model containing 300 decision trees, and the R... 2 The model with the highest performance index was selected as the optimal predictor. Simultaneously, the Kriging interpolation model was trained in parallel, calculating the distance matrix and semi-variogram between samples using vectorization, fitting the parameters of the spherical variogram, and pre-calculating the inverse matrix of the equation system. In the model application phase, the prediction accuracy of the random forest model and the Kriging model was evaluated using a validation set. Based on this, the optimal random forest model was applied to a raster stack for cell-level prediction, and batch vectorized computation was used to implement Kriging spatial interpolation. Finally, the prediction results of both models were saved to a geographic database, and a detailed report including accuracy evaluation indicators was generated. The results show that both random forest and Kriging interpolation can effectively predict residual agricultural film, with the random forest reaching its optimal performance on the 50th training iteration. , Fitting pairs, for example Figure 5 As shown; Kriging interpolation in The value above is 4.4414. The value is 0.7121. In the global residual analysis, the random forest has a relatively uniform distribution of positive and negative errors across different residual agricultural film content ranges, and its overall performance is better than that of Kriging interpolation.

[0160] (3) Dielectric constant of residual agricultural film in soil obtained by GPR inversion

[0161] In the process of inverting the dielectric constant of residual agricultural film in soil using ground penetrating radar, this invention does not adopt the traditional analytical inversion method based on a single hyperbola fitting. Instead, it constructs a data-driven inversion framework that integrates "physical forward modeling constraints and deep learning inversion" to enhance the stability and accuracy of the inversion results under complex soil medium conditions and weak scattering target scenarios.

[0162] Since the energy of direct waves from ground-penetrating radar (GPR) is significantly higher than that of reflected echoes from residual agricultural film underground, the graph-based window function is first improved based on the distribution characteristics of the direct waves along the time axis. This improved window function is then used to filter the original GPR data, effectively suppressing direct wave interference while ensuring the energy integrity of the effective reflected waves within the time window. The filtered GPR amplitude data is further subjected to uniform amplitude normalization and time alignment, and used as a one-dimensional amplitude sequence as input features for the deep learning model.

[0163] According to the method described in this invention, since the ground-penetrating radar (GPR) observation data is located in the time-depth domain, to ensure scale consistency between the neural network input and output, the velocity model in the spatial depth domain needs to be converted into a time-depth domain velocity model, so that the amplitude sequence of the network input and the output radar wave velocity form a one-to-one correspondence on the time scale. Based on this, the two-dimensional finite-difference time-domain (FDTD) method is used to perform forward modeling on each velocity model, generating the corresponding GPR amplitude response. Finally, approximately 9000 sets of "amplitude-wave velocity" sample pairs are constructed as the training and testing sets for the neural network.

[0164] After the convolutional neural network is trained, four sets of one-dimensional velocity model samples are randomly selected from the test set. The corresponding ground-penetrating radar amplitude data are preprocessed and input into the trained neural network model, which outputs a continuous one-dimensional radar wave velocity profile to predict the wave velocity value. All values ​​were above 95%, demonstrating good accuracy. Subsequently, the predicted electromagnetic wave velocity was converted into a relative permittivity distribution map of different media at different depths. By comparing this map with the permittivity of residual agricultural film measured in the field, the spatial distribution of residual agricultural film in the soil at different underground depths can be obtained, as illustrated below. Figure 6 As shown.

[0165] (4) Construction of a spectral-dielectric joint inversion model based on multi-head attention mechanism

[0166] Multispectral satellite data and ground-penetrating radar (GPR) data differ significantly in observation scale, data structure, and physical mechanism. There is no stable analytical mapping relationship between them, and GPR data exhibits highly nonlinear characteristics due to the combined influence of multiple factors such as soil type, moisture content, and the burial state of agricultural film. To effectively constrain the response of underground radar using surface remote sensing information, this invention constructs a Transformer-based multispectral satellite-GPR parameter correlation model to establish a conditional mapping relationship between satellite environmental characteristics and key radar observation parameters.

[0167] Using GPR survey line sampling points as a reference, spatial registration and interpolation methods are used to establish a one-to-one correspondence between them and multispectral satellite pixels, constructing a collaborative sample set containing 11-dimensional satellite environmental features and corresponding GPR amplitudes and two-way times. The input features are standardized and mapped to a high-dimensional embedding space through linear projection. At the same time, learnable spatial location codes are introduced to preserve the spatial constraint information of the remote sensing data.

[0168] The correlation model adopts a Transformer architecture with an Encoder–Decoder structure. The Encoder uses a multi-head self-attention mechanism to globally model the dependencies of multi-source satellite environmental features, extracting conditional feature representations that constrain the radar response. The Decoder uses amplitude and two-way time as query targets, and selectively extracts relevant information from the Encoder output through self-attention and cross-attention mechanisms to generate corresponding continuous radar parameter predictions. The initial hyperparameters are as follows: , , Training uses the AdamW optimizer (weight_decay). ), initial learning rate Ir A linear warmup of 2000 steps is used, followed by cosine decay scheduling; the batch size can be set to 512 based on the available GPU memory, and global norm clipping is used. To ensure stable training.

[0169] For data preparation, first calculate and save the mean and standard deviation of each input feature and label on the training set (scaler), then standardize the input features and labels (zero-mean unit-variance). During training, a small amount of Gaussian noise can be added to the input for enhancement. and The samples were stratified by time and scenario into training, validation, and testing, with a typical ratio of 80:10:10.

[0170] Model training employs a standard supervised learning process. For each training batch, the model first performs forward computation to obtain predicted values ​​for radar amplitude and two-way time. and Then, a weighted loss function is calculated based on the true labels. The network gradient is updated via backpropagation. Gradients are pruned before parameter updates to suppress gradient explosion, followed by optimizer parameter updates and learning rate scheduling. The model is trained iteratively in epochs. After each epoch, the root mean square error (RMSE) of amplitude and two-way time, RMSE(T), and the overall coefficient of determination are calculated on the validation set. The algorithm employs an early stopping strategy based on validation performance. Training terminates when the validation metrics no longer improve after several consecutive rounds, while simultaneously saving the optimal model parameters, optimizer state, and feature normalizer from the validation set. The initial training rounds are set to a finite number; if the validation error continues to decrease, the training period can be extended appropriately. To ensure the reproducibility of the algorithm results, a fixed random seed (seed = 42) and uniform computational environment settings are implemented before training begins. After training, the optimal model and its corresponding standardized parameters are exported to maintain consistency in input-output processing during the inference phase.

[0171] The final model uses weighted mean square error as the training objective and introduces regularization constraints to improve generalization ability. After training, it can predict radar amplitude and two-way time parameters of the target area, relying only on multispectral satellites and environmental variables, even in the absence of measured GPR data. These parameters can then be used as input to the dielectric constant inversion model, realizing cross-scale transfer from surface remote sensing information to subsurface physical parameters.

[0172] (5) Three-dimensional inversion of residual film and verification of results

[0173] The distribution of agricultural film content in hot zones on the land surface was initially obtained through satellite remote sensing inversion. Then, using a spectral-dielectric joint inversion model, environmental variables from satellite remote sensing were mapped to the amplitude and two-way travel time of the Gaussian spectral density (GPR). Subsequently, a model for inverting dielectric constant and depth using GPR was used to obtain a three-dimensional spatial distribution model of residual agricultural film in farmland. This model generated profiles of main survey lines and connecting survey lines supporting arbitrary regions, allowing for the viewing of three-dimensional rendering images, such as... Figure 7 As shown in the figure, the inversion results were verified using field sampling data. Within the depth range of 0-30 cm, the coefficient of determination (R2) between the model inversion values ​​and the measured values ​​reached 0.83, and the root mean square error (RMSE) was 0.12 g / kg. This confirms the effectiveness and reliability of this method for high-precision and visual calibration of agricultural film residue in three-dimensional space.

Claims

1. A residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration, characterized in that, The method comprises the following steps: S1, satellite remote sensing and ground penetrating radar data extraction: extracting feature factors from multispectral satellite images, the feature factors including terrain features, climate features, spectral principal component variables and plastic indexes, which together constitute a satellite remote sensing environmental variable set; obtaining A-scan and B-scan signals of ground penetrating radar along the survey line, calculating electromagnetic wave phase velocity, target depth and relative permittivity through hyperbolic fitting, performing anti-gain and geometric compensation on the echo envelope signal, and calculating the attenuation coefficient by means of sliding window logarithmic linear fitting, constructing depth-sensitive dielectric features based on the depth weighting principle, and forming a radar feature subset; S2, satellite remote sensing inversion of surface residual agricultural film hot area: based on the satellite remote sensing environmental variable set, using random forest model and Kriging interpolation model to spatially predict and map the content of surface residual agricultural film, and generating a predicted grid of surface agricultural film content; S3, ground penetrating radar inversion of dielectric constant of residual agricultural film in soil: constructing and training a deep learning model, inverting the ground penetrating radar amplitude data into electromagnetic wave velocity, and then calculating the relative permittivity distribution of residual agricultural film in soil as the underground dielectric constant inversion result; S4, spectral and dielectric combined inversion model construction based on multi-head attention mechanism: constructing a Transformer encoder-decoder model, taking satellite remote sensing environmental variables as input, learning the nonlinear mapping relationship between them and ground penetrating radar amplitude and two-way time through multi-head attention mechanism, realizing cross-modal combined inversion from satellite features to radar parameters; S5, three-dimensional residual agricultural film pixel-level visualization output: fusing the predicted grid of surface agricultural film content and the underground dielectric constant inversion result to generate and output a three-dimensional spatial distribution map of residual agricultural film.

2. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, In the step S1, the terrain features include elevation, slope and aspect; the climate features include precipitation, land surface temperature and temperature vegetation drought index; the spectral principal component variables include absorption peak depth, band ratio and derivative spectral features; the plastic indexes include advanced plastic greenhouse index (APGI) and modified plastic greenhouse index (MPGI); The formula for calculating the advanced plastic greenhouse index (APGI) is: ； wherein, is the surface reflectance in the coastal band, is the surface reflectance in the red band, is the surface reflectance in the near infrared band, is the surface reflectance in the short wave infrared second band; The formula for calculating the modified plastic greenhouse index (MPGI) is: ； wherein is the surface reflectance in the ultra-blue band, is the surface reflectance in the short-wave infrared first band; The formula for the temperature vegetation drought index (TVDI) is: ； wherein, is the surface temperature of the target pixel; and are the minimum and maximum surface temperatures of the wet edge in the large triangle feature space, respectively.

3. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, The step S1 further comprises a step of preprocessing ground penetrating radar amplitude data: using a graph-based window function to filter radar raw data to remove direct wave interference, and the formula of the graph-based window function is: ； wherein is an index of a discrete sample point of the window function, M = (N - 1) / 2, N is a length of the window function, is a constant between 0 and 1.

4. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, The step S2 comprises: S21, reading sample data containing geographic coordinates, agricultural film content target values and multiple environmental variables from a geographic database, and reading corresponding environmental variable grid data, filtering invalid values, unifying spatial reference and resolution, and converting data types to form a model input data set; S22, dividing the sample data into a training set and a validation set, training the random forest model on the training set multiple times, and selecting the model with the highest coefficient of determination on the internal validation set as the optimal model; S23, synchronously construct a vectorized optimized Kriging interpolation model, calculate a semi-variogram based on sample points, and fit a spherical model to batch calculate prediction weights through matrix operation; S24, respectively use the optimal random forest model and the vectorized Kriging interpolation model to predict the environmental variable grid to obtain two types of prediction results; compare the accuracy of the two types of prediction results based on the validation set, and generate a surface agricultural film content prediction grid and store it in a geographic database.

5. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, The step S3 comprises: S31, a training data set covering the common electromagnetic wave velocity range in ground penetrating radar application is constructed, the electromagnetic wave velocity interval of 0.054-0.212 m / ns is selected, the corresponding relative permittivity range is 2-30, and a velocity model containing 3-10 random layers is generated; S32, the velocity model in the spatial depth domain is converted into a time-depth domain model, and the formula is: ； Wherein, t is the double-path propagation time of the ground penetrating radar electromagnetic wave from emission to target reflection and return to the receiving end, d is the one-way propagation distance of the electromagnetic wave in the medium from the antenna to the reflection interface, and v is the propagation velocity of the electromagnetic wave in the corresponding medium; S33, a two-dimensional time domain finite difference method is used to forward simulate the converted time-depth domain velocity model to generate corresponding ground penetrating radar amplitude data; S34, an inversion target function based on deep learning is constructed, and a nonlinear mapping function obtained by training is used to convert the observed ground penetrating radar amplitude data into corresponding radar wave velocity attributes, thereby realizing indirect inversion of the dielectric constant, and the formula is: ； wherein is the phase velocity of the electromagnetic wave in the target medium, c is the speed of light , is the relative permittivity of the residual agricultural film; ； where m is the data tag, i.e. the electromagnetic wave velocity, is the mapping function to convert the observed data into radar wave velocity, d is the measured ground penetrating radar data, is the set of parameters to be optimized.

6. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, The step S4 comprises: S41, cooperatively sample construction: point-pixel mapping is performed on the multispectral satellite grid data and the ground penetrating radar survey line data, the satellite pixel position where each ground penetrating radar sampling point is located is determined through bilinear interpolation, for each ground penetrating radar sampling point, an 11-dimensional satellite feature vector corresponding to the ground penetrating radar sampling point is extracted, and a training sample is constructed, wherein the input features are , and the output label is ; S42, preprocessing, normalization and embedding representation of satellite variables: The input 11 satellite characteristic variables are respectively processed by single variable normalization, and each variable Normalization is performed according to the following formula: ； wherein, , are the mean and standard deviation, respectively, calculated from the training set. The 11-dimensional normalized input vector of each pixel is mapped to an embedding vector with a dimension of d_model through linear projection, and the formula is: ； wherein is a projection matrix , is the original input ; Introduce spatial position encoding and temporal encoding, respectively, to the image element center coordinates and observation time Construct a learnable position embedding: ； The final encoder input is: ； S43, construction of a Transformer model: The Transformer model adopts an encoder-decoder architecture; the encoder is stacked by multiple layers, each layer containing a multi-head self-attention mechanism, a feedforward network and a residual connection structure, for encoding input features and outputting satellite condition features ​ The decoder takes two learnable query vectors As input, it is updated by a self-attention mechanism: ； where, are learnable linear projection parameters in the attention mechanism, map the input representation to query, key and value vectors, respectively, is the feature dimension of the key vectors, as the normalization scale in the scaled dot-product attention to control the dot-product magnitude, stabilize the softmax distribution; get updated query representation , as query, encoder output perform cross-attention computation: ； Obtaining conditioned representation ; into the feedforward network: ； where x denotes each row vector of C, and is the learnable weight matrix of the first linear transformation; is the learnable weight matrix of the second linear transformation; and are the learnable bias vectors of the first and second linear transformations, respectively. and the final output of the decoder is obtained through residual connection and layer normalization ; Two continuous prediction values are obtained through independent linear mapping heads: ； wherein , are two row vectors of matrix , respectively, corresponding to two queries (Q ) updated task-specific representations; , denote the learnable weights of the two regression heads, respectively; , denote the learnable bias terms of the two regression heads, respectively; , denote the two consecutive prediction values output by the two independent linear regression heads, respectively; S44, loss function and training strategy: The model training uses a weighted mean squared error loss function and introduces an L2 weight regularization term. For the There are n samples, and the sample weights are 1. The real label is The model predicts that The weighted mean square error of the two objectives is defined as follows: ； wherein, the weighted mean squared error for the target variable A, the weighted mean squared error for the target variable T, N being the total number of samples participating in the loss calculation; The overall loss function is: ； wherein, , is a loss weight coefficient, is a regularization coefficient, denotes a set of model parameters.

7. The residual agricultural film three-dimensional distribution inversion method based on satellite radar cooperative calibration according to claim 1, characterized in that, The step S5 comprises: S51, a visualization result containing a two-dimensional view and a three-dimensional view is generated; the two-dimensional view supports the superimposed display of remote sensing base maps, agricultural film distribution maps, ground penetrating radar survey line positions and sampling point information, and supports the synchronous display of ground penetrating radar original profiles, processed profiles and interpretation result profiles; the three-dimensional view fuses satellite remote sensing classification results and depth information obtained by ground penetrating radar inversion to generate a three-dimensional spatial distribution map of agricultural residual agricultural film content, and supports profile viewing of any regional main survey line and tie survey line and interactive browsing of three-dimensional rendering results; S52, the agricultural film distribution map, the burial depth map and the risk assessment map are exported in raster or vector format; the spatial distribution characteristic information of residual agricultural film is counted and output, the information includes residual area of different grades, average burial depth and estimated total mass; a standardized detection report is automatically generated, and the report content covers project overview, technical method, result chart and comprehensive conclusion.

Images