Land engineering terrain rapid surveying method based on unmanned aerial vehicle

Abstract

The invention relates to the technical field of unmanned aerial vehicle surveying, in particular to a land engineering terrain rapid surveying method based on an unmanned aerial vehicle, which comprises the following steps: S1, acquiring the boundary of a target area, performing grid division, planning a cruising path of the unmanned aerial vehicle, controlling the unmanned aerial vehicle to fly to each grid point in sequence, and at each grid point, determining the cruising path of the unmanned aerial vehicle; acquiring a depth feedback signal generated by excitation of the first sound wave signal, and acquiring an earth surface vibration response signal generated by excitation of the second sound wave signal; s2, determining a continuous and stable underground reference interface in the target area based on the depth feedback signals collected by all the grid points; according to the invention, full-automatic gridding cruising and data acquisition and processing are realized based on the unmanned aerial vehicle, rapid general survey of a large-area area is realized, and the comprehensive survey cost and period are reduced.

Description

Technical Field

[0001] This invention relates to the field of unmanned aerial vehicle (UAV) surveying technology, and specifically to a rapid land engineering terrain surveying method based on UAVs. Background Technology

[0002] In recent years, with the development of drone technology, surveying methods based on drone platforms equipped with optical cameras, multispectral sensors or lidar have become the mainstream means of obtaining surface morphology data. However, these technologies can only achieve high-precision measurement of the surface geometry of landforms, vegetation or structures, and belong to the category of surface sensing. They cannot penetrate the surface to obtain information on shallow underground geological structures that are crucial for engineering construction, nor can they directly detect the engineering mechanical properties of surface rock and soil.

[0003] Traditional land engineering topographic surveys mainly rely on manual field investigations, drilling and sampling, and large-scale physical exploration equipment. Manual investigations are inefficient, have limited coverage, and are subject to subjective errors. Drilling methods are expensive, destructive, and can only obtain point information, making it difficult to reflect continuous regional changes. Large-scale physical equipment is often cumbersome to deploy, has long operation cycles, and is poorly adaptable to the environment. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs).

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: A rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs), comprising the following steps: S1. Obtain the boundary of the target area, divide it into grids and plan the drone's cruise path, control the drone to fly to each grid point in sequence, and at each grid point, obtain the depth feedback signal generated by the first acoustic signal and the ground vibration response signal generated by the second acoustic signal. S2. Based on the depth feedback signals collected from all grid points, determine a continuous and stable underground benchmark interface within the target area; S3. For each grid point, based on the depth feedback signal and surface vibration response signal of that grid point, and with the underground benchmark interface as a reference benchmark, determine the relative geological structure information and relative engineering material parameters of that grid point. S4. Based on the spatial coordinates of all grid points, relative geological structure information, and relative engineering material parameters, generate a comprehensive land engineering terrain model.

[0006] In a preferred embodiment, in step S1, boundary information of the target area is obtained, the target area is divided into grids based on the boundary information to form several grid points, and the cruise path of the UAV is planned. The cruise path is used to control the UAV to fly to each grid point in sequence. For each grid point, the drone is controlled to hover above it and vertically emit a first frequency band sound wave with a first penetration characteristic toward the ground, and receive a first vibration signal sequence generated by the first frequency band sound wave and containing reflection information from different underground interfaces, as a depth feedback signal for that grid point. For each grid point, a second-band sound wave with a second penetration characteristic is emitted to the ground to excite the ground surface, and a second vibration signal characterized by the micro-vibration of the ground surface generated by the excitation of the second-band sound wave is received as the ground vibration response signal of that grid point. The first penetration characteristic is used to enable the first frequency band sound waves to penetrate the earth's surface and reach the underground interface, and the second penetration characteristic is used to enable the second frequency band sound waves to couple with the shallow material of the earth's surface to excite micro-vibrations on the earth's surface.

[0007] In a preferred embodiment, in S2, based on the depth feedback signals collected at all grid points, the first vibration signal sequence contained in the depth feedback signal of each grid point is time aligned. In each first vibration signal sequence, the peak point of the characteristic waveform is selected as the reference zero time of the sequence. Then, the reference zero times of all sequences are aligned on the time coordinate. In each time-aligned first vibration signal sequence, the following parallel sub-steps are performed: Detect all local amplitude extreme points that exceed the preset amplitude threshold, and record the time corresponding to each extreme point as the arrival time of the reflected wave; For each detected reflected wave event, the integral value of the absolute amplitude or the peak amplitude of the waveform within a fixed-width time window centered on its arrival time is calculated as the waveform amplitude of the reflected wave. For the waveform segment of the same reflected wave event, its dominant frequency or the proportion of its energy distribution in multiple preset frequency bands is calculated through spectrum analysis, which is used as the frequency distribution of the reflected wave. And from each aligned first vibration signal sequence, a first feature set containing the arrival time of the reflected wave, waveform amplitude, and frequency distribution is extracted; The arrival time, waveform amplitude, and frequency distribution of the same reflected wave event are bound together as a feature tuple. The set of feature tuples of all reflected wave events in the sequence constitutes the first feature set of the grid point. Based on the first feature set, the adjacent grid points of each grid point are determined according to the predetermined spatial adjacency relationship. For each pair of adjacent grid points, the two sets of feature tuples with the closest arrival time of the reflected wave are selected from their respective first feature sets for pairing. Calculate the normalized cross-correlation value of the paired tuples on the waveform amplitude, and slide the calculation within a time window centered on the time difference of arrival. Take the maximum value as the first correlation coefficient of the reflected wave pair, and record the time difference of arrival corresponding to the maximum correlation coefficient as the spatial time difference. Based on the first correlation coefficient and spatial time difference calculated for all adjacent point pairs, the reflected wave group that maintains a high first correlation coefficient and a gradual change in spatial time difference between adjacent points is tracked to form a spatially traceable in-phase axis of reflected waves. For each tracked in-phase axis of the reflected wave, spatial continuity and characteristic consistency tests are performed. The spatial continuity test involves counting the total number of grid points traversed by the phase axis of the reflected wave and calculating its percentage of the total number of grid points in the target area. If the percentage is greater than or equal to a first preset ratio, the spatial continuity requirement is deemed to be met. The feature consistency test involves checking the first correlation coefficient between each grid point traversed by the phase axis and at least one adjacent grid point on the phase axis. It requires that the first correlation coefficient of all checked points be higher than a first preset threshold. At the same time, the waveform amplitude and main frequency data of the phase axis on all grid points are obtained, and the standard deviation of the maximum and minimum values ​​in the entire region is calculated. If the standard deviation is less than a second preset threshold, the feature consistency requirement is met. The in-phase axes of reflected waves that pass the spatial continuity test and the characteristic consistency test are determined to be stable reflected wave groups. The arrival time of the reflected waves at each grid point of the stable reflected wave group is converted into depth information according to a preset wave velocity model. The specific steps of the preset wave velocity model are as follows: S211. The preset wave velocity model is a structured parameter list. Based on the order from the surface downwards, the list divides the underground space from the surface to the expected depth into several hypothetical horizontal layered medium units. Two key parameters are assigned to each medium unit, including the layer number of the unit and a constant wave velocity value that uniquely corresponds to the unit. S212. When performing depth conversion on the reflection event of a stable reflection wave group, input the two-way travel time T of the reflected wave at a specific grid point of the event, where T represents the total time it takes for the wave to travel vertically down the surface to the target reflection interface and back to the surface. S213. Starting from the first layer, perform time allocation and thickness calculation layer by layer. Based on the wave velocity value of the first layer, calculate the time required for the wave to propagate vertically to one unit thickness in the medium of that layer. Allocate a portion of the total time T as the time consumed by the wave to propagate back and forth in the first layer. Based on this allocated time and the wave velocity of the first layer, the thickness of the first layer can be derived. After completing the calculation of the first layer, deduct the time allocated to the first layer from the total time T. The remaining time is the total time consumed by the wave to propagate in deeper layers. S214. Apply the process of S213 to each layer in the wave velocity model. Each layer uses its own constant wave velocity value and the time allocated to the layer from the remaining travel time in the previous step to calculate the thickness of the layer. This iterative process continues until the round-trip propagation time allocated to the latest layer corresponds exactly to the vertical distance from the top of the layer to the target reflection interface. At this point, the target reflection interface is located inside the layer. S215. The calculated thicknesses of all traversed layers in the wave velocity model are accumulated. The sum is the total vertical depth from the surface to the target reflection interface. For each grid point covered by the stable reflection wave group, the logical process from S212 to S215 is repeated to generate a depth value for each point. All depth values ​​together constitute a spatially continuously distributed dataset. The three-dimensional surface defined by this dataset is determined as the underground reference interface.

[0008] In a preferred embodiment, in step S3, based on the depth feedback signal of the current grid point, at least one local reflection interface is parsed out. First, the first vibration signal sequence of the current grid point is read. Then, in the sequence, the reflected wave groups known to correspond to the underground reference interface are identified and excluded. Then, in the remaining signal sequence, the first vibration signal sequence of the current grid point after time pair processing is read, and the preset amplitude threshold parameter is called. Starting from the second data point in the sequence and traversing to the second-to-last data point, positions with amplitude values ​​greater than both the preceding and following data points are marked as candidate local peaks, and positions with amplitude values ​​less than both the preceding and following data points are marked as candidate local troughs. Then, the absolute amplitude value of each candidate point is calculated and compared with a preset amplitude threshold. Candidate points with amplitude values ​​greater than or equal to the preset amplitude threshold are determined as valid reflected wave events and their time indices are recorded. Adjacent valid events with time intervals less than a preset minimum time window are compared, and only events with larger absolute amplitude values ​​are retained. This results in the final list of valid reflected wave events, where each event corresponds to a local reflection interface, and the arrival time of the corresponding reflected wave is recorded. The spatial relationship between the local reflection interface and the underground reference interface is determined. Based on the preset wave velocity model, the arrival time of the reflected wave of the local reflection interface is converted into a preliminary depth estimate. The known depth value of the underground reference interface at the current grid point is used as the reference reference depth. The difference between the preliminary depth estimate of each local reflection interface and the reference reference depth is calculated to obtain the preliminary relative depth relationship. The spectrum of the surface vibration response signal at the current grid point is analyzed to check whether there are abnormal resonance peaks, absorption valleys or obvious spectral distortions in the frequency components corresponding to the preliminary depth estimates of each local reflection interface. This phenomenon is defined as a characteristic response. Based on the presence and morphology of the characteristic responses, the preliminary relative depth relationship is verified and corrected. If there are clear characteristic responses, the existence of the local reflection interface is confirmed, and its attitude is further constrained by the frequency characteristics of the characteristic responses. If there are no obvious characteristic responses, the preliminary identification result may be determined to be unreliable or the interface property is weak. Finally, the accurate depth difference and spatial orientation relationship of each verified local reflection interface relative to the underground reference interface are output as the relative geological structure information of the current grid point. To separate the reference reflection signal component and extract its features from the depth feedback signal, firstly, in the first vibration signal sequence of the current grid point, based on the known arrival time window of the underground reference interface reflection wave group, a signal segment within that time window is extracted. This signal segment is the reference reflection signal component. Subsequently, feature extraction is performed on the extracted signal segment: spectral analysis is performed to obtain its amplitude spectrum, and the frequency point with the strongest energy in the amplitude spectrum is determined as its reference dominant frequency to extract the reference spectral features. At the same time, the envelope attenuation curve of the signal segment in the time domain is calculated to measure the time required for its amplitude to attenuate to a certain preset proportion of the initial value or to calculate the slope of its attenuation curve to extract the reference energy attenuation features. The complete second vibration signal of the current grid point is directly read as input, and spectrum analysis is performed to obtain its amplitude spectrum. The frequency point with the strongest energy in the amplitude spectrum is determined as its surface main frequency to extract the surface spectrum characteristics. At the same time, the envelope decay curve of the second vibration signal in the time domain is calculated to measure the time required for its amplitude to decay to the same preset ratio as the initial value or to calculate the slope of its decay curve to extract the surface energy decay characteristics. Subtracting the reference frequency from the surface main frequency gives the main frequency offset, which is the offset of the surface spectral characteristics relative to the reference spectral characteristics. Dividing the time required for the surface signal amplitude to decay to a preset ratio by the decay time of the reference signal gives a ratio, which is the ratio of the surface energy decay characteristics to the reference energy decay characteristics. If the calculated main frequency offset is positive and greater than the first set margin, the relative hardness parameter of the current grid point surface is determined to be higher than the reference. If the offset is negative and its absolute value is greater than the first set margin, the relative hardness parameter is determined to be lower than the reference. If the absolute value of the offset is less than or equal to the first set margin, the relative hardness parameter is determined to be close to the reference. If the calculated attenuation ratio is less than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter of the current grid point is determined to be higher than the benchmark. If the ratio is greater than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter is determined to be lower than the benchmark. If the difference between the ratio and 1 is less than or equal to the second set margin, then the relative density parameter is determined to be close to the benchmark. The first set margin and the second set margin are thresholds preset to distinguish between measurement noise and the actual difference.

[0009] In a preferred embodiment, in step S4, the spatial coordinates of all grid points, the relative geological structure information corresponding to each grid point, and the relative engineering material parameters corresponding to each grid point are obtained. The relative geological structure information includes at least the depth of the underground reference interface at each grid point, and the depth and spatial morphological relationship of the local reflective interface relative to the underground reference interface. The relative engineering material parameters include at least relative hardness and relative density parameters. Based on the spatial coordinates of all grid points and their corresponding depths of the underground reference interface, a first continuous surface is generated through spatial interpolation. This surface represents the three-dimensional shape of the underground reference interface. Subsequently, for each local reflection interface, the depth of the underground reference interface at each grid point is superimposed with the relative depth relationship of the interface at each grid point to obtain the absolute depth of the local reflection interface at each grid point. Based on the absolute depth and spatial coordinates, a subsequent continuous surface representing the local reflection interface is generated through spatial interpolation. Finally, the first continuous surface and all the subsequent continuous surfaces are combined and visualized according to their vertical spatial relationship to form a three-dimensional geological structure layer. Using relative hardness and relative density parameters as attribute values ​​and the spatial coordinates of grid points as the positioning basis, continuous attribute fields corresponding to hardness and density distributions are generated by two-dimensional spatial interpolation. Then, according to the preset color mapping rules, the continuous attribute fields are rendered to generate the two-dimensional engineering attribute layer expressed in the form of a planar diagram. The two-dimensional engineering attribute layer is used as the base plane map, and the planar projection of the three-dimensional geological structure layer is superimposed and associated with the base plane map in a unified geographic coordinate system to form an integrated model that supports synchronous query and visualization of spatial location.

[0010] The beneficial effects of the present invention are: by actively emitting first and second sound waves with different physical characteristics and receiving the depth feedback signal and surface vibration response signal generated by them, the present invention can simultaneously acquire structural information reflecting the underground stratum interface and attribute information reflecting the mechanical properties of the surface medium in a single flight operation. By utilizing the penetrating characteristics of the first sound wave, it is possible to acquire reflected signals from different underground interfaces, thereby constructing a three-dimensional model of the underground geological structure. By analyzing the differences between the surface vibration response signal and the vibration characteristics of the stable reference layer, the relative hardness and relative density parameters of the surface can be directly estimated. These parameters are directly related to key engineering indicators such as the bearing capacity, deformation characteristics, and compaction quality of the soil and rock mass, providing an intuitive and quantitative basis for decision-making in earthwork engineering, foundation treatment, and selection of construction machinery. By leveraging drones to achieve fully automated grid-based patrols, data collection, and processing, rapid surveys of large areas have been realized, reducing overall survey costs and timelines. Attached Figure Description

[0011] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0012] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0013] In the description of this application, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of the stated features. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0014] In the description of this application, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this application is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed in this application.

[0015] like Figure 1 This embodiment provides a rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs), comprising the following steps: S1. Obtain the boundary of the target area, divide it into grids and plan the drone's cruise path, control the drone to fly to each grid point in sequence, and at each grid point, obtain the depth feedback signal generated by the first acoustic signal and the ground vibration response signal generated by the second acoustic signal. Furthermore, in S1, the boundary information of the target area is obtained, and the target area is divided into grids based on the boundary information to form several grid points. The cruise path of the UAV is planned, and the cruise path is used to control the UAV to fly to each grid point in sequence. For each grid point, the drone is controlled to hover above it and vertically emit a first frequency band sound wave with a first penetration characteristic toward the ground, and receive a first vibration signal sequence generated by the first frequency band sound wave and containing reflection information from different underground interfaces, as a depth feedback signal for that grid point. For each grid point, a second-band sound wave with a second penetration characteristic is emitted to the ground to excite the ground surface, and a second vibration signal characterized by the micro-vibration of the ground surface generated by the excitation of the second-band sound wave is received as the ground vibration response signal of that grid point. The first penetration characteristic is used to enable the first frequency band sound waves to penetrate the earth's surface and reach the underground interface, and the second penetration characteristic is used to enable the second frequency band sound waves to couple with the shallow material of the earth's surface to excite micro-vibrations on the earth's surface.

[0016] S2. Based on the depth feedback signals collected from all grid points, determine a continuous and stable underground benchmark interface within the target area; Furthermore, in S2, based on the depth feedback signals collected at all grid points, the first vibration signal sequence contained in the depth feedback signal of each grid point is time aligned. In each first vibration signal sequence, the peak point of the characteristic waveform is selected as the reference zero time of the sequence. Subsequently, the reference zero times of all sequences are aligned on the time coordinate. In each time-aligned first vibration signal sequence, the following parallel sub-steps are performed: Detect all local amplitude extreme points that exceed the preset amplitude threshold, and record the time corresponding to each extreme point as the arrival time of the reflected wave; The preset amplitude threshold is used to distinguish meaningful reflected wave events from background noise or small fluctuations from the vibration signal sequence. Only peaks or troughs in the signal with an absolute amplitude value exceeding this threshold are considered valid reflected events. This threshold is preset through statistical analysis of typical background noise levels in the target area or through prior experimental data.

[0017] For each detected reflected wave event, the integral value of the absolute amplitude or the peak amplitude of the waveform within a fixed-width time window centered on its arrival time is calculated as the waveform amplitude of the reflected wave. For the waveform segment of the same reflected wave event, its dominant frequency or the proportion of its energy distribution in multiple preset frequency bands is calculated through spectrum analysis, which is used as the frequency distribution of the reflected wave. The preset frequency band refers to the specific frequency range used to define energy distribution calculations when extracting frequency distribution characteristics. It includes a set of predefined frequency intervals (e.g., low frequency band, mid frequency band, high frequency band) that are pre-divided based on the emission spectrum of the acoustic equipment used in the survey and the expected frequency response characteristics of the surface and underground media.

[0018] And from each aligned first vibration signal sequence, a first feature set containing the arrival time of the reflected wave, waveform amplitude, and frequency distribution is extracted; The arrival time, waveform amplitude, and frequency distribution of the same reflected wave event are bound together as a feature tuple. The set of feature tuples of all reflected wave events in the sequence constitutes the first feature set of the grid point. Based on the first feature set, the adjacent grid points of each grid point are determined according to the predetermined spatial adjacency relationship. For each pair of adjacent grid points, the two sets of feature tuples with the closest arrival time of the reflected wave are selected from their respective first feature sets for pairing. Calculate the normalized cross-correlation value of the paired tuples on the waveform amplitude, and slide the calculation within a time window centered on the time difference of arrival. Take the maximum value as the first correlation coefficient of the reflected wave pair, and record the time difference of arrival corresponding to the maximum correlation coefficient as the spatial time difference. Based on the first correlation coefficient and spatial time difference calculated for all adjacent point pairs, the reflected wave group that maintains a high first correlation coefficient and a gradual change in spatial time difference between adjacent points is tracked to form a spatially traceable in-phase axis of reflected waves. For each tracked in-phase axis of the reflected wave, spatial continuity and characteristic consistency tests are performed. The spatial continuity test involves counting the total number of grid points traversed by the phase axis of the reflected wave and calculating its percentage of the total number of grid points in the target area. If the percentage is greater than or equal to a first preset ratio, the spatial continuity requirement is deemed to be met. The feature consistency test involves checking the first correlation coefficient between each grid point traversed by the phase axis and at least one adjacent grid point on the phase axis. It requires that the first correlation coefficient of all checked points be higher than a first preset threshold. At the same time, the waveform amplitude and main frequency data of the phase axis on all grid points are obtained, and the standard deviation of the maximum and minimum values ​​in the entire region is calculated. If the standard deviation is less than a second preset threshold, the feature consistency requirement is met. The first preset ratio is a key quantitative standard for determining whether the reflection in-phase axes have a regional distribution. This represents the minimum proportion of grid points that the reflecting phase axis must cover. It is preset according to the continuity of the regional geological structure and the requirements of survey accuracy. For example, if it is expected to identify stable interfaces throughout the region, this value is high, allowing for local variations. If this value is low, the first preset threshold is used in the characteristic consistency test of stable reflection wave groups to determine whether the similarity of adjacent point signals is high enough. It is a correlation coefficient threshold between 0 and 1. Only when the first correlation coefficient of adjacent point signals is higher than this value is the waveform characteristic considered consistent.

[0019] The in-phase axes of reflected waves that pass the spatial continuity test and the characteristic consistency test are determined to be stable reflected wave groups. The arrival time of the reflected waves at each grid point of the stable reflected wave group is converted into depth information according to a preset wave velocity model. The specific steps of the preset wave velocity model are as follows: S211. The preset wave velocity model is a structured parameter list. Based on the order from the surface downwards, the underground space from the surface to the expected depth is divided into several hypothetical horizontal layered medium units. Two key parameters are assigned to each medium unit, including the layer number of the unit (such as the first layer, the second layer), and a constant wave velocity value that uniquely corresponds to the unit. S212. When performing depth transformation on the reflection event of a stable reflection wave group, input the two-way travel time T of the reflected wave at a specific grid point of the event. T represents the total time taken for the wave to travel vertically down from the ground surface to the target reflection interface and back to the ground surface. According to the layered structure defined by the wave velocity model, the total time T is regarded as the sum of the round-trip time taken by the wave when traversing each layer of medium in the model. S213. Starting from the first layer, perform time allocation and thickness calculation layer by layer. Based on the wave velocity value of the first layer, calculate the time required for the wave to propagate vertically to one unit thickness in the medium of that layer. Allocate a portion of the total time T as the time consumed by the wave to propagate back and forth in the first layer. Based on this allocated time and the wave velocity of the first layer, the thickness of the first layer can be derived. After completing the calculation of the first layer, deduct the time allocated to the first layer from the total time T. The remaining time is the total time consumed by the wave to propagate in deeper layers. S214. Apply the process of S213 to each layer in the wave velocity model. Each layer uses its own constant wave velocity value and the time allocated to the layer from the remaining travel time in the previous step to calculate the thickness of the layer. This iterative process continues until the round-trip propagation time allocated to the latest layer corresponds exactly to the vertical distance from the top of the layer to the target reflection interface. At this point, the target reflection interface is located inside the layer. S215. The calculated thicknesses of all traversed layers in the wave velocity model (from the first layer to the layer where the target interface is located) are accumulated. The sum is the total vertical depth from the surface to the target reflection interface. For each grid point covered by the stable reflection wave group, the logical process from S212 to S215 is repeated to generate a depth value for each point. All depth values ​​together constitute a spatially continuously distributed dataset. The three-dimensional surface defined by this dataset is determined as the underground reference interface.

[0020] S3. For each grid point, based on the depth feedback signal and surface vibration response signal of that grid point, and with the underground benchmark interface as a reference benchmark, determine the relative geological structure information and relative engineering material parameters of that grid point. Furthermore, in S3, based on the depth feedback signal of the current grid point, at least one local reflection interface is parsed out. First, the first vibration signal sequence of the current grid point is read. Then, in the sequence, the reflected wave groups known to correspond to the underground reference interface are identified and excluded. Then, in the remaining signal sequence, the first vibration signal sequence of the current grid point after time pair processing is read, and the preset amplitude threshold parameter is called. Starting from the second data point in the sequence and traversing to the second-to-last data point, positions with amplitude values ​​greater than both the preceding and following data points are marked as candidate local peaks, and positions with amplitude values ​​less than both the preceding and following data points are marked as candidate local troughs. Then, the absolute amplitude value of each candidate point is calculated and compared with a preset amplitude threshold. Candidate points with amplitude values ​​greater than or equal to the preset amplitude threshold are determined as valid reflected wave events and their time indices are recorded. Adjacent valid events with time intervals less than a preset minimum time window are compared, and only events with larger absolute amplitude values ​​are retained. This results in the final list of valid reflected wave events, where each event corresponds to a local reflection interface, and the arrival time of the corresponding reflected wave is recorded. The preset minimum time window is used in the merging process after reflection event detection to determine whether two adjacent reflection events may originate from the same physical interface, in order to avoid redundant records caused by waveform oscillation. It is a minimum time interval. If the time difference between two valid reflection events is less than this window, they are considered to be possible duplicates and need to be merged. It is pre-calculated and set according to the dominant wavelength (or wavelet length) of seismic waves or sound waves in shallow media and the time resolution of the system.

[0021] The spatial relationship between the local reflection interface and the underground reference interface is determined. Based on the preset wave velocity model, the arrival time of the reflected wave of the local reflection interface is converted into a preliminary depth estimate. The known depth value of the underground reference interface at the current grid point is used as the reference reference depth. The difference between the preliminary depth estimate of each local reflection interface and the reference reference depth is calculated to obtain the preliminary relative depth relationship. The spectrum of the surface vibration response signal at the current grid point is analyzed to check whether there are abnormal resonance peaks, absorption valleys or obvious spectral distortions in the frequency components corresponding to the preliminary depth estimates of each local reflection interface. This phenomenon is defined as a characteristic response. Based on the presence and morphology of the characteristic responses, the preliminary relative depth relationship is verified and corrected. If there are clear characteristic responses, the existence of the local reflection interface is confirmed, and its attitude is further constrained by the frequency characteristics of the characteristic responses. If there are no obvious characteristic responses, the preliminary identification result may be determined to be unreliable or the interface property is weak. Finally, the accurate depth difference and spatial orientation relationship of each verified local reflection interface relative to the underground reference interface are output as the relative geological structure information of the current grid point. The occurrence is a geological term used to describe the spatial morphology of a local interface. It is part of the relative geological structure information and refers to the direction and inclination of the geological interface in three-dimensional space, expressed by strike, dip and dip angle. To separate the reference reflection signal component and extract its features from the depth feedback signal, firstly, in the first vibration signal sequence of the current grid point, based on the known arrival time window of the underground reference interface reflection wave group, a signal segment within that time window is extracted. This signal segment is the reference reflection signal component. Subsequently, feature extraction is performed on the extracted signal segment: spectral analysis is performed to obtain its amplitude spectrum, and the frequency point with the strongest energy in the amplitude spectrum is determined as its reference dominant frequency to extract the reference spectral features. At the same time, the envelope attenuation curve of the signal segment in the time domain is calculated to measure the time required for its amplitude to attenuate to a certain preset proportion of the initial value or to calculate the slope of its attenuation curve to extract the reference energy attenuation features. It should be noted that the specific logical steps for calculating the envelope attenuation curve of the signal segment in the time domain and measuring the attenuation time or calculating the slope include: taking the absolute value of the input signal segment to obtain an absolute value signal, then performing low-pass filtering on the absolute value signal to extract the envelope, wherein the cutoff frequency of the low-pass filter is preset according to the main frequency components and attenuation characteristics of the signal, then determining the initial amplitude value of the envelope, specifically by taking the average amplitude of a preset number of data points in the beginning part of the envelope, and if it is necessary to measure the time required for the amplitude to attenuate to a certain preset proportion of the initial value, then based on the preset proportion parameter, the target amplitude value is calculated to be equal to the initial amplitude value multiplied by the proportion parameter; Starting from the beginning of the envelope, compare each data point along the positive time axis and record the time corresponding to the first data point whose amplitude value is less than or equal to the target amplitude value. The time difference between this time and the beginning point is the time required to decay to the preset ratio. The slope of the decay curve is calculated based on another preset proportional parameter. The corresponding target amplitude value is calculated. Starting from the starting point of the envelope, the first data point with an amplitude value less than or equal to the target amplitude value is found along the positive time axis. This point is taken as the end point of the decay segment, and the starting point of the envelope is taken as the starting point of the decay segment. The average rate of change of the amplitude value of the envelope with time in the decay segment is calculated. That is, the amplitude value at the end point is subtracted from the amplitude value at the starting point, and then divided by the difference between the time at the end point and the time at the starting point. The resulting value is the slope of the decay curve.

[0022] The complete second vibration signal of the current grid point is directly read as input, and spectrum analysis is performed to obtain its amplitude spectrum. The frequency point with the strongest energy in the amplitude spectrum is determined as its surface main frequency to extract the surface spectrum characteristics. At the same time, the envelope decay curve of the second vibration signal in the time domain is calculated to measure the time required for its amplitude to decay to the same preset ratio as the initial value or to calculate the slope of its decay curve to extract the surface energy decay characteristics. The same preset ratio provides a unified and standardized measurement point when calculating signal attenuation time, making the attenuation rate of different signal segments comparable. The first preset margin is used to distinguish between the actual material hardness difference and the frequency shift caused by measurement noise or small changes when determining the relative hardness parameter. It is an allowable range threshold for the main frequency offset. If the absolute value of the main frequency offset is less than or equal to this margin, it is judged to be close to the benchmark. It is preset based on the statistical analysis of the repeatability of the main frequency measurement of the benchmark layer signal or the engineering requirement of not being sensitive to small changes in the hardness of the ground material. The second preset margin is used to distinguish between the actual density difference and the attenuation ratio fluctuation caused by measurement noise when determining the relative density parameter. If the absolute value of the difference between the attenuation ratio and 1 is less than or equal to this margin, it is judged to be close to the benchmark. It is preset based on the statistical analysis of the repeatability of the measurement of the attenuation characteristics of the benchmark layer signal or the engineering tolerance of not being sensitive to small changes in density.

[0023] Subtracting the reference frequency from the surface main frequency gives the main frequency offset, which is the offset of the surface spectral characteristics relative to the reference spectral characteristics. Dividing the time required for the surface signal amplitude to decay to a preset ratio by the decay time of the reference signal gives a ratio, which is the ratio of the surface energy decay characteristics to the reference energy decay characteristics. If the calculated main frequency offset is positive and greater than the first set margin, the relative hardness parameter of the current grid point surface is determined to be higher than the reference. If the offset is negative and its absolute value is greater than the first set margin, the relative hardness parameter is determined to be lower than the reference. If the absolute value of the offset is less than or equal to the first set margin, the relative hardness parameter is determined to be close to the reference. If the calculated attenuation ratio is less than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter of the current grid point is determined to be higher than the benchmark. If the ratio is greater than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter is determined to be lower than the benchmark. If the difference between the ratio and 1 is less than or equal to the second set margin, then the relative density parameter is determined to be close to the benchmark. The first set margin and the second set margin are thresholds preset to distinguish between measurement noise and the actual difference.

[0024] S4. Based on the spatial coordinates of all grid points, relative geological structure information, and relative engineering material parameters, generate a comprehensive land engineering terrain model.

[0025] Furthermore, in step S4, the spatial coordinates of all grid points, the relative geological structure information corresponding to each grid point, and the relative engineering material parameters corresponding to each grid point are obtained. The relative geological structure information includes at least the depth of the underground reference interface at each grid point, and the depth and spatial morphology relationship of the local reflection interface relative to the underground reference interface. The relative engineering material parameters include at least the relative hardness parameter and the relative density parameter. Based on the spatial coordinates of all grid points and their corresponding depths of the underground reference interface, a first continuous surface is generated through spatial interpolation. This surface represents the three-dimensional shape of the underground reference interface. Subsequently, for each local reflection interface, the depth of the underground reference interface at each grid point is superimposed with the relative depth relationship of the interface at each grid point to obtain the absolute depth of the local reflection interface at each grid point. Based on the absolute depth and spatial coordinates, a subsequent continuous surface representing the local reflection interface is generated through spatial interpolation. Finally, the first continuous surface and all the subsequent continuous surfaces are combined and visualized according to their vertical spatial relationship to form a three-dimensional geological structure layer. Using relative hardness and relative density parameters as attribute values ​​and the spatial coordinates of grid points as the positioning basis, continuous attribute fields corresponding to hardness and density distributions are generated by two-dimensional spatial interpolation. Then, according to the preset color mapping rules, the continuous attribute fields are rendered to generate the two-dimensional engineering attribute layer expressed in the form of a planar diagram. The two-dimensional engineering attribute layer is used as the base plane map, and the planar projection of the three-dimensional geological structure layer is superimposed and associated with the base plane map in a unified geographic coordinate system to form an integrated model that supports synchronous query and visualization of spatial location.

[0026] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0027] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0028] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0029] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0030] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0031] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0032] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A rapid land engineering topography survey method based on unmanned aerial vehicles (UAVs), characterized in that, Includes the following steps: S1. Obtain the boundary of the target area, divide it into grids, plan the drone's cruise path, and control the drone to fly to each grid point in sequence. At each grid point, obtain the depth feedback signal generated by the first acoustic signal and the ground vibration response signal generated by the second acoustic signal. S2. Based on the depth feedback signals collected from all grid points, determine a continuous and stable subsurface reference interface within the target area; S3. For each grid point, based on the depth feedback signal and surface vibration response signal of that grid point, and using the underground benchmark interface as a reference, determine the relative geological structure information and relative engineering material parameters of that grid point; S4. Based on the spatial coordinates of all grid points, relative geological structure information, and relative engineering material parameters, generate a comprehensive land engineering terrain model.

2. The rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, In step S1, the boundary information of the target area is obtained, and the target area is divided into grids based on the boundary information to form several grid points. The cruise path of the UAV is planned, and the cruise path is used to control the UAV to fly to each grid point in sequence. For each grid point, the UAV is controlled to hover above it and vertically emit a first-band sound wave with a first penetration characteristic toward the ground, and receive a first vibration signal sequence generated by the first-band sound wave and containing reflection information from different underground interfaces, as the depth feedback signal for that grid point; for each grid point, a second-band sound wave with a second penetration characteristic is emitted toward the ground to excite the surface, and a second vibration signal characterized by the micro-vibration features of the surface generated by the excitation of the second-band sound wave is received, as the surface vibration response signal for that grid point; The first penetration characteristic is used to enable the first frequency band sound waves to penetrate the earth's surface and reach the underground interface, and the second penetration characteristic is used to enable the second frequency band sound waves to couple with the shallow material of the earth's surface to excite micro-vibrations on the earth's surface.

3. The rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, In step S2, based on the depth feedback signals collected at all grid points, the first vibration signal sequence contained in the depth feedback signal of each grid point is time-aligned. In each first vibration signal sequence, the peak point of the characteristic waveform is selected as the reference zero time of the sequence, and then the reference zero times of all sequences are aligned on the time coordinate.

4. The rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 3, characterized in that, In each time-aligned first vibration signal sequence, the following parallel sub-steps are performed: Detect all local amplitude extremes that exceed the preset amplitude threshold, and record the time corresponding to each extreme point as the arrival time of the reflected wave. For each detected reflected wave event, the integral value of the absolute amplitude or the peak amplitude of the waveform within a fixed-width time window centered on its arrival time is calculated as the waveform amplitude of the reflected wave; for the waveform segment of the same reflected wave event, its dominant frequency or the proportion of its energy distribution in multiple preset frequency bands is calculated through spectrum analysis as the frequency distribution of the reflected wave. And from each aligned first vibration signal sequence, a first feature set containing the arrival time of the reflected wave, waveform amplitude, and frequency distribution is extracted; The arrival time, waveform amplitude, and frequency distribution of the same reflected wave event are bound together as a feature tuple. The set of feature tuples of all reflected wave events in the sequence constitutes the first feature set of the grid point. Based on the first feature set, the adjacent grid points of each grid point are determined according to the predetermined spatial adjacency relationship. For each pair of adjacent grid points, the two sets of feature tuples with the closest arrival time of the reflected wave are selected from their respective first feature sets for pairing. Calculate the normalized cross-correlation value of the paired tuples on the waveform amplitude, and slide the calculation within a time window centered on the time difference of arrival. Take the maximum value as the first correlation coefficient of the reflected wave pair, and record the time difference of arrival corresponding to the maximum correlation coefficient as the spatial time difference. Based on the first correlation coefficient and spatial time difference calculated for all adjacent point pairs, the reflected wave group that maintains a high first correlation coefficient and a gradual change in spatial time difference between adjacent points is tracked to form a spatially traceable in-phase axis of reflected waves. For each tracked phase axis of the reflected wave, a spatial continuity test and a feature consistency test are performed. The spatial continuity test is performed by counting the total number of grid points traversed by the phase axis of the reflected wave and calculating its percentage of the total number of grid points in the target area. If the percentage is greater than or equal to a first preset ratio, the spatial continuity requirement is determined to be met. Among them, the feature consistency test checks the first correlation coefficient between each grid point through which the phase axis passes and at least one adjacent grid point on the phase axis. It requires that the first correlation coefficient of all checked points is higher than the first preset threshold. At the same time, the waveform amplitude and main frequency data of the phase axis on all grid points are obtained, and the standard deviation of the maximum and minimum values ​​in the whole region is calculated. When the standard deviation is less than the second preset threshold, the feature consistency requirement is met. The in-phase axes of reflected waves that pass the spatial continuity test and the characteristic consistency test are determined to be stable reflected wave groups. The arrival time of the reflected waves at each grid point of the stable reflected wave group is converted into depth information according to the preset wave velocity model.

5. A rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 4, characterized in that, The specific steps of the preset wave velocity model are as follows: S211. The preset wave velocity model is a structured parameter list. Based on the order from the surface downwards, the list divides the underground space from the surface to the expected depth into several hypothetical horizontal layered medium units. Two key parameters are assigned to each medium unit, including the layer number of the unit and a constant wave velocity value that uniquely corresponds to the unit. S212. When performing depth conversion on the reflection event of a stable reflection wave group, input the two-way travel time T of the reflected wave at a specific grid point of the event, where T represents the total time it takes for the wave to travel vertically down the surface to the target reflection interface and back to the surface. S213. Starting from the first layer, perform time allocation and thickness calculation layer by layer. Based on the wave velocity value of the first layer, calculate the time required for the wave to propagate vertically to one unit thickness in the medium of that layer. Allocate a portion of the total time T as the time consumed by the wave to propagate back and forth in the first layer. Based on this allocated time and the wave velocity of the first layer, the thickness of the first layer can be derived. After completing the calculation of the first layer, deduct the time allocated to the first layer from the total time T. The remaining time is the total time consumed by the wave to propagate in deeper layers. S214. Apply the process of S213 to each layer in the wave velocity model. Each layer uses its own constant wave velocity value and the time allocated to the layer from the remaining travel time in the previous step to calculate the thickness of the layer. This iterative process continues until the round-trip propagation time allocated to the latest layer corresponds exactly to the vertical distance from the top of the layer to the target reflection interface. At this point, the target reflection interface is located inside the layer. S215. The calculated thicknesses of all traversed layers in the wave velocity model are accumulated. The sum is the total vertical depth from the surface to the target reflection interface. For each grid point covered by the stable reflection wave group, the logical process from S212 to S215 is repeated to generate a depth value for each point. All depth values ​​together constitute a spatially continuously distributed dataset. The three-dimensional surface defined by this dataset is determined as the underground reference interface.

6. The rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, Based on the depth feedback signal of the current grid point, at least one local reflection interface is analyzed. First, the first vibration signal sequence of the current grid point is read. Then, in this sequence, the reflected wave groups known to correspond to the underground reference interface are identified and excluded. Then, in the remaining signal sequence, the first vibration signal sequence of the current grid point after time pairing processing is read, and the preset amplitude threshold parameter is called. Starting from the second data point in the sequence and traversing to the second-to-last data point, positions with amplitude values ​​greater than both their immediate and adjacent data points are marked as candidate local peaks, and positions with amplitude values ​​less than their immediate and adjacent data points are marked as candidate local troughs. Next, the absolute amplitude value of each candidate point is calculated and compared with a preset amplitude threshold. Candidate points with amplitude values ​​greater than or equal to the preset threshold are determined as valid reflected wave events, and their time indices are recorded. Adjacent valid events with time intervals less than a preset minimum time window are compared, and only events with larger absolute amplitude values ​​are retained. This results in the final list of valid reflected wave events, where each event corresponds to a local reflection interface, and the arrival time of the corresponding reflected wave is recorded. The spatial relationship between the local reflecting interface and the underground reference interface is determined. Based on the preset wave velocity model, the arrival time of the reflected wave at the local reflecting interface is converted into a preliminary depth estimate. The known depth value of the underground reference interface at the current grid point is used as the reference depth. The difference between the preliminary depth estimate of each local reflecting interface and the reference depth is calculated to obtain the preliminary relative depth relationship. The spectrum of the surface vibration response signal at the current grid point is analyzed to check whether there are abnormal resonance peaks, absorption valleys or obvious spectral distortions in the frequency components corresponding to the preliminary depth estimates of each local reflecting interface. This phenomenon is defined as the characteristic response.

7. A rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 6, characterized in that, Based on the presence and morphology of the characteristic responses, the preliminary relative depth relationship is verified and corrected. If there are clear characteristic responses, the existence of the local reflection interface is confirmed, and its attitude is further constrained by the frequency characteristics of the characteristic responses. If there are no obvious characteristic responses, the preliminary identification result may be determined to be unreliable or the interface property is weak. Finally, the accurate depth difference and spatial orientation relationship of each verified local reflection interface relative to the underground reference interface are output as the relative geological structure information of the current grid point. To separate the reference reflection signal component and extract its features from the depth feedback signal, firstly, in the first vibration signal sequence of the current grid point, based on the known arrival time window of the underground reference interface reflection wave group, a signal segment within that time window is extracted. This signal segment is the reference reflection signal component. Subsequently, feature extraction is performed on the extracted signal segment: spectral analysis is performed to obtain its amplitude spectrum, and the frequency point with the strongest energy in the amplitude spectrum is determined as its reference dominant frequency to extract the reference spectral features. At the same time, the envelope attenuation curve of the signal segment in the time domain is calculated to measure the time required for its amplitude to decay to a certain preset proportion of the initial value or to calculate the slope of its attenuation curve to extract the reference energy attenuation features. The complete second vibration signal of the current grid point is directly read as input, and spectrum analysis is performed to obtain its amplitude spectrum. The frequency point with the strongest energy in the amplitude spectrum is determined as its surface main frequency to extract the surface spectrum characteristics. At the same time, the envelope decay curve of the second vibration signal in the time domain is calculated to measure the time required for its amplitude to decay to the same preset ratio as the initial value or to calculate the slope of its decay curve to extract the surface energy decay characteristics. Subtracting the reference frequency from the surface dominant frequency yields the dominant frequency offset, which is the offset of the surface spectral characteristics relative to the reference spectral characteristics. Dividing the time required for the surface signal amplitude to decay to a preset ratio by the decay time of the reference signal yields a ratio, which is the ratio of the surface energy decay characteristics to the reference energy decay characteristics.

8. A rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 7, characterized in that, If the calculated main frequency offset is positive and greater than the first set margin, the relative hardness parameter of the current grid point surface is determined to be higher than the reference. If the offset is negative and its absolute value is greater than the first set margin, the relative hardness parameter is determined to be lower than the reference. If the absolute value of the offset is less than or equal to the first set margin, the relative hardness parameter is determined to be close to the reference. If the calculated attenuation ratio is less than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter of the current grid point is determined to be higher than the benchmark. If the ratio is greater than 1 and the difference between it and 1 is greater than the second set margin, then the relative density parameter is determined to be lower than the benchmark. If the difference between the ratio and 1 is less than or equal to the second set margin, then the relative density parameter is determined to be close to the benchmark. The first set margin and the second set margin are thresholds preset to distinguish between measurement noise and actual differences.

9. A rapid land engineering terrain survey method based on unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, In step S4, the spatial coordinates of all grid points, the relative geological structure information corresponding to each grid point, and the relative engineering material parameters corresponding to each grid point are obtained. The relative geological structure information includes at least the depth of the underground reference interface at each grid point, and the depth relationship and spatial morphology relationship of the local reflection interface relative to the underground reference interface. The relative engineering material parameters include at least the relative hardness parameter and the relative density parameter. Based on the spatial coordinates of all grid points and their corresponding depths of the underground reference interface, a first continuous surface is generated through spatial interpolation. This surface represents the three-dimensional shape of the underground reference interface. Subsequently, for each local reflection interface, the depth of the underground reference interface at each grid point is superimposed with the relative depth relationship of the interface at each grid point to obtain the absolute depth of the local reflection interface at each grid point. Based on the absolute depth and spatial coordinates, a subsequent continuous surface representing the local reflection interface is generated through spatial interpolation. Finally, the first continuous surface and all the subsequent continuous surfaces are combined and visualized according to their vertical spatial relationship to form a three-dimensional geological structure layer. Using relative hardness and relative density parameters as attribute values ​​and the spatial coordinates of grid points as the positioning basis, continuous attribute fields corresponding to hardness and density distributions are generated through two-dimensional spatial interpolation. Then, according to the preset color mapping rules, the continuous attribute fields are rendered to generate a two-dimensional engineering attribute layer expressed in the form of a planar diagram. The two-dimensional engineering attribute layer is used as the base plane map, and the planar projection of the three-dimensional geological structure layer is superimposed and associated with the base plane map in a unified geographic coordinate system to form an integrated model that supports synchronous query and visualization of spatial location.

Images