Exploration surveying and mapping method based on remote sensing and surveying and mapping instrument thereof

Abstract

The invention relates to the technical field of geophysical exploration, in particular to an exploration surveying and mapping method based on remote sensing and a surveying and mapping instrument thereof. The exploration surveying and mapping method based on remote sensing comprises the following steps: step 1, multi-platform collaborative data dynamic acquisition; step 2, multi-modal data fusion and feature enhancement: vegetation noise filtering is performed on the laser radar point cloud data, and a surface terrain model under vegetation coverage is reconstructed; seismic wave speed data and ground penetrating radar dielectric constant data are combined to construct an underground medium physical property combined distribution model; associating the reflectivity characteristics of the multispectral image with the underground physical property parameters, and identifying a geological abnormal region; and step 3, hidden geologic body dynamic modeling and risk assessment. By means of the method, real-time monitoring and risk assessment of the complex geologic body can be achieved, and the geological disaster prediction and prevention and control capacity is greatly improved.

Description

Technical Field

[0001] The present invention relates to the field of geophysical exploration technology, and in particular to a remote sensing-based surveying and mapping method and a surveying and mapping instrument thereof. Background Art

[0002] In karst-developed areas, hidden cavities within 5-10 meters below the surface account for 83% of geological disaster sources, but their spatial distribution is affected by multiple factors such as lithology, hydrology, and vegetation coverage. The existing single technical means have a detection miss rate of up to 40%-60%. In addition, 65% of road collapse accidents caused by leakage in urban underground pipe networks occurred in sections with pipe diameters less than 0.8 meters and burial depths greater than 4 meters. Traditional detection methods are unable to achieve large-scale rapid screening. In the fields of underground cavity detection in karst areas, identification of hidden dangers in urban underground pipe networks, and early warning of shallow geological disasters in vegetation-covered areas, existing technologies have the following significant defects:

[0003] Optical satellite imagery is limited by vegetation cover and weather conditions, making it ineffective for capturing topographic information beneath the vegetation canopy and shallow geological depths. While synthetic aperture radar (SAR) has some penetration, its resolution (>5 meters) for small-scale geological structures within 10 meters below the surface (such as karst conduits with a diameter of less than 3 meters) is insufficient, and it cannot distinguish differences in the physical properties of the medium. While multispectral satellite data can capture the spectral characteristics of surface lithologies, atmospheric attenuation significantly reduces the signal-to-noise ratio in the shortwave infrared band (1500-2500nm).

[0004] Traditional seismic exploration relies on manually deployed dense detector arrays (spacing <10 meters), which takes up to 2-3 weeks to deploy in complex terrain areas. In addition, data inversion relies on empirical velocity models, making it difficult to accurately characterize shallow heterogeneous geological structures (such as fracture zones and water-filled cavities). Summary of the Invention

[0005] Based on the above objectives, the present invention provides a surveying and mapping method based on remote sensing and a surveying and mapping instrument thereof, wherein a surveying and mapping method based on remote sensing includes the following steps:

[0006] Step 1: Dynamic data collection through multi-platform collaboration:

[0007] Using a rotary-wing drone equipped with a penetrating lidar and a multispectral camera, a 3D scan of the vegetation-covered area is performed, simultaneously acquiring vegetation penetration point cloud data and multispectral image data.

[0008] Deploy seismic wave sensor arrays and mobile ground-penetrating radar nodes in the target area to collect seismic wave propagation velocity data and dielectric constant distribution data of the underground medium, respectively;

[0009] Collect rainfall, soil moisture and surface temperature data of the target area in real time;

[0010] Step 2: Multimodal data fusion and feature enhancement:

[0011] filtering vegetation noise from the lidar point cloud data to reconstruct a surface terrain model under vegetation cover;

[0012] Combine seismic wave velocity data with ground-penetrating radar dielectric constant data to construct a joint distribution model of underground medium physical properties;

[0013] Correlate multispectral image reflectance characteristics with underground physical parameters to identify geological anomalies;

[0014] Step 3: Dynamic modeling and risk assessment of hidden geological bodies:

[0015] Generate a 3D geological body probability model based on the fused multi-source data;

[0016] Dynamically integrate new sensor data and meteorological factors to update model confidence;

[0017] Output multi-scale risk assessment results, including regional risk heat maps and critical collapse warning parameters.

[0018] Preferably, the three-dimensional scanning flight path planning in step 1 includes:

[0019] Calculate vegetation index through multispectral imagery and dynamically divide the scanned area into grades based on vegetation index;

[0020] For areas with high vegetation coverage, a spiral intensive route is adopted, and the flight altitude is adjusted according to the calculated results of the vegetation index:

[0021] When the vegetation index is higher than a first threshold, the flight altitude is lowered to a preset minimum safe altitude;

[0022] The sampling density of the ground penetrating radar mobile node is increased synchronously, and the sampling interval is inversely proportional to the flight altitude.

[0023] Preferably, the acquisition of seismic wave propagation velocity data in step 1 includes:

[0024] A triangular grid sensor array is deployed, and the node spacing is set according to the complexity of the terrain:

[0025] The spacing is 20 meters in flat areas and 10 meters in hilly areas;

[0026] The elastic waves are excited by artificial seismic sources, and the arrival time difference between the direct wave and the reflected wave is recorded;

[0027] The velocity distribution of the underground medium is inverted based on the velocity tomography algorithm. The constraints in the velocity inversion process include:

[0028] Dielectric density range mapped by GPR dielectric constant data;

[0029] Indicators of shallow rock mass integrity as reflected by lidar terrain curvature data.

[0030] Preferably, the vegetation noise filtering in step 2 includes:

[0031] Calculate the reflection intensity distribution histogram of the laser point cloud and set a dynamic threshold to distinguish vegetation from surface point clouds:

[0032] The reflection intensity threshold is dynamically modified according to the atmospheric transmittance during the scanning period, and the atmospheric transmittance is obtained by inverting the water vapor absorption band of the multispectral image;

[0033] Perform local surface fitting on the initially screened surface point cloud and remove abnormal points that deviate from the fitted surface by more than an adaptive threshold. The adaptive threshold is calculated based on the terrain roughness:

[0034] The roughness index is determined by combining the standard deviation of the point cloud elevation and the neighborhood slope.

[0035] Preferably, the construction of the physical property joint distribution model in step 2 includes:

[0036] Establish a joint objective function of seismic wave travel time residual and ground penetrating radar reflection amplitude;

[0037] The regional geological map is introduced to constrain the lithologic categories of each grid point, and the lithologic categories are converted into prior distributions of physical property parameters;

[0038] An improved particle swarm optimization algorithm is used to solve the objective function. During the optimization process, the weight coefficients of seismic data and ground penetrating radar data are dynamically adjusted. The weight coefficients are calculated based on the following factors:

[0039] Seismic wave signal-to-noise ratio;

[0040] GPR sampling density;

[0041] The influence factor of surface temperature on the propagation speed of electromagnetic waves.

[0042] Preferably, the generation of the three-dimensional geological body probability model in step 3 includes:

[0043] Mapping seismic wave velocity, dielectric constant, terrain curvature and spectral anomaly index to a unified spatial grid;

[0044] The non-parametric Bayesian method is used to calculate the probability value of each grid point being a gap, where the prior distribution is constructed based on the following data:

[0045] Karst development patterns in historical disaster databases;

[0046] Groundwater level gradient estimated from real-time rainfall data;

[0047] Optimize probability distribution surfaces via Markov chain Monte Carlo sampling.

[0048] Preferably, the model confidence update in step 3 includes:

[0049] When the deviation between the newly added seismic velocity data and the model prediction exceeds the historical fluctuation range, the following update mechanism is triggered:

[0050] Calculate the correlation coefficient between the deviation value and meteorological data (rainfall, soil moisture);

[0051] Adjust the weight of each data source in the physical property joint distribution model according to the correlation coefficient;

[0052] The updated model outputs include:

[0053] Correction value of probability distribution of geological body boundary;

[0054] Heatmap of the spatial distribution of model uncertainty.

[0055] Preferably, the generation of the critical collapse warning parameters in step 3 includes:

[0056] Extracting three-dimensional geometric parameters of potential cavities based on a probabilistic model, including roof thickness, span ratio, and surrounding rock fracture development index;

[0057] A coupled model of rainfall infiltration, groundwater level, and mechanical stability was established to calculate the critical collapse conditions:

[0058] Input parameters include real-time rainfall intensity, soil permeability coefficient and rock mass tensile strength;

[0059] The permeability coefficient is obtained by inverting the mapping relationship between the dielectric constant of ground penetrating radar and soil moisture content;

[0060] The output parameters include the curve of safety factor changing over time and the threshold value of warning level division.

[0061] Accordingly, an embodiment of the present invention further provides a remote sensing-based surveying and mapping instrument for executing a remote sensing-based surveying and mapping method according to an embodiment of the present invention, comprising:

[0062] Data acquisition module, data processing module, dynamic modeling module and risk assessment module, each module is connected through a high-speed data bus;

[0063] The data acquisition module includes:

[0064] The rotary-wing UAV unit integrates a penetrating laser radar, a multispectral imaging unit, and a navigation control unit. The laser radar's transmitting wavelength is configured to penetrate vegetation, and the multispectral imaging unit covers the visible to short-wave infrared band.

[0065] A ground sensing unit, comprising a seismic sensor array deployed in a triangular grid topology and a ground-penetrating radar mobile node equipped with an autonomous navigation chassis;

[0066] The meteorological monitoring unit, including a distributed rain gauge, soil moisture probe, and infrared thermometer, communicates with the main control unit via the LoRa wireless protocol;

[0067] The data processing module includes:

[0068] Point cloud filtering unit, connected to the laser radar data output terminal, with built-in reflection intensity analysis subunit and elevation surface fitting subunit;

[0069] A joint inversion unit receives output data from the seismic wave sensor array and the ground penetrating radar node and is equipped with a parallel computing accelerator;

[0070] The spectral correlation unit is connected to the data interface of the multispectral imaging unit and has a built-in band ratio calculator and anomaly feature extractor;

[0071] The dynamic modeling module includes:

[0072] The three-dimensional gridding unit receives the physical parameter matrix output by the data processing module and generates a voxel grid in a unified spatial coordinate system;

[0073] The probability calculation unit integrates the Bayesian update engine and Monte Carlo sampler, and is connected to the real-time data stream of the meteorological monitoring unit;

[0074] The risk assessment module includes:

[0075] The heat map generation unit receives the output data of the probability model and interacts with the geographic information system database;

[0076] Mechanical coupling analysis unit, with built-in finite element solver, and the output end connected to the early warning signal transmitter.

[0077] Preferably, in the data acquisition module, the rotor UAV unit and the ground sensor unit establish a collaborative control link, specifically including:

[0078] The UAV navigation control unit receives the position data of the ground penetrating radar mobile node in real time and dynamically adjusts the flight path to maintain complementary spatial coverage between the two;

[0079] The scanning frequency of the lidar and the sampling timing of the ground penetrating radar are aligned through hardware synchronization signals, with a time error of ≤1 millisecond;

[0080] The joint inversion unit of the data processing module is provided with a weight dynamic allocator, the input end of which is connected to:

[0081] The output of the seismic wave signal-to-noise ratio detector, where the signal-to-noise ratio detection is based on the power spectral density ratio of the direct wave to the background noise;

[0082] Ground penetrating radar sampling density counter, which counts the number of effective detection points within a unit area;

[0083] Surface temperature compensator, which corrects the electromagnetic wave propagation velocity model based on infrared thermometer data;

[0084] A spatial interpolation acceleration card is configured in the three-dimensional gridding unit of the dynamic modeling module, and its interpolation algorithm selection logic includes:

[0085] When the meteorological monitoring unit detects that the rainfall is greater than the threshold, it automatically switches to the Kriging interpolation mode with rainfall drift term;

[0086] In vegetation-covered areas, an anisotropic interpolation kernel function based on lidar terrain curvature is enabled;

[0087] A two-way feedback channel is provided between the mechanical coupling analysis unit and the dynamic modeling module of the risk assessment module:

[0088] The forward channel transmits the geometric parameters of the geological body to the finite element solver;

[0089] The reverse channel transmits the mechanical stability calculation results back to the probability calculation unit, triggering the model confidence correction.

[0090] Beneficial effects of the present invention:

[0091] This multi-platform collaborative data collection, fusion, and analysis approach enables accurate and efficient geological surveys. First, the combined operation of rotary-wing drones and ground sensors ensures comprehensive coverage of both surface and subsurface data, making data collection more precise. Second, the fusion of multimodal data (such as lidar point clouds, seismic wave data, ground-penetrating radar data, and multispectral imagery) enhances information accuracy, enabling a comprehensive analysis of underground structures from multiple dimensions. Dynamically updated geological probabilistic models and risk assessment results enable this system to monitor potential geological risks in real time and issue timely warnings, providing important decision-making support for related fields (such as civil engineering and mineral exploration), improving exploration efficiency, and reducing disaster risks. This approach helps reduce human error, enhances surveying accuracy, and offers strong real-time and adaptability, making it a beneficial innovation in modern surveying and mapping. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only for the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0093] Figure 1 is a flow chart of the steps of the method of the present invention;

[0094] Figure 2 A flowchart showing the steps for establishing a rainfall infiltration-groundwater level-mechanical stability coupling model and calculating critical collapse conditions for the method of the present invention;

[0095] Figure 3 It is a structural block diagram of the system of the present invention. DETAILED DESCRIPTION

[0096] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments. It is also noted that, to provide a more detailed description, the following embodiments are best and preferred embodiments, and those skilled in the art may employ alternative methods for implementing certain known technologies. Furthermore, the accompanying drawings are intended only to provide a more detailed description of the embodiments and are not intended to limit the present invention.

[0097] See Figure 1-Figure 3 An embodiment of the present invention provides a remote sensing-based survey and mapping method. In step 1, a three-dimensional scan is performed in a vegetation-covered area using a rotary-wing drone equipped with a penetrating laser radar and a multispectral camera to obtain point cloud data and image data in real time. These data will be used for further processing and analysis. The deployment of seismic wave sensor arrays and mobile ground-penetrating radar nodes can respectively collect seismic wave propagation velocity data and dielectric constant distribution data of the underground medium. The data collected by these sensors can help understand the characteristics of the underground structure. Meteorological monitoring equipment (such as rain gauges, soil moisture probes, and infrared thermometers) is responsible for obtaining surface environmental data, such as rainfall, soil moisture, and surface temperature.

[0098] In step 2, vegetation noise is first filtered from the lidar point cloud data. This process effectively removes noise caused by vegetation interference by analyzing reflection intensity and atmospheric transmittance to correct the reflection intensity threshold, thereby accurately reconstructing the surface terrain model under vegetation cover. Subsequently, by combining seismic wave velocity data with dielectric constant data from ground-penetrating radar, a joint distribution model of underground medium physical properties is constructed, further enhancing the understanding of underground physical parameters. At the same time, the multispectral image reflectance characteristics are correlated with underground physical parameters to identify potential geological anomalies, helping surveyors identify risky geological areas.

[0099] In step 3, a 3D probabilistic model of the geological volume is generated by fusing multi-source data. This model is constructed based on multimodal data such as lidar point clouds, seismic wave propagation velocity, and ground-penetrating radar dielectric constant, forming a multidimensional, comprehensive model of the underground geological structure. With the real-time integration of new sensor data and meteorological factors, the model's confidence level is updated, ensuring that the assessment results remain current and accurate. By generating regional risk heat maps and critical collapse warning parameters, the system can help users monitor geological risks in real time, provide timely warning information, and prevent potential geological disasters.

[0100] This multi-platform collaborative data collection, fusion, and analysis approach enables accurate and efficient geological surveys. First, the combined operation of rotary-wing drones and ground sensors ensures comprehensive coverage of both surface and subsurface data, making data collection more precise. Second, the fusion of multimodal data (such as lidar point clouds, seismic wave data, ground-penetrating radar data, and multispectral imagery) enhances information accuracy, enabling a comprehensive analysis of underground structures from multiple dimensions. Dynamically updated geological probabilistic models and risk assessment results enable this system to monitor potential geological risks in real time and issue timely warnings, providing important decision-making support for related fields (such as civil engineering and mineral exploration), improving exploration efficiency, and reducing disaster risks. This approach helps reduce human error, enhances surveying accuracy, and offers strong real-time and adaptability, making it a beneficial innovation in modern surveying and mapping.

[0101] In one possible implementation, prior to 3D scanning, a multispectral camera acquires image data of the target area and then calculates a vegetation index (such as NDVI). The vegetation index is a key parameter for measuring vegetation growth, reflecting its health and coverage. During this step, the vegetation index dynamically categorizes the scan area into low, medium, and high vegetation coverage zones. This classification allows for differentiated management of scanning requirements across different areas, ensuring optimized and efficient flight paths.

[0102] Based on the calculated vegetation index, a corresponding flight path is planned. Areas with a higher vegetation index (i.e., areas with high vegetation coverage) are scanned using a dense spiral flight path. This spiral flight path provides more detailed ground coverage, ensuring high 3D scanning density in high-coverage areas, thereby improving data quality and accuracy. During this process, the flight altitude is dynamically adjusted based on the vegetation index. When the vegetation index exceeds the first threshold, the flight altitude is adjusted to a preset minimum safe altitude. This adjustment is intended to minimize LiDAR signal attenuation caused by dense vegetation and to improve data acquisition accuracy.

[0103] During flight path planning, when the flight altitude is adjusted to the minimum safe altitude, the sampling density of the GPR mobile node must be increased simultaneously. Specifically, the sampling spacing of the GPR is inversely proportional to the flight altitude. That is, as the flight altitude decreases, the sampling spacing of the GPR decreases to maintain the fineness and high-density acquisition of radar data, thereby improving the spatial resolution and accuracy of the subsurface data.

[0104] This strategy of path planning and sampling density adjustment is an effective optimization of traditional exploration methods. It not only improves flight efficiency and data accuracy, but also provides a more solid data foundation for subsequent multimodal data fusion and risk assessment.

[0105] In one possible implementation, when collecting seismic wave propagation velocity data, it is first necessary to lay out a triangular grid sensor array. This grid layout can provide uniform coverage of the entire measurement area and improve the spatial resolution of the data. The setting of the node spacing is graded according to the complexity of the terrain. In flat areas, due to the small terrain undulations, the node spacing is set to 20 meters; in hilly areas, due to the high complexity of the terrain, higher precision is required, so the node spacing is set to 10 meters. This spacing setting can ensure the accurate collection of seismic wave data under different terrain conditions, while avoiding the waste of resources caused by oversampling.

[0106] After the sensor array is deployed, an artificial seismic source excites elastic waves (such as longitudinal or shear waves). These waves propagate through the subsurface and are reflected and refracted at different levels. The propagation speed of these waves is affected by the properties and structure of the subsurface medium. By placing multiple sensors in the sensor array, the propagation of seismic waves from the source to each sensor can be recorded, specifically the arrival time difference between the direct wave and the reflected wave. This arrival time difference data can be used to infer the propagation speed of the seismic waves.

[0107] After acquiring source-induced seismic wave propagation data, a velocity tomography algorithm is used to invert the velocity distribution of the subsurface medium. This inversion algorithm simulates the seismic wave propagation path and its reflection characteristics, combining it with actual arrival time difference data to calculate the wave velocities at different levels of the subsurface medium. During the inversion process, the algorithm optimizes the inversion results based on pre-defined constraints to improve the accuracy of the velocity distribution.

[0108] In the process of velocity tomography inversion, specific constraints need to be introduced to ensure the rationality and accuracy of the results. Specific constraints include:

[0109] Medium density range mapped by GPR dielectric constant data: The measurement results of GPR dielectric constant can provide a preliminary estimate of the density of the underground medium. Density is closely related to wave velocity, and this factor needs to be considered in the inversion process to ensure the physical rationality of wave velocity inversion.

[0110] LiDAR terrain curvature data reveals shallow rock mass integrity: LiDAR terrain data helps reveal the integrity of shallow rock masses. Changes in terrain curvature often correlate with changes in subsurface structure. This information provides a link between surface morphology and subsurface structure for velocity inversion, further improving the reliability of inversion results.

[0111] In one possible implementation, the first step in filtering out vegetation noise is to calculate the reflection intensity distribution histogram of the laser point cloud. The point cloud data collected by the laser scanner contains reflection signals from different objects, typically including reflections from multiple objects such as the ground and plants. Therefore, using the distribution histogram of reflection intensity can help us distinguish the reflection characteristics of different objects. Generally speaking, there is a difference in the reflection intensity of surface points and vegetation points. The reflection intensity of ground point clouds is usually higher, while the reflection intensity of vegetation points is weaker or more irregular.

[0112] Based on the reflection intensity distribution histogram, a dynamic threshold can be set to distinguish vegetation from surface point clouds. This threshold is not fixed but is dynamically adjusted based on atmospheric transmittance during the scanning period. Atmospheric transmittance affects the propagation of laser signals, which in turn affects the intensity of the reflected signal. To accurately distinguish vegetation from surface point clouds, the threshold must be adjusted in real time to ensure accurate filtering.

[0113] Dynamic correction of atmospheric transmittance is achieved by inverting the water vapor absorption band of multispectral imagery. This band provides information about atmospheric humidity and water vapor content, factors that directly affect laser signal propagation. By inverting atmospheric transmittance, the laser reflection intensity threshold can be precisely adjusted, thereby improving the accuracy of distinguishing vegetation from surface point clouds.

[0114] After initially screening the surface point cloud, the next step is to perform local surface fitting on these surface point clouds. The goal of local surface fitting is to describe the continuity of the surface through a mathematical model and remove outliers that do not conform to surface patterns. Specifically, the deviation between each point and the fitted surface is calculated. If the deviation exceeds a set adaptive threshold, the point is considered an outlier and removed.

[0115] The adaptive threshold is set by calculating terrain roughness. Terrain roughness reflects the undulation of the surface and is typically determined by the standard deviation of the point cloud elevation and the slope of the neighborhood. A larger elevation standard deviation indicates greater terrain variation; the slope of the neighborhood directly influences the distribution of ground points. Based on these terrain parameters, the adaptive threshold is dynamically adjusted to accurately remove outliers under varying terrain conditions.

[0116] In one possible implementation, constructing a joint physical property distribution model first requires establishing a joint objective function composed of seismic traveltime residuals and ground-penetrating radar reflection amplitudes. Seismic traveltime residuals are calculated by calculating the difference between the seismic wave propagation time and a predetermined seismic model, while ground-penetrating radar reflection amplitudes represent the intensity of the radar signal reflected after propagating through the subsurface medium. These two variables reflect the characteristics of different subsurface materials, and by combining them, a more comprehensive description of the physical properties of the subsurface structure can be achieved.

[0117] When establishing the joint objective function, a regional geological map is also required as a constraint to constrain the lithologic classification of each grid point. Lithologic classification represents the type of subsurface material, such as sandstone, limestone, and clay. These lithologic types appear differently in seismic and ground-penetrating radar signals. By converting lithologic classification into prior distributions of physical property parameters, more geological context can be incorporated into the joint model, thereby improving its accuracy.

[0118] Lithology categories are converted into probability distributions of physical parameters (such as density and elastic modulus) using prior knowledge. These physical parameters serve as constraints in the joint objective function, helping guide the particle swarm optimization algorithm's solution process and ensuring that the final optimization results conform to geological characteristics and physical laws.

[0119] To optimize the joint objective function, a modified particle swarm optimization (PSO) algorithm was employed. PSO, an optimization method that mimics the behavior of natural particle swarms, minimizes or maximizes the objective function by searching for the optimal solution in the solution space. This process specifically incorporates dynamic adjustment of the weighting coefficients for seismic and ground-penetrating radar data. This optimization optimizes the model by varying the contribution of these data to the objective function. The weighting coefficients are dynamically calculated based on the following factors:

[0120] Seismic wave signal-to-noise ratio: The signal-to-noise ratio affects the quality of seismic wave data. Data with a high signal-to-noise ratio can provide more accurate information. Therefore, the weight coefficient will increase the proportion of seismic data under high signal-to-noise ratio conditions.

[0121] GPR sampling density: Sampling density determines the resolution of GPR data. High-density GPR data usually provides more details, so the weight coefficient of the data will be increased in areas with high density.

[0122] The influence of surface temperature on electromagnetic wave propagation speed: The propagation speed of electromagnetic waves is affected by temperature. In ground-penetrating radar (GPR) applications, changes in surface temperature can affect the propagation characteristics of radar signals. Therefore, the model dynamically adjusts the weight of GPR data based on the current surface temperature.

[0123] In one possible implementation, when generating a probabilistic model of a 3D geological volume, various types of geological and remote sensing data, such as seismic velocity, dielectric constant, terrain curvature, and spectral anomaly index, are first mapped onto a unified spatial grid. These data represent the physical properties of the subsurface medium. Seismic velocity and dielectric constant reflect the density and dielectric properties of the subsurface material, terrain curvature relates to surface morphology, and spectral anomaly index reflects the specific characteristics of the subsurface minerals or structures. By mapping this information onto the same spatial grid, disparate data sources can be combined within a unified coordinate system, providing a unified foundation for subsequent modeling and analysis.

[0124] After obtaining the physical parameters on the unified grid, the next step is to calculate the void probability value corresponding to each grid point using nonparametric Bayesian methods. Nonparametric Bayesian methods are powerful statistical inference tools that automatically adjust the model complexity based on the data distribution. In this method, the void probability value reflects the presence of a void or cavity at a certain location (such as a karst area). The generation of voids is closely related to physical properties, and Bayesian methods gradually optimize the probability prediction by updating the prior distribution.

[0125] To accurately calculate the probability of a gap, it is necessary to construct a reasonable prior distribution. The key to this step is to build a model using the following data sources:

[0126] Karst is an important cause of the formation of underground voids. The historical disaster database records the development patterns of karst in different regions and can provide important prior information for the probability prediction of voids.

[0127] Changes in groundwater levels have a direct impact on the formation and evolution of voids. Real-time rainfall data can help estimate changes in groundwater levels, which is crucial for predicting the filling or dissolution of voids by groundwater.

[0128] During the calculation process, Markov Chain Monte Carlo (MCMC) sampling was used to optimize the probability distribution. MCMC is a statistical method that uses iterative sampling to generate samples that conform to the target distribution, thereby optimizing the probability distribution surface. Here, MCMC was used to optimize the distribution of void probability values on the grid, enabling the generated 3D geological model to more accurately reflect the actual underground void distribution.

[0129] In one possible implementation, during model operation, if the deviation between newly added seismic velocity data and the existing model's prediction exceeds a predefined historical fluctuation range—that is, the deviation falls outside the normal error range—the system automatically triggers a model confidence update. This indicates that the current model may not accurately reflect recent geological changes or that unaccounted-for factors are affecting the prediction results.

[0130] If the deviation exceeds a preset range, the system further analyzes the correlation between the deviation and relevant meteorological data, such as rainfall and soil moisture. This step aims to identify potential factors, particularly how meteorological conditions interact with subsurface geological properties, such as seismic wave velocity. By calculating the correlation coefficient between these data, it is possible to identify which meteorological factors significantly influence the model's forecast accuracy.

[0131] Based on the calculated correlation coefficients, the system automatically adjusts the weights of various data sources (such as seismic velocity, dielectric constant, and rainfall) in the joint distribution model of physical properties. This adjustment mechanism dynamically optimizes the contributions of different data sources, ensuring that the model more accurately reflects the actual underground environment in the face of new meteorological changes. Specifically, if the correlation between rainfall or soil moisture and seismic velocity is high, the weight of the relevant data will be increased, and vice versa.

[0132] After adjusting the data weights in the model, the system generates a new probability distribution of geological body boundaries, reflecting the model's revised predictions of underground structure boundaries. This revised value more accurately represents the actual distribution of geological bodies and improves the accuracy of predictions.

[0133] To further assess the reliability of the model, the system also generates a heat map of the spatial distribution of model uncertainty. This heat map shows the degree of uncertainty in the model's predictions across different regions, helping decision makers identify areas with more reliable predictions and those that may require more data or further analysis.

[0134] In one possible implementation, in the preliminary stage of the method, the three-dimensional geometric parameters of potential cavities are extracted from remote sensing data by using a probabilistic model. This mainly includes:

[0135] Roof Thickness: The thickness of the roof is a key factor in cavity stability. Thin roofs may be more susceptible to collapse and therefore need to be accurately assessed in the model.

[0136] Span ratio: The span ratio is the ratio of the horizontal span of a cavity to its vertical height, reflecting the morphological characteristics of the cavity. A larger span ratio indicates a cavity is more likely to become unstable, increasing the risk of collapse.

[0137] Surrounding rock crack development index: This parameter describes the degree of development of cracks in the surrounding rock. The number and distribution of cracks will significantly affect the mechanical properties of the rock mass, thereby affecting the stability of the cavity.

[0138] Based on the extraction of 3D geometric parameters, a comprehensive rainfall infiltration-groundwater level-mechanical stability coupling model was developed. This model simulates the rainfall process, groundwater level changes, and geotechnical behavior to calculate the critical conditions for potential collapse.

[0139] Rainfall infiltration: Rainfall intensity directly affects the infiltration of water into the soil, which in turn affects changes in the groundwater level. Real-time rainfall intensity information is fed into the model to predict the impact of precipitation on soil and rock.

[0140] Groundwater level: As rainfall infiltrates, the groundwater level may rise, increasing the moisture content of soil and rock, changing their mechanical properties, especially in hollow areas, which may cause further slip or collapse of unstable rock masses.

[0141] Mechanical stability: A coupled model considers the tensile strength of the rock mass and the interaction of groundwater with the rock mass to calculate the mechanical stability of the cavity. This stability can be affected by precipitation and rising groundwater levels.

[0142] The permeability coefficient is a key parameter in the model, directly affecting the rate of water infiltration and groundwater level fluctuations. To obtain this parameter, the dielectric constant of the soil layer is measured using ground penetrating radar (GPR) technology. The permeability coefficient is then inverted by mapping the dielectric constant to soil moisture content. Key to this process is the ability to non-invasively measure soil moisture content and calculate the permeability coefficient by linking it to soil permeability.

[0143] Through the calculation of the coupled model, the following parameters are finally output:

[0144] Safety Factor Over Time Curve: This curve reflects the evolution of cavity stability over time. With changes in rainfall and groundwater levels, the safety factor may decrease. When the safety factor falls below a critical value, a collapse risk may occur.

[0145] Warning level thresholds: Based on the calculated safety factor, the system automatically divides warning levels. For example, when the safety factor falls below a certain threshold, the system will issue a red warning, indicating a possible collapse and recommending emergency measures.

[0146] Through precise parameter extraction, dynamic model calculation and real-time data input, the prediction accuracy of potential collapse risks can be greatly improved, and strong support can be provided for emergency response.

[0147] Accordingly, an embodiment of the present invention further provides a remote sensing-based surveying and mapping instrument for executing a remote sensing-based surveying and mapping method described in an embodiment of the present invention, comprising:

[0148] Data acquisition module, data processing module, dynamic modeling module and risk assessment module, each module is connected through a high-speed data bus;

[0149] The data acquisition module includes:

[0150] The rotary-wing UAV unit integrates a penetrating laser radar, a multispectral imaging unit, and a navigation control unit. The laser radar's transmitting wavelength is configured to penetrate vegetation, and the multispectral imaging unit covers the visible to short-wave infrared band.

[0151] A ground sensing unit, comprising a seismic sensor array deployed in a triangular grid topology and a ground-penetrating radar mobile node equipped with an autonomous navigation chassis;

[0152] The meteorological monitoring unit, including a distributed rain gauge, soil moisture probe, and infrared thermometer, communicates with the main control unit via the LoRa wireless protocol;

[0153] The data processing module includes:

[0154] Point cloud filtering unit, connected to the laser radar data output terminal, with built-in reflection intensity analysis subunit and elevation surface fitting subunit;

[0155] A joint inversion unit receives output data from the seismic wave sensor array and the ground penetrating radar node and is equipped with a parallel computing accelerator;

[0156] The spectral correlation unit is connected to the data interface of the multispectral imaging unit and has a built-in band ratio calculator and anomaly feature extractor;

[0157] The dynamic modeling module includes:

[0158] The three-dimensional gridding unit receives the physical parameter matrix output by the data processing module and generates a voxel grid in a unified spatial coordinate system;

[0159] The probability calculation unit integrates the Bayesian update engine and Monte Carlo sampler, and is connected to the real-time data stream of the meteorological monitoring unit;

[0160] The risk assessment module includes:

[0161] The heat map generation unit receives the output data of the probability model and interacts with the geographic information system database;

[0162] Mechanical coupling analysis unit, with built-in finite element solver, and the output end connected to the early warning signal transmitter.

[0163] Preferably, in the data acquisition module, the rotor UAV unit and the ground sensor unit establish a collaborative control link, specifically including:

[0164] The UAV navigation control unit receives the position data of the ground penetrating radar mobile node in real time and dynamically adjusts the flight path to maintain complementary spatial coverage between the two;

[0165] The scanning frequency of the lidar and the sampling timing of the ground penetrating radar are aligned through hardware synchronization signals, with a time error of ≤1 millisecond;

[0166] The joint inversion unit of the data processing module is provided with a weight dynamic allocator, the input end of which is connected to:

[0167] The output of the seismic wave signal-to-noise ratio detector, where the signal-to-noise ratio detection is based on the power spectral density ratio of the direct wave to the background noise;

[0168] Ground penetrating radar sampling density counter, which counts the number of effective detection points within a unit area;

[0169] Surface temperature compensator, which corrects the electromagnetic wave propagation velocity model based on infrared thermometer data;

[0170] A spatial interpolation acceleration card is configured in the three-dimensional gridding unit of the dynamic modeling module, and its interpolation algorithm selection logic includes:

[0171] When the meteorological monitoring unit detects that the rainfall is greater than the threshold, it automatically switches to the Kriging interpolation mode with rainfall drift term;

[0172] In vegetation-covered areas, an anisotropic interpolation kernel function based on lidar terrain curvature is enabled;

[0173] A two-way feedback channel is provided between the mechanical coupling analysis unit and the dynamic modeling module of the risk assessment module:

[0174] The forward channel transmits the geometric parameters of the geological body to the finite element solver;

[0175] The reverse channel transmits the mechanical stability calculation results back to the probability calculation unit, triggering the model confidence correction.

[0176] The present invention encompasses any alternatives, modifications, equivalents, and solutions that fall within the spirit and scope of the present invention. To provide a thorough understanding of the present invention, specific details are described in detail below in connection with the preferred embodiments of the present invention, but those skilled in the art will be able to fully understand the present invention without these detailed descriptions. Furthermore, to avoid unnecessary confusion regarding the essence of the present invention, well-known methods, processes, procedures, components, and circuits have not been described in detail.

[0177] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A survey and mapping method based on remote sensing, characterized in that: The following steps are involved: Step 1: Dynamic data collection through multi-platform collaboration: Using a rotary-wing drone equipped with a penetrating lidar and a multispectral camera, a 3D scan of the vegetation-covered area is performed, simultaneously acquiring vegetation penetration point cloud data and multispectral image data. Deploy seismic wave sensor arrays and mobile ground-penetrating radar nodes in the target area to collect seismic wave propagation velocity data and dielectric constant distribution data of the underground medium, respectively; Real-time collection of rainfall, soil moisture and surface temperature data in the target area; Step 2: Multimodal data fusion and feature enhancement: filtering vegetation noise from the lidar point cloud data to reconstruct a surface terrain model under vegetation cover; Combine seismic wave velocity data with ground-penetrating radar dielectric constant data to construct a joint distribution model of underground medium physical properties; Correlate multispectral image reflectance characteristics with underground physical parameters to identify geological anomalies; Step 3: Dynamic modeling and risk assessment of hidden geological bodies: Generate a 3D geological body probability model based on the fused multi-source data; Dynamically integrate new sensor data and meteorological factors to update model confidence; Output multi-scale risk assessment results, including regional risk heat maps and critical collapse warning parameters.

2. The remote sensing-based surveying and mapping method according to claim 1, characterized in that: The flight path planning for the 3D scan described in step 1 includes: Calculate vegetation index through multispectral imagery and dynamically divide the scanned area into grades based on vegetation index; For areas with high vegetation coverage, a spiral intensive route is adopted, and the flight altitude is adjusted according to the calculated results of the vegetation index: When the vegetation index is higher than a first threshold, the flight altitude is lowered to a preset minimum safe altitude; The sampling density of the ground penetrating radar mobile node is increased synchronously, and the sampling interval is inversely proportional to the flight altitude.

3. The remote sensing-based surveying and mapping method according to claim 2, characterized in that: The acquisition of seismic wave propagation velocity data in step 1 includes: A triangular grid sensor array is deployed, and the node spacing is set according to the complexity of the terrain: The spacing is 20 meters in flat areas and 10 meters in hilly areas; The elastic waves are excited by artificial seismic sources, and the arrival time difference between the direct wave and the reflected wave is recorded; The velocity distribution of the underground medium is inverted based on the velocity tomography algorithm. The constraints in the velocity inversion process include: Dielectric density range mapped by GPR dielectric constant data; Indicators of shallow rock mass integrity as reflected by lidar terrain curvature data.

4. The remote sensing-based surveying and mapping method according to claim 3, characterized in that: The vegetation noise filtering in step 2 includes: Calculate the reflection intensity distribution histogram of the laser point cloud and set a dynamic threshold to distinguish vegetation from surface point clouds: The reflection intensity threshold is dynamically modified according to the atmospheric transmittance during the scanning period, and the atmospheric transmittance is obtained by inverting the water vapor absorption band of the multispectral image; Perform local surface fitting on the initially screened surface point cloud and remove abnormal points that deviate from the fitted surface by more than an adaptive threshold. The adaptive threshold is calculated based on the terrain roughness: The roughness index is determined by combining the standard deviation of the point cloud elevation and the neighborhood slope.

5. The remote sensing-based surveying and mapping method according to claim 4, characterized in that: The construction of the physical property joint distribution model in step 2 includes: Establish a joint objective function of seismic wave travel time residual and ground penetrating radar reflection amplitude; The regional geological map is introduced to constrain the lithologic categories of each grid point, and the lithologic categories are converted into prior distributions of physical property parameters; An improved particle swarm optimization algorithm is used to solve the objective function. During the optimization process, the weight coefficients of seismic data and ground penetrating radar data are dynamically adjusted. The weight coefficients are calculated based on the following factors: Seismic wave signal-to-noise ratio; GPR sampling density; The influence factor of surface temperature on the propagation speed of electromagnetic waves.

6. The remote sensing-based surveying and mapping method according to claim 5, characterized in that: The generation of the three-dimensional geological body probability model in step 3 includes: Mapping seismic wave velocity, dielectric constant, terrain curvature and spectral anomaly index to a unified spatial grid; The non-parametric Bayesian method is used to calculate the probability value of each grid point being a gap, where the prior distribution is constructed based on the following data: Karst development patterns in historical disaster databases; Groundwater level gradient estimated from real-time rainfall data; Optimize probability distribution surfaces via Markov chain Monte Carlo sampling.

7. The remote sensing-based surveying and mapping method according to claim 6, characterized in that: The model confidence update in step 3 includes: When the deviation between the newly added seismic velocity data and the model prediction exceeds the historical fluctuation range, the following update mechanism is triggered: Calculate the correlation coefficient between the deviation value and meteorological data, including rainfall and soil moisture content; Adjust the weight of each data source in the physical property joint distribution model according to the correlation coefficient; The updated model outputs include: Correction value of probability distribution of geological body boundary; Heatmap of the spatial distribution of model uncertainty.

8. The remote sensing-based surveying and mapping method according to claim 7, characterized in that: The generation of critical collapse warning parameters in step 3 includes: Extracting three-dimensional geometric parameters of potential cavities based on a probabilistic model, including roof thickness, span ratio, and surrounding rock fracture development index; A coupled model of rainfall infiltration, groundwater level, and mechanical stability was established to calculate the critical collapse conditions: Input parameters include real-time rainfall intensity, soil permeability coefficient and rock mass tensile strength; The permeability coefficient is obtained by inverting the mapping relationship between the dielectric constant of ground penetrating radar and soil moisture content; The output parameters include the curve of safety factor changing over time and the threshold value of warning level division.

9. A surveying and mapping instrument based on remote sensing, used for running a surveying and mapping method based on remote sensing according to any one of claims 1 to 8, characterized in that: include: Data acquisition module, data processing module, dynamic modeling module and risk assessment module, each module is connected through a high-speed data bus; The data acquisition module includes: The rotary-wing UAV unit integrates a penetrating laser radar, a multispectral imaging unit, and a navigation control unit. The laser radar's transmitting wavelength is configured to penetrate vegetation, and the multispectral imaging unit covers the visible to short-wave infrared band. A ground sensing unit, comprising a seismic sensor array deployed in a triangular grid topology and a ground-penetrating radar mobile node equipped with an autonomous navigation chassis; The meteorological monitoring unit, including a distributed rain gauge, soil moisture probe, and infrared thermometer, communicates with the main control unit via the LoRa wireless protocol; The data processing module includes: Point cloud filtering unit, connected to the laser radar data output terminal, with built-in reflection intensity analysis subunit and elevation surface fitting subunit; A joint inversion unit receives output data from the seismic wave sensor array and the ground penetrating radar node and is equipped with a parallel computing accelerator; The spectral correlation unit is connected to the data interface of the multispectral imaging unit and has a built-in band ratio calculator and anomaly feature extractor; The dynamic modeling module includes: The three-dimensional gridding unit receives the physical parameter matrix output by the data processing module and generates a voxel grid in a unified spatial coordinate system; The probability calculation unit integrates the Bayesian update engine and Monte Carlo sampler, and is connected to the real-time data stream of the meteorological monitoring unit; The risk assessment module includes: The heat map generation unit receives the output data of the probability model and interacts with the geographic information system database; Mechanical coupling analysis unit, with built-in finite element solver, and the output end connected to the early warning signal transmitter.

10. The remote sensing-based surveying and mapping instrument according to claim 9, characterized in that: In the data acquisition module, the rotor UAV unit and the ground sensor unit establish a collaborative control link, specifically including: The UAV navigation control unit receives the position data of the ground penetrating radar mobile node in real time and dynamically adjusts the flight path to maintain complementary spatial coverage between the two; The scanning frequency of the lidar and the sampling timing of the ground penetrating radar are aligned through hardware synchronization signals, with a time error of ≤1 millisecond; A weight dynamic allocator is provided in the joint inversion unit of the data processing module, the input end of which is connected to the output end of the seismic wave signal-to-noise ratio detector. The signal-to-noise ratio detection is based on the power spectrum density ratio of the direct wave and the background noise.

Images