Terrain surveying and mapping method and system based on GNSS-RTK

Abstract

The invention provides a GNSS-RTK-based topographic mapping method and system, and the method comprises the steps: obtaining multi-source positioning signal data in a target region, the multi-source positioning signal data comprising a satellite original observation value, a receiver antenna phase center deviation parameter, and real-time dynamic differential correction information; performing multi-path interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression; inputting the phase observation sequence after interference suppression into a topographic feature calculation model, and extracting topographic elevation change feature parameters of the target area; performing spatial fusion processing based on the terrain elevation change characteristic parameters and reference elevation parameters in a preset terrain database to generate three-dimensional terrain curved surface data of the target area; and dynamically adjusting the topographic map updating frequency according to the difference degree between the three-dimensional topographic surface data and the historical topographic surveying and mapping result, and outputting a real-time topographic surveying and mapping result. According to the invention, precision, real-time performance and resource efficiency can be considered in complex terrain environment application.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing, and in particular, to a topographic surveying method and system based on GNSS-RTK. Background Art

[0002] In the prior art, a topographic surveying system based on GNSS-RTK usually adopts a processing mode of a single data source to obtain topographic information. For example, it relies on the carrier phase solution of satellite observations to generate elevation data, or post-processes and corrects the positioning results through static differential correction parameters. Such methods can achieve centimeter-level positioning accuracy in open and flat areas, but in environments with significant multipath effects such as urban building-dense areas and vegetation-covered areas, due to the accumulation of pseudorange noise caused by satellite signal reflection and diffraction, the elevation solution results show centimeter-level to decimeter-level fluctuation errors. Some improvement schemes attempt to suppress multipath interference by increasing the number of observed satellites or extending the observation time, but it is difficult to meet the timeliness requirements of real-time topographic monitoring. Another technical solution performs loose coupling processing on GNSS-RTK data and inertial navigation system data. Although it can briefly improve the positioning stability in a dynamic environment, due to the characteristic that the sensor errors diverge over time, it cannot guarantee the reliability of long-term continuous surveying and mapping. At the data fusion level, traditional methods linearly superimpose GNSS-RTK elevation data and historical topographic databases with fixed weights, without considering the influence of terrain slope changes and surface coverage types on the data fusion rules, resulting in elevation jumps at the junction of vegetation-covered areas and bare ground surfaces in the fused topographic surface. In addition, the existing topographic map update mechanism generally adopts a global update strategy with a preset time interval, without distinguishing the data change sensitivity of geologically stable areas and active areas, resulting in the coexistence of waste of computing resources and delay in updating key areas. The above technical defects make it difficult for traditional GNSS-RTK surveying and mapping methods to balance accuracy, real-time performance, and resource efficiency in complex terrain environments, restricting their popularization and application in high-dynamic scenarios such as disaster warning and engineering monitoring. Summary of the Invention

[0003] In view of this, the present invention provides a topographic surveying method and system based on GNSS-RTK. The technical solution of the embodiment of the present invention is realized as follows:

[0004] On the one hand, an embodiment of the present invention provides a terrain mapping method based on GNSS-RTK. The method includes: obtaining multi-source positioning signal data within a target area, where the multi-source positioning signal data includes satellite raw observations, receiver antenna phase center deviation parameters, and real-time kinematic differential correction information; performing multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression; inputting the phase observation sequence after interference suppression into a terrain feature calculation model to extract terrain elevation change feature parameters of the target area; performing spatial fusion processing based on the terrain elevation change feature parameters and reference elevation parameters in a preset terrain database to generate three-dimensional terrain surface data of the target area; dynamically adjusting the topographic map update frequency according to the difference degree between the three-dimensional terrain surface data and historical terrain mapping results, and outputting real-time terrain mapping results.

[0005] On the other hand, the present invention provides a terrain mapping system, including a memory and a processor. The memory stores a computer program that can run on the processor, and when the processor executes the program, the steps in the above method are implemented.

[0006] The terrain mapping method based on GNSS-RTK provided by the present invention can comprehensively utilize the spatial distribution characteristics of satellite signals, receiver hardware error parameters, and the correction ability of real-time differential data by obtaining multi-source positioning signal data within a target area, where the multi-source positioning signal data includes satellite raw observations, receiver antenna phase center deviation parameters, and real-time kinematic differential correction information; perform multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression, effectively eliminating the problem of multipath signal aliasing caused by urban building reflections and vegetation occlusion; input the phase observation sequence into a terrain feature calculation model to extract terrain elevation change feature parameters, and use a deep learning model to adaptively analyze the correlation between signal propagation laws and elevation changes under complex terrains; perform spatial fusion processing based on the terrain elevation change feature parameters and the reference elevation parameters of a preset terrain database to generate three-dimensional terrain surface data, and achieve spatial consistency fusion of multi-source terrain data through contour distribution correction, slope continuity optimization, and vegetation coverage compensation; dynamically adjust the update frequency according to the difference between the three-dimensional terrain surface data and historical results and output real-time terrain mapping results, and can intelligently switch between timed polling, event triggering, and real-time streaming update modes according to the actual terrain change rate, optimizing system resource allocation while ensuring mapping accuracy. In this way, satellite raw observations can provide high-frequency positioning reference information, receiver phase center deviation parameters can compensate for inherent hardware errors, and real-time kinematic differential correction information can eliminate environmental interferences such as atmospheric delays. The three work together to form a basis for high-precision data acquisition; multipath interference suppression processing significantly improves the signal-to-noise ratio in complex environments through joint carrier phase compensation and pseudorange noise dynamic segmentation; the terrain feature calculation model can accurately extract the microscopic change trend of surface elevation from phase fluctuation features through strongly correlated learning with lidar verification data during the training phase; spatial fusion processing solves the inherent defect of insufficient GNSS signal penetrability through a vegetation coverage area compensation mechanism, ensuring the continuity of elevation data in bare ground and vegetation areas; the difference-driven dynamic update strategy incorporates terrain stability indices and change sensitivity coefficients into the threshold determination logic, enabling the system to achieve millimeter-level change capture and minute-level data refresh in scenarios such as geological disaster early warning and engineering construction monitoring, ultimately achieving comprehensive improvements in mapping accuracy, real-time performance, and resource efficiency in complex terrain environments and expanding the practical boundaries of GNSS-RTK technology in the field of dynamic terrain monitoring. Description of the Drawings

[0007] Figure 1 It is a schematic flowchart of the implementation of a terrain mapping method based on GNSS-RTK provided by an embodiment of the present invention.

[0008] Figure 2 It is a schematic diagram of the hardware entity of a terrain mapping system provided by an embodiment of the present invention. Detailed Implementation Modes

[0009] An embodiment of the present invention provides a topographic surveying method based on GNSS-RTK, which can be executed by a processor of a topographic surveying system. Among them, the topographic surveying system may refer to a device with data processing capabilities such as a server, a laptop computer, a tablet computer, and a desktop computer.

[0010] Figure 1 It is a schematic flowchart of the implementation of a topographic surveying method based on GNSS-RTK provided by an embodiment of the present invention. As Figure 1 shown, the method includes the following steps:

[0011] Step S100: Obtain multi-source positioning signal data within the target area. The multi-source positioning signal data includes satellite raw observations, receiver antenna phase center deviation parameters, and real-time kinematic differential correction information.

[0012] In step S100, omnidirectional data collection and integration of multiple signal sources within the target area are performed. The multi-source positioning signal data refers to a set of positioning information with complementary characteristics obtained through different technical means. Its core components include satellite raw observations, receiver antenna phase center deviation parameters, and real-time kinematic differential correction information. Satellite raw observations refer to the raw signal measurement data directly captured by a Global Navigation Satellite System (GNSS) receiver, specifically including carrier phase observations, pseudorange measurements, Doppler frequency shift data, and satellite ephemeris parameters transmitted by the satellite. Among them, the carrier phase observation is the sum of the integer part and the fractional part of the satellite signal carrier cycle recorded by the receiver. The pseudorange measurement is the geometric distance between the satellite and the receiver calculated through the signal propagation time. The Doppler frequency shift data reflects the relative motion state between the receiver and the satellite. The satellite ephemeris parameters include satellite orbital position, clock bias, and health status information. The receiver antenna phase center deviation parameter refers to the spatial offset between the physical center and the signal phase center of the receiver antenna. This parameter is determined by the antenna design characteristics and is specifically manifested as the position correction value of the antenna in a three-dimensional coordinate system, which needs to be obtained through a precise calibration experiment. The real-time kinematic differential correction information refers to the real-time error correction data provided by a reference station network, specifically including atmospheric delay correction parameters (ionospheric and tropospheric delays), satellite orbit error correction values, and satellite clock bias correction values. These data are transmitted to the rover receiver through a wireless communication link to eliminate the influence of common error sources on the positioning accuracy.

[0013] Step S200: Perform multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after suppressing interference.

[0014] The implementation process of step S200 is used to eliminate the interference of multipath effects caused by signal reflection on positioning data. Multipath interference refers to the superposition of multiple signals formed by the reflection of GNSS signals from obstacles such as the ground and buildings during the propagation process, which will cause periodic fluctuation errors in the carrier phase and pseudorange measured by the receiver. Multipath interference suppression processing adopts dynamic signal separation and noise suppression technology, which specifically includes the following core operations: First, the carrier phase measurement sequence in the original satellite observation value is analyzed in time and frequency to identify the periodic fluctuation pattern caused by the multipath effect; secondly, a multipath noise model is constructed based on the receiver motion state parameters, and the real signal and the reflected signal components are separated by an adaptive filtering algorithm; finally, the residual noise is corrected for spatial correlation in combination with the horizontal positioning error parameters in the real-time dynamic differential correction information.

[0015] In the specific implementation, the generation of the phase observation sequence after interference suppression requires the processing of multi-source positioning signal data in stages. First, the carrier phase measurement value is compensated for the antenna phase center deviation, and the original phase observation value is corrected in three dimensions using the receiver antenna phase center deviation parameter to eliminate the systematic deviation caused by the physical characteristics of the antenna. Secondly, for the multipath noise in the pseudorange measurement value, the sliding time window analysis method is used to divide the pseudorange sequence into multiple time segments. By calculating the standard deviation and mean of the pseudorange in each segment, the mutation noise points exceeding the dynamic threshold are identified and marked as outliers for elimination. Subsequently, the compensated carrier phase sequence and the filtered pseudorange measurement value are jointly adjusted. The adjustment model introduces the receiver motion state parameters (such as velocity and acceleration) as constraints, and the optimal phase observation solution is solved by the least squares optimization algorithm. In this process, the horizontal direction error parameters in the real-time dynamic differential correction information are used to construct the correlation matrix of the multipath effect and pseudorange noise. The noise components related to the multipath are extracted by matrix decomposition technology and stripped from the original observation sequence. The final output interference-suppressed phase observation sequence is a high-precision carrier phase data set that has undergone bias compensation, noise removal, and adjustment optimization. Its data format is a continuous phase observation value sequence with timestamp alignment.

[0016] Step S300: Input the interference-suppressed phase observation sequence into the terrain feature solution model to extract the terrain elevation change characteristic parameters of the target area.

[0017] Step S300 converts the phase observation data into elevation feature parameters characterizing the terrain undulation through a mathematical model. The terrain feature resolution model is a multi-modal data fusion architecture based on deep learning. Its input is the phase observation sequence after interference suppression, and the output is the elevation change gradient, surface curvature, and slope distribution parameters within the target area. This model first analyzes the spatial correlation features in the phase observation sequence through the signal feature extraction layer, including the amplitude of carrier phase fluctuations, the distribution pattern of satellite elevation angles, and the spatial distribution characteristics of the receiver's movement trajectory. Subsequently, the elevation mapping layer models the non-linear relationship between the signal features and the terrain elevation to generate a preliminary elevation estimate. Finally, the terrain surface generation layer smooths and optimizes the preliminary estimate to eliminate local noise and enhance the expression of terrain continuity. In a specific implementation, the extraction of terrain elevation change feature parameters requires processing the phase observation sequence in stages. First, the signal feature extraction layer of the terrain feature resolution model uses a Convolutional Neural Network (CNN) structure to extract time-domain and spatial-domain features from the input phase observation sequence. The time-domain features include the short-term fluctuation trend and long-term change period of the phase observation values, which are realized by sliding a one-dimensional convolutional kernel. The spatial-domain features are extracted by analyzing the spatial geometric distribution of different satellite signals (such as the combined relationship between satellite elevation angles and azimuth angles), and a Graph Convolutional Network (GCN) is specifically used to model the spatial topology graph formed by the satellite-receiver. Secondly, the elevation mapping layer fuses the extracted signal features with the preset terrain prior knowledge (such as the elevation distribution law of typical landforms) and generates preliminary elevation change parameters through a fully connected neural network, including elevation increment, elevation change rate, and local curvature. Finally, the terrain surface generation layer uses a surface optimization algorithm based on physical constraints to jointly optimize the preliminary elevation parameters and the terrain continuity prior (such as the consistency of elevation gradients in adjacent areas) and outputs the terrain elevation change feature parameters. The feature parameters are stored in a rasterized data format, and each raster cell contains an elevation value, an elevation change direction, and a confidence index, which are used for subsequent 3D terrain reconstruction.

[0018] Step S400: Based on the terrain elevation change feature parameters and the reference elevation parameters in the preset terrain database, perform spatial fusion processing to generate the 3D terrain surface data of the target area.

[0019] Step S400 is used to fuse the elevation change features calculated in real time with the historical reference terrain data to construct a high-precision three-dimensional terrain model. The preset terrain database is a standardized data set containing the historical surveying and mapping results of the target area, and its reference elevation parameters include contour distribution data, slope parameters, and surface cover type parameters. The contour distribution data is a contour vector map obtained through aerial photogrammetry or lidar scanning. The slope parameter is slope raster data calculated based on a Digital Elevation Model (DEM). The surface cover type parameter is a land cover class label (such as vegetation, water area, bare soil) obtained through remote sensing image classification. The core technologies for spatial fusion processing include elevation data superposition, terrain continuity correction, and surface cover compensation.

[0020] In specific implementation, generating three-dimensional terrain surface data requires performing multi-level spatial data fusion operations. First, the elevation increment data in the terrain elevation change feature parameters is spatially superimposed with the contour distribution data in the preset terrain database. The superposition process uses raster-vector conversion technology to convert the contour vector data into a raster format with the same resolution as the elevation increment data, and generates the superimposed elevation distribution map through pixel-by-pixel addition operations. Second, based on the slope parameter, the terrain continuity of the superimposed elevation distribution map is corrected: by calculating the elevation gradient of adjacent raster cells, elevation mutation regions that do not conform to the natural terrain evolution law (such as steep cliff illusions caused by data noise) are identified, and the anisotropic diffusion algorithm is used to smooth the mutation regions. Subsequently, according to the surface cover type parameter, elevation compensation for the vegetation-covered areas in the corrected elevation distribution map is performed: for areas with dense vegetation such as forests and shrubs, a preset vegetation layer thickness compensation value (such as 0.5 meters to 3 meters) is added to the original elevation value, and this compensation value is dynamically adjusted through an empirical relationship library of vegetation type and height. Finally, the compensated elevation distribution data is converted into a gridded elevation point set, and a continuous three-dimensional terrain surface data is generated through a cubic spline interpolation algorithm.

[0021] Step S500: Dynamically adjust the topographic map update frequency according to the difference between the three-dimensional terrain surface data and the historical topographic survey results, and output the real-time topographic survey results.

[0022] The difference degree refers to the quantitative index of the elevation deviation at the same spatial position between the current three-dimensional terrain surface data and the historical topographic survey results. Its calculation process includes the statistics of the absolute value of the elevation deviation, the weighting of terrain stability, and the correction of the sensitivity coefficient. The strategy for dynamically adjusting the topographic map update frequency is based on a difference degree threshold classification mechanism: when the difference degree is lower than the first threshold, a timed polling update mode is adopted to update the data at a fixed time interval (such as 24 hours); when the difference degree is between the first threshold and the second threshold, it switches to an event-triggered update mode, and only starts a local update when significant terrain changes are detected; when the difference degree exceeds the second threshold, a real-time streaming update mode is enabled, and the entire region data is processed and verified in real time through parallel computing nodes.

[0023] In specific implementation, multiple levels of data processing and scheduling operations need to be performed to output the real-time topographic survey results. First, when calculating the difference degree, the three-dimensional terrain surface data and the historical topographic survey results need to be spatially matched: the elevation point sets of the two are unified to the same coordinate system through coordinate transformation, and the nearest neighbor interpolation method is used to align the spatial resolution. Secondly, for the matched elevation point pairs, the absolute deviation between the current elevation value and the historical value is calculated point by point, and weighted averaging is performed according to the terrain stability coefficient of the region where the point is located (for example, the coefficient of the geologically active region is 0.3, and the coefficient of the flat region is 0.7) to obtain a preliminary difference degree index. Subsequently, the difference degree is corrected by the terrain change sensitivity coefficient: the sensitivity coefficient is calculated from the geological structure activity monitoring data of the target region (such as the surface displacement rate, earthquake frequency), and its function is to amplify the difference degree weight of the tectonically active region. Based on the corrected difference degree value, the system automatically selects the update mode and triggers the corresponding processing flow. In the real-time streaming update mode, the three-dimensional terrain surface data is divided into multiple spatial data blocks, and each data block is assigned to an independent parallel computing node for outlier detection and interpolation compensation. The processed data blocks are generated into the final result with spatio-temporal consistency markers through spatial stitching and timestamp alignment. The spatio-temporal consistency markers include the data acquisition time, spatial range, and processing pipeline number, and are bound to the terrain data through a hash algorithm to ensure data integrity and traceability.

[0024] As an implementation method, in step S200, multi-path interference suppression processing is performed on the multi-source positioning signal data to generate a phase observation sequence after suppressing interference, which may specifically include:

[0025] Step S210: Obtain the carrier phase measurement value and the pseudorange measurement value in the satellite raw observation values, and construct a carrier phase fluctuation sequence and a pseudorange noise distribution sequence respectively.

[0026] The satellite raw observations are the unprocessed raw signal measurement results output by a global navigation satellite system receiver, specifically including two types of core data: carrier phase measurements and pseudorange measurements. The carrier phase measurement refers to the phase difference of the satellite carrier signal propagation path recorded by the receiver. Its value consists of an integer cycle part and a fractional cycle part, reflecting the precise length change of the signal propagation path, with millimeter-level accuracy but having the problem of integer ambiguity. The pseudorange measurement refers to the geometric distance between the satellite and the receiver calculated based on the signal propagation time. Its value is obtained through modulation code phase measurement, with the characteristic of being unambiguous but vulnerable to multipath effects and ionospheric delays. The process of constructing the carrier phase fluctuation sequence is as follows: extract the carrier phase measurements of each satellite from the satellite raw observations in a continuous time period, arrange them in chronological order to form a time series, and eliminate the common errors of the receiver clock error and the satellite clock error through differential operation to generate sequence data reflecting the short-term fluctuations of the carrier phase. The process of constructing the pseudorange noise distribution sequence is as follows: classify the pseudorange measurements in the same time period according to the satellite number, perform time serialization processing on the pseudorange measurements of each satellite, and separate the low-frequency trend term and the high-frequency noise term through sliding mean filtering, and extract the high-frequency noise term to construct the pseudorange noise distribution sequence. During this process, it is necessary to ensure that the timestamps of the carrier phase fluctuation sequence and the pseudorange noise distribution sequence are strictly aligned to support subsequent joint analysis.

[0027] Step S220: Perform antenna phase deviation compensation on the carrier phase fluctuation sequence based on the receiver antenna phase center deviation parameter to generate the compensated carrier phase sequence.

[0028] Step S220 is used to eliminate the systematic carrier phase deviation introduced by the hardware characteristics of the receiver antenna. The receiver antenna phase center deviation parameter refers to the three-dimensional spatial offset between the physical geometric center of the antenna and its electromagnetic wave phase center. This parameter is determined by the radiation characteristics of the antenna model and needs to be obtained through laboratory calibration or calibration files provided by the manufacturer. The operation process of antenna phase deviation compensation is as follows: First, call the corresponding phase center deviation parameter from the pre-stored antenna parameter database according to the receiver antenna model. This parameter is usually stored in the form of three-dimensional offsets in the east-north-up (ENU) coordinate system, denoted as ΔX, ΔY, and ΔZ. Second, combine each measurement value in the carrier phase fluctuation sequence with its corresponding satellite azimuth and elevation angle information, and project the phase center deviation parameter to the satellite signal propagation direction through a coordinate transformation model to calculate the direction-dependent phase deviation correction amount. Finally, subtract the calculated phase deviation correction amount from the original carrier phase fluctuation sequence one by one to generate the compensated carrier phase sequence. For example, for the carrier phase fluctuation value L1(ti) of satellite G01 at time ti, its phase deviation correction amount ΔL1(ti) can be calculated by the formula ΔL1(ti)=(ΔX·cosα + ΔY·sinα)·sinβ + ΔZ·cosβ, where α is the azimuth angle of satellite G01 at time ti and β is the elevation angle. Through this compensation operation, the hardware-related errors in the carrier phase fluctuation sequence are effectively suppressed, and the compensated carrier phase sequence can more truly reflect the physical changes in the signal propagation path.

[0029] Step S230: Establish an association model between the multipath effect influence factor and the pseudorange noise distribution sequence according to the horizontal direction positioning error parameter in the real-time kinematic differential correction information.

[0030] The horizontal direction positioning error parameter in the real-time kinematic differential correction information refers to the horizontal direction (eastward and northward) positioning residuals calculated by the reference station network through real-time differential technology. This parameter reflects the remaining positioning deviation of the reference station receiver after eliminating the common errors, and its value is affected by both the multipath effect and the residual error of atmospheric delay. The multipath effect influence factor is defined as a quantitative index describing the intensity of multipath interference, and its value is related to factors such as the length of the reflected signal path and the material of the reflecting surface. The process of establishing the association model includes the following steps: First, extract the time series of the horizontal direction positioning error parameter from the real-time kinematic differential correction information, and perform a synthesis calculation on the eastward error sequence E(t) and the northward error sequence N(t) to generate the horizontal direction error modulus sequence M(t)=√(E(t) 2 +N(t) 2);Secondly, synchronize the pseudorange noise distribution sequence P(t) and the horizontal direction error modulus sequence M(t) in time. Finally, use the multiple regression analysis method to construct the correlation model between P(t) and M(t), and its specific form is P(t)=k·M(t)+C+ε, where k is the proportionality coefficient, C is the constant term, and ε is the random noise. Through this model, the multipath effect influence factor k can quantitatively characterize the contribution degree of the horizontal direction positioning error to the pseudorange noise.

[0031] Step S240: Perform dynamic threshold segmentation on the correlation model through a sliding time window to identify the mutation noise points in the pseudorange noise distribution sequence.

[0032] Step S240 is used to detect the abnormal mutation points in the pseudorange noise distribution sequence to support noise rejection. The sliding time window refers to a data segmentation mechanism that slides along the time axis at a fixed step length, and its window length is dynamically set according to the signal sampling rate and the time-varying characteristics of multipath interference. For example, a 5-second window length can be used at a 1Hz sampling rate. The operation process of dynamic threshold segmentation is as follows: First, divide the pseudorange noise distribution sequence into multiple subsequences according to the sliding time window, analyze the statistical characteristics of the noise data in each subsequence, and calculate its mean μ and standard deviation σ; Secondly, set dynamic thresholds for each subsequence based on the statistical characteristics. For example, set the threshold to μ±3σ; Finally, detect the noise points in the subsequence that exceed the dynamic threshold and mark them as mutation noise points. In this process, it is necessary to optimize the threshold setting strategy in combination with the output results of the correlation model: for the subsequences with a higher multipath effect influence factor k, appropriately reduce the threshold to enhance the detection sensitivity of mutation points; conversely, for the subsequences with a lower k value, increase the threshold to avoid false detection.

[0033] Step S250: Perform joint adjustment calculation on the compensated carrier phase sequence and the pseudorange measurement values after removing the mutation noise points to generate a phase observation sequence with suppressed interference.

[0034] Step S250 aims to fuse the compensated high-precision carrier phase data and the filtered pseudorange observation data, and generate an anti-interference optimized phase observation sequence through the joint adjustment algorithm. The joint adjustment calculation refers to incorporating the measurement values of different observation types into a unified mathematical model for overall solution, and its advantage is that it can use the high-precision characteristics of the carrier phase and the ambiguity-free characteristics of the pseudorange to complement and improve the solution accuracy. The specific implementation process includes the following steps: First, establish a joint adjustment model of the carrier phase observation equation and the pseudorange observation equation, where the carrier phase equation is expressed as λ·Φ=ρ+c·(dT-dt)+Trop+Iono+ε Φ , and the pseudorange equation is expressed as P=ρ+c·(dT-dt)+Trop+Iono+ε P, where λ is the carrier wavelength, Φ is the measured carrier phase, ρ is the geometric distance, dT and dt are the receiver clock error and satellite clock error respectively, Trop is the tropospheric delay, Iono is the ionospheric delay, and ε Φ and ε P are the observation noises of the carrier phase and pseudorange respectively; Secondly, substitute the compensated carrier phase sequence and the pseudorange measurement value after removing the mutation noise points into the combined adjustment model, introduce the receiver motion state parameters (such as velocity, acceleration) as constraint conditions, and use the least squares estimation algorithm to solve the optimal integer ambiguity combination and position correction amount; Finally, substitute the fixed value of the integer ambiguity obtained by the solution back into the carrier phase observation equation to generate the phase observation sequence after suppressing interference.

[0035] As an implementation manner, the method further includes the training process of the terrain feature solution model, which may specifically include the following steps:

[0036] Step S301: Collect multiple groups of training data in historical topographic surveying tasks. Each group of training data includes a phase observation sequence sample after multipath interference suppression processing and the corresponding lidar elevation verification data of the phase observation sequence sample.

[0037] Step S301 is used to construct a high-quality training dataset, and its core lies in obtaining multi-source terrain observation data pairs with spatio-temporal consistency. Historical terrain mapping tasks refer to past completed and verified mapping projects. The data contains the complete original observation data stream of the global navigation satellite system, and is accompanied by high-precision lidar elevation verification data. At the same time, multi-path interference suppression processing has been completed and a standardized phase observation sequence has been generated. The phase observation sequence sample after multi-path interference suppression processing refers to the carrier phase time series data processed by the method of Step S200, and its data format is a matrix of phase observation values arranged by time stamps. The matrix dimensions include satellite number, observation epoch, and phase compensation value. The lidar elevation verification data refers to the digital elevation model (DEM) generated after processing the three-dimensional point cloud data collected by an airborne or ground lidar system. Its spatial resolution needs to reach the sub-meter level (such as 0.5 meters), the elevation accuracy is better than ±5 cm, and it is strictly aligned with the phase observation sequence sample in time and space. In specific implementation, the following operations need to be performed for data collection: First, screen task records that meet the time span requirements from the historical mapping project database to ensure that the selected tasks cover diverse terrain types (such as mountains, plains, urban construction areas); Second, perform multi-path interference suppression processing of Step S200 on the original observation data of the global navigation satellite system in each task to generate phase observation sequence samples; At the same time, call the lidar point cloud data of the corresponding task, and generate standardized lidar elevation verification data through point cloud filtering, ground point classification, and rasterization processing; Finally, perform spatio-temporal alignment processing on the phase observation sequence sample and the lidar elevation verification data, specifically including timestamp matching (aligning the lidar acquisition time with the epoch time of the phase observation sequence) and coordinate system unification (converting both to the same plane coordinate system and elevation datum).

[0038] Step S302: Construct a deep neural network model. The deep neural network model includes a signal feature extraction layer, an elevation mapping layer, and a terrain surface generation layer.

[0039] The deep neural network model is, for example, a three - stage cascaded structure. Among them, the signal feature extraction layer is responsible for parsing multi - level spatio - temporal features from the phase observation sequence. Its input is the phase observation matrix of multiple satellites and multiple epochs, and the output is a feature vector containing phase fluctuation patterns, satellite geometric distribution characteristics, and receiver motion states; the elevation mapping layer is used to establish a non - linear mapping relationship between signal features and terrain elevation, and it realizes the conversion from the feature space to elevation parameters through a fully - connected network and an attention mechanism; the terrain surface generation layer focuses on reconstructing discrete elevation estimates into a continuous terrain surface, and uses a physically - constrained interpolation algorithm and a smoothing filter to eliminate local errors. In specific implementation, model construction needs to be realized in modules. The signal feature extraction layer can adopt a hybrid architecture of a convolutional neural network and a long - short - term memory network. Among them, a one - dimensional convolutional kernel (with a size of 3×1) is used to extract local fluctuation features of the phase observation sequence, and the LSTM unit is used to capture the dynamic pattern of the satellite elevation distribution changing with time; secondly, the elevation mapping layer can be composed of four fully - connected networks, and the number of neurons in each layer decreases sequentially (such as 512→256→128→64). The activation function uses the Rectified Linear Unit (ReLU), and a residual connection is introduced in the output layer to accelerate convergence; finally, the terrain surface generation layer integrates a de - convolutional network and a Thin Plate Spline (TPS) interpolation algorithm. The de - convolutional network is responsible for up - sampling the feature vector to the target spatial resolution, and the TPS algorithm generates a smooth surface according to the control point constraints.

[0040] Step S303: Input the phase observation sequence sample into the signal feature extraction layer to extract phase fluctuation features, satellite elevation distribution features, and receiver motion state features.

[0041] The phase fluctuation feature refers to the short-term fluctuation pattern of the carrier phase observation value changing with time. Its value reflects the surface reflection characteristics and the minute changes in the signal propagation path, and is obtained by sliding a one-dimensional convolution kernel in the time dimension. The satellite elevation distribution feature refers to the statistical characteristics of the elevations of the satellites participating in the solution changing with time in the sky. Its value affects the multipath effect intensity and the signal-to-noise ratio, and the spatial correlation feature is extracted by constructing a satellite-epoch elevation matrix and applying a Graph Convolutional Network (GCN). The receiver motion state features include the receiver speed, acceleration, and the curvature parameter of the motion trajectory. Its value is derived from the dynamic positioning result of the global navigation satellite system, and the time evolution law of the motion state is captured by an LSTM network. In the specific implementation, feature extraction needs to be processed channel by channel: First, the phase observation sequence sample is reshaped into a three-dimensional tensor form [number of satellites × number of epochs × 1 channel], and is input to the CNN module of the signal feature extraction layer. The multi-scale phase fluctuation feature is extracted through three convolutional layers (the number of filters is 16, 32, and 64 respectively), and a feature map with a size of [number of satellites × number of epochs × 64 channels] is output; Second, the satellite elevation distribution feature is obtained by constructing an elevation-time matrix (with a size of [number of satellites × number of epochs]), applying the GCN model to analyze the spatial geometric relationship between satellites, and generating a 32-dimensional elevation distribution feature vector; At the same time, the receiver motion state features (including the eastward speed, northward speed, elevation direction speed, and three-axis acceleration) are arranged in chronological order of epochs as a time series, input to the LSTM network (the number of hidden layer units is 64), and a 128-dimensional motion state feature vector is output.

[0042] Step S304: Fuse the phase fluctuation feature and the satellite elevation distribution feature to generate a spatial correlation feature vector, and input the spatial correlation feature vector into the elevation mapping layer for non-linear transformation to obtain the non-linear transformation result.

[0043] Step S304 integrates the time-domain characteristics of phase fluctuations and the spatial distribution characteristics of satellite elevation angles to construct spatial correlation features reflecting terrain elevation changes. The feature fusion adopts a strategy combining cascading and attention mechanism: First, the phase fluctuation feature matrix (size [number of satellites × number of epochs × 64 channels]) is subjected to max pooling along the satellite dimension and compressed into a time-domain feature sequence of [number of epochs × 64 channels]; Second, the satellite elevation distribution feature vector is extended into a sequence form of [number of epochs × 32 channels], and is concatenated with the phase fluctuation time-domain feature sequence in the channel dimension to generate a preliminary fusion feature of [number of epochs × 96 channels]; Finally, a self-attention mechanism is introduced to calculate the correlation weights between different epoch features, and a spatial correlation feature vector of [number of epochs × 96 channels] is generated through weighted summation. The non-linear transformation process of the elevation mapping layer is as follows: The spatial correlation feature vector is input into a four-layer fully-connected network in epoch order. The first layer maps the 96-channel feature to a 512-dimensional high-dimensional space, the second layer reduces it to 256 dimensions, the third layer further compresses it to 128 dimensions, and the fourth layer outputs a 64-dimensional non-linear transformation result. During this process, batch normalization and ReLU activation functions are applied after each fully-connected network, and a residual connection is added between the third and fourth layers to suppress gradient vanishing.

[0044] Step S305: Weightedly superimpose the receiver motion state features and the non-linear transformation results to generate preliminary terrain elevation parameters.

[0045] Step S305 is used to fuse dynamic motion information and spatial elevation features to improve the robustness of terrain calculation. The weighted superposition adopts a method combining feature concatenation and adaptive weighting: First, the receiver motion state feature vector (size

[128] ) is replicated and extended along the epoch dimension to the same size as the non-linear transformation result ([100 × 128]); Second, the extended motion state features and the non-linear transformation results (size [100 × 64]) are concatenated in the feature dimension to generate a fusion feature matrix of [100 × 192]; Subsequently, the fusion features are weighted in channels by a learnable weight matrix, which is dynamically generated by a fully-connected network. Its input is the motion state feature and the non-linear transformation result of the current epoch, and the output is a 192-dimensional weight vector; Finally, the weighted fusion features are input into a fully-connected layer (output dimension is 1) to generate the elevation change parameters corresponding to each epoch, and after aggregating by spatial position, a preliminary terrain elevation parameter matrix is formed.

[0046] Step S306: Perform surface smoothing on the preliminary terrain elevation parameters through the terrain surface generation layer to output the predicted terrain surface data.

[0047] Step S306 is used to reconstruct the discrete elevation parameters into surface data that conforms to the constraints of terrain continuity and smoothness. The surface smoothing process adopts a combined algorithm of Thin Plate Spline (TPS) and Anisotropic Diffusion Filter: First, each grid cell in the preliminary terrain elevation parameter matrix is regarded as a control point, and the TPS algorithm is applied to generate a continuous surface, whose energy function minimizes the bending energy to achieve global smoothness; Second, an anisotropic diffusion coefficient matrix is constructed based on the land cover type data (such as vegetation, bare soil, water area), and edge-preserving smoothing is performed on the basis of the TPS surface, that is, the smoothing intensity is enhanced in homogeneous regions (such as flat farmland), and details are retained in heterogeneous regions (such as building edges). During this process, the terrain surface generation layer also integrates a deconvolution network to upsample the low-resolution elevation parameters. For example, the initial elevation matrix with a resolution of 0.5 meters is upsampled to 0.2 meters, and then high-resolution terrain surface data is output through TPS and diffusion filtering.

[0048] Step S307: Calculate the elevation difference loss value between the predicted terrain surface data and the lidar elevation verification data, and use the backpropagation algorithm to optimize the parameters of the deep neural network model until the elevation difference loss value is less than the preset threshold.

[0049] Step S307 optimizes the model parameters through supervised learning to minimize the solution error. The elevation difference loss value is defined as the mean of the squared elevation differences between the predicted terrain surface data and the lidar elevation verification data at the same spatial position, and its calculation formula is: ; where N is the number of valid grid points, is the predicted value, is the lidar verification value. The backpropagation algorithm uses the Adaptive Moment Estimation (Adam) optimizer, the initial learning rate is set to 1e-4, and a learning rate decay strategy (decay factor 0.5 every 10 epochs) is applied. During the training process, the following operations need to be performed: First, input the batch data (such as 32 groups of training samples) into the deep neural network model, and generate predicted terrain surface data through forward propagation; Second, calculate the elevation difference loss value of the current batch, and calculate the gradients of the parameters of each layer through backpropagation; Finally, update the model parameters and repeat the iteration until the loss value converges to the preset threshold (such as 1e-3 meters 2 ).

[0050] As an implementation, step S400, based on the terrain elevation change characteristic parameters and the reference elevation parameters in the preset terrain database, performs spatial fusion processing to generate the three-dimensional terrain surface data of the target area, which may specifically include:

[0051] Step S410: Extract the reference elevation parameters corresponding to the target area from the preset terrain database. The reference elevation parameters include contour line distribution data, slope parameters, and surface cover type parameters.

[0052] Step S410 is used to call the basic geographic information data in the preset terrain database to support terrain fusion calculation. The preset terrain database is a standardized geospatial database containing the historical surveying and mapping results of the target area. Its reference elevation parameters are generated by standardizing multi-source geographic data, and specifically may include contour line distribution data, slope parameters, and surface cover type parameters. The contour line distribution data is a vector data set of contour lines obtained through aerial photogrammetry or lidar scanning technology. Each contour line represents the continuous spatial distribution of the same elevation value, and its data attributes include elevation value, coordinate point sequence, and accuracy level information. The storage format is a Shapefile or GeoJSON file conforming to the standards of the Open Geospatial Consortium (OGC). The slope parameter is a raster data set generated based on a Digital Elevation Model (DEM) through a slope calculation algorithm. Each raster cell stores the slope value (in degrees or percentage), reflecting the surface inclination degree. The third-order inverse distance squared difference method is used in the calculation to balance the calculation efficiency and accuracy. The surface cover type parameter is a raster data of land cover classes obtained through multi-spectral remote sensing image classification. Its classification system follows international geographic information standards (such as ISO 19144), including classes such as vegetation, water area, bare soil, and buildings. Each raster cell stores the class code and classification confidence. In specific implementation, the extraction operation needs to execute the following process: perform a spatial query from the preset terrain database according to the boundary coordinates of the target area (such as the minimum bounding rectangle or the geographic fence polygon), and filter out the reference elevation parameter data blocks that are completely contained or partially overlapped; secondly, perform format conversion and coordinate system unification processing on the filtered data blocks, convert the contour line distribution data from the vector format to the same raster format as the slope parameter and the surface cover type parameter (such as GeoTIFF), and convert all data to the same planar coordinate system (such as WGS84 UTM) and elevation datum (such as the EGM2008 geoid); finally, perform resolution alignment on the converted data, and adjust the spatial resolution of the slope parameter and the surface cover type parameter to be the same as that of the contour line distribution data (such as 0.5 meters) through bilinear interpolation.

[0053] Step S420: Perform spatial overlay on the elevation increment data in the terrain elevation change feature parameters and the contour line distribution data to generate an overlaid elevation distribution map.

[0054] Step S420 is used to integrate the real-time calculated elevation change data and the historical reference terrain data to generate an updated elevation distribution map. The elevation increment data in the terrain elevation change characteristic parameters refers to the rasterized elevation change amount output by the terrain feature calculation model in step S300. Its value represents the elevation change value of each position in the target area relative to the historical reference during the current calculation period (positive value for uplift and negative value for subsidence). The data format is a floating-point raster consistent with the resolution of the reference elevation parameter. The spatial overlay operation uses the raster algebra operation method. The specific process is as follows: First, convert the contour line distribution data from the vector format to the same raster format as the elevation increment data. During the conversion process, assign the elevation value of the corresponding contour line to each raster cell, and fill the area not directly covered by the contour line through linear interpolation; Second, perform a pixel-by-pixel addition operation on the elevation increment data and the rasterized contour line data to generate the overlaid elevation distribution map. Its mathematical expression can be, for example, H new(x,y) =H base(x,y) +ΔH(x,y), where H base(x,y) is the reference elevation value and ΔH(x,y) is the elevation increment value; Finally, perform data validity verification on the overlay result to eliminate outliers (such as mutation values outside the reasonable elevation range) caused by coordinate offset or resolution mismatch.

[0055] Step S430: Correct the terrain continuity of the overlaid elevation distribution map according to the slope parameter to eliminate the elevation mutation area and obtain the corrected elevation distribution map.

[0056] Step S430 is used to eliminate the elevation mutation phenomenon caused by data noise or fusion error, ensuring that the terrain surface conforms to the continuity law of natural landforms. The terrain continuity correction adopts the Slope-Constrained Anisotropic Diffusion Filter algorithm. Its principle is to dynamically adjust the filtering intensity according to the slope parameter, retain the detail features in the steep slope area, and enhance the smoothing effect in the gentle slope area. The specific process includes: First, calculate the elevation gradient of each grid cell in the superimposed elevation distribution map (using the Sobel operator) to generate the horizontal gradient Gx and the vertical gradient Gy; Second, calculate the anisotropic diffusion coefficient matrix according to the slope parameter. The diffusion coefficient D(x,y) is defined as D(x,y)=exp(-|S(x,y)| / k), where S(x,y) is the slope value and k is the slope sensitivity factor (usually set to 30 degrees). This formula ensures that the diffusion coefficient approaches zero in the area where the slope is greater than k (such as cliffs) to suppress smoothing, while the diffusion coefficient increases in the area where the slope is less than k to promote smoothing; Finally, apply the explicit iteration method to perform diffusion filtering on the elevation distribution map. The number of iterations is dynamically set according to the severity of the elevation mutation (such as 5 to 20 times). The formula for updating the elevation value in each iteration is H t+1 (x,y)=H t (x,y)+λ·[D(x,y)·(Gx 2 +Gy 2 )], where λ is the iteration step size (usually set to 0.25).

[0057] Step S440: Perform elevation compensation for the vegetation-covered area on the corrected elevation distribution map based on the land surface cover type parameter to generate the compensated elevation distribution data.

[0058] Step S440 is used to correct the elevation underestimation error caused by the signal penetration effect in the vegetation-covered area. The elevation compensation for the vegetation-covered area adopts a category-related compensation strategy, and its compensation value is dynamically determined according to the empirical relationship library of vegetation type and height. The specific implementation steps are: First, extract the vegetation category subclasses (such as coniferous forest, broad-leaved forest, shrub) from the land surface cover type parameter, and configure a preset vegetation layer thickness compensation value for each subclass (such as 0.8 meters for coniferous forest compensation, 1.2 meters for broad-leaved forest compensation, 0.3 meters for shrub compensation); Second, spatially align the corrected elevation distribution map with the vegetation cover type grid, and perform elevation compensation operations on each grid cell in the vegetation-covered area; Finally, perform edge transition processing on the compensated elevation data, and use the Gaussian filtering method to smoothly transition the elevation values at the boundary between the vegetation and non-vegetation areas to avoid sudden changes in the compensation value.

[0059] Step S450: Convert the compensated elevation distribution data into a gridded elevation point set, and generate three-dimensional terrain surface data through the cubic spline interpolation algorithm.

[0060] Step S450 is used to convert the elevation data in raster format into a continuous surface model to support three-dimensional visualization and analysis. The gridded elevation point set refers to a set of discrete elevation points arranged at regular intervals, and its data format is a sequence of triples containing longitude, latitude, and elevation values. The grid interval is the same as the resolution of the compensated elevation distribution data (such as 0.5 meters). The conversion process includes the following operations: First, traverse the compensated elevation distribution raster data row by row and column, generate corresponding elevation points for each valid raster cell, and determine the coordinates based on the raster origin coordinates, pixel size, and row and column indices; Second, perform a topological check on the elevation point set, remove isolated points caused by data missing or invalid values, and fill in the void areas through the nearest neighbor interpolation method. The cubic spline interpolation algorithm uses a bivariate cubic polynomial function to fit the discrete elevation points. Boundary conditions (such as natural spline boundary conditions) need to be set during the interpolation process to ensure smooth transition at the surface edges. The finally output three-dimensional terrain surface data is stored in the format of a Triangulated Irregular Network (TIN), and each triangle vertex contains accurate coordinates and elevation values, and is accompanied by interpolation accuracy metadata, which can be directly imported into a Geographic Information System (GIS) platform for three-dimensional rendering or terrain analysis.

[0061] As an implementation, in step S500, dynamically adjust the topographic map update frequency according to the difference between the three-dimensional terrain surface data and the historical topographic survey results, and output the real-time topographic survey results, which may specifically include:

[0062] Step S510: Calculate the absolute value of the elevation deviation at the corresponding positions in the three-dimensional terrain surface data and the historical topographic survey results.

[0063] Step S510 is used to quantify the spatial difference between the current terrain data and the historical benchmark. The core lies in establishing an accurate elevation deviation assessment system. The absolute value of elevation deviation refers to the absolute value of the difference between the elevation value in the three-dimensional terrain surface data and the elevation value in the historical topographic survey results at the same spatial position. Its calculation needs to meet the spatio-temporal consistency constraints: 1) The spatial coordinate system and elevation benchmark are strictly unified; 2) The historical survey results corresponding to the time stamp should be comparable to the acquisition period of the current terrain data. The specific implementation process includes: First, extract the historical elevation data set of the target area from the historical topographic survey results. This data set contains the elevation point coordinates and corresponding elevation values collected at multiple historical time points, and the storage format is a spatio-temporal four-dimensional array (longitude, latitude, elevation, time); Second, perform spatial matching on the elevation point coordinates in the three-dimensional terrain surface data and the coordinates in the historical elevation data set. The matching method uses bilinear interpolation to resample the historical data to the resolution grid of the current terrain data; Subsequently, perform elevation difference calculation for each successfully matched coordinate pair. The formula is Δh = |h current(x,y) -h historical(x,y) |, where h current(x,y) is the current elevation value and h historical(x,y) is the historical elevation value; Finally, perform weighted average on the absolute value of the difference based on the terrain stability coefficient. The terrain stability coefficient is set according to the geological activity classification. For example, a lower weight of 0.3 is given in the landslide-prone area (geologically active area) to weaken the impact of instantaneous changes, and a higher weight of 0.7 is given in the plain area (terrain-stable area) to enhance the significance of long-term trends. The weighted average formula is Δh avg =Σ(w i ·Δh i ) / Σw i , where w i is the terrain stability coefficient of the area where the coordinate point i is located, and Δh i is the elevation difference. In addition, by introducing a terrain change sensitivity coefficient (calculated from the surface displacement rate, seismic activity frequency, and groundwater level change parameters), the weighted average value is non-linearly amplified. The final expression for the absolute value of elevation deviation is Δh abs =Δh avg ×k sense , where k sense is the sensitivity coefficient.

[0064] Step S520: Select the topographic map update mode according to the preset range where the absolute value of elevation deviation is located. Among them, when the absolute value of elevation deviation is less than the first threshold, the timed polling update mode is adopted; when the absolute value of elevation deviation is between the first threshold and the second threshold, the event-triggered update mode is adopted; when the absolute value of elevation deviation is greater than the second threshold, the real-time streaming update mode is started.

[0065] Step S520 aims to achieve dynamic switching of the topographic map update strategy to balance data freshness and computational resource consumption. The preset range is defined by the first threshold α1 and the second threshold α2, and their values are set according to the terrain change sensitivity classification of the target area. For example, in ordinary areas, α1 = 5 cm and α2 = 15 cm are set, and in disaster-sensitive areas, α1 = 2 cm and α2 = 8 cm are set. The timed polling update mode means full-scale updating of topographic data at fixed time intervals (such as 24 hours), which is applicable to scenarios where the terrain changes slowly or stably; the event-triggered update mode means starting incremental updates when the detected local elevation deviation exceeds α1 but is lower than α2, and only processing data in the changed area; the real-time streaming update mode means performing millisecond-level continuous processing and publishing of all-region data, which is applicable to scenarios such as geological disaster emergency response. In specific implementation, the update mode selection logic is as follows: First, compare the calculated absolute value of the elevation deviation Δh abs with α1 and α2; if Δh abs < α1, activate the timed polling update mode, and the system calls the preset timed task scheduler (such as Cron job) to trigger the data update process at a period T (such as T = 86400 seconds); if α1 ≤ Δh abs < α2, switch to the event-triggered update mode, and the system registers an elevation change event listener. When the Δh abs of a specific grid exceeds α1 three times continuously, generate an event message and trigger local data re-sampling and fusion; if Δh abs ≥ α2, immediately start the real-time streaming update mode, enable the high-frequency data output of the GNSS-RTK receiver (such as increasing from 1 Hz to 10 Hz), and allocate parallel computing resources to perform streaming processing.

[0066] Step S530: In the event-triggered update mode, dynamically adjust the data sampling interval according to the elevation deviation change rate. Among them, when the elevation deviation change rate exceeds the preset rate threshold, shorten the data sampling interval and increase the data output frequency of the GNSS-RTK receiver.

[0067] Step S530 is used to achieve adaptive optimization of the data acquisition frequency in the event-triggered update mode. The elevation deviation change rate is defined as the time derivative of the current absolute value of the elevation deviation relative to the previous measurement, and the calculation formula is v = Δh abs (t2) - Δh abs (t1) / (t2 - t1), and its unit is cm / hour. t2 is the current measurement time, and t1 is the previous measurement time. The preset rate threshold β is set according to the terrain change risk level. For example, in general monitoring scenarios, β = 0.5 cm / hour is set, and in landslide warning scenarios, β = 2 cm / hour is set. The dynamic adjustment strategy includes: when v ≤ β, maintain the basic sampling interval T base(e.g., 2 hours); when v > β, the sampling interval is shortened according to the rate overrun ratio, and the formula is T new =T base ×(β / v), and at the same time, the data output frequency of the GNSS-RTK receiver is increased from the conventional 1 Hz to ceil(v / β) × 1 Hz (ceil is the ceiling function). During the implementation process, it is necessary to synchronously adjust the data transmission bandwidth and the size of the storage buffer to ensure the stable processing of the high-frequency data stream. In addition, the system continuously monitors the rate change trend. If v drops below β and remains for more than 3 sampling periods, it will gradually return to the basic sampling parameters.

[0068] Step S540: In the real-time streaming update mode, the three-dimensional terrain surface data is divided into multiple data blocks, and each data block is independently verified and interpolated by a parallel computing node to generate processed data blocks.

[0069] Step S540 aims to achieve real-time distributed processing of large-scale terrain data. The data division adopts a spatial block strategy. The specific operation is as follows: First, according to the spatial range of the target area (such as the longitude span Δλ and the latitude span Δφ) and the number of computing nodes N, the three-dimensional terrain surface data is divided into N sub-blocks, and each sub-block contains a longitude range of [λ start +(i - 1)×Δλ / N, λ start +i×Δλ / N] and a latitude range of [φ start +(j - 1)×Δφ / M, φ start +j×Δφ / M] (i = 1, 2, … N; j = 1, 2, … M), where M×N is the total number of nodes, λ start is the starting value of longitude, φ startis the starting value of the latitude; secondly, a unique spatial identification code is assigned to each data block, and the coding rule is "lower longitude limit_upper longitude limit_lower latitude limit_upper latitude limit". For example, the identification code "116.300E_116.305E_39.900N_39.905N" represents the data block with longitude from 116.300° to 116.305° and latitude from 39.900° to 39.905°. The parallel computing nodes load a preset terrain data verification rule library, which includes elevation mutation detection rules (such as the single-point elevation change exceeding 3 times the standard deviation of the average value of adjacent points), contour continuity rules (such as the elevation difference between adjacent contour lines not exceeding the preset interval), and slope consistency rules (such as the deviation between the local slope and the regional average slope not exceeding 20%). The verification and processing process is as follows: each node uses a streaming data processing engine (such as Apache Flink) to scan each elevation point in the data block point by point. When an elevation mutation violation point is detected, it is marked as abnormal and the interpolation compensation mechanism is triggered. The interpolation method uses the sliding average value of the surrounding 8-neighborhood elevation points to calculate the compensation value; after the processing is completed, the node sends the verification log and the corrected data block to the central node.

[0070] Step S550: Perform spatial stitching and timestamp alignment on the processed data blocks, and output the real-time terrain mapping result with spatio-temporal consistency marks.

[0071] Step S550 is used to integrate the distributed processing results and ensure the consistency of the spatio-temporal attributes of the data. The spatial stitching operation includes: first, sort the processed data blocks according to the spatial identification codes of the data blocks, and arrange them in ascending order of longitude and descending order of latitude; secondly, use a seamless stitching algorithm to eliminate the boundary differences between blocks. The specific method is to extract the overlapping area with a width of 5 pixels at the boundary of adjacent data blocks and calculate the elevation mean difference Δh overlap , if Δh overlap ≤ the stitching threshold γ (such as γ = 0.1 m), directly stitch; if Δh overlap >γ, then perform Gaussian weighted average fusion on the overlapping area, and the weight is inversely proportional to the distance from the boundary. The timestamp alignment needs to ensure that the processing time deviation of all data blocks does not exceed the system clock accuracy (such as 1 ms). The implementation method is to record the start time t start and the end time t end in the metadata of the data block, and use t endAs the unified timestamp of this data block. The generation process of the spatio-temporal consistency mark is as follows: Combine the spatial identification code, timestamp, and processing pipeline number (such as "Node03_Phase2") of each data block into a globally unique identifier, for example, "116.300E-116.305E_39.900N-39.905N_20231012T153045Z_Node03_Phase2", and generate a fixed-length digest value (such as a 64-bit hexadecimal string) through the SHA-256 hash algorithm and embed it into the metadata segment of the data file. The format of the finally output real-time topographic mapping result is GeoPackage or LASer (LAS) point cloud file, which contains complete three-dimensional coordinates, elevation values, timestamps, and hash marks, and supports direct loading and integrity verification on the GIS platform.

[0072] As an implementation method, step S510: Calculate the absolute value of the elevation deviation at the corresponding position between the three-dimensional terrain surface data and the historical topographic mapping result, which may specifically include:

[0073] Step S511: Extract the historical elevation data set corresponding to the target area from the historical topographic mapping result. The historical elevation data set contains the elevation point coordinates and corresponding elevation values collected at multiple historical time points.

[0074] The spatio-temporal retrieval and standardization processing of the historical topographic data in step S511. The historical topographic mapping result is a verified topographic data set obtained in the past through lidar, photogrammetry, or GNSS-RTK technology, and its storage form is a spatio-temporal four-dimensional database (longitude, latitude, elevation, timestamp). Each data entry contains the elevation value of a unique coordinate point and its collection time. The extraction operation needs to execute the following process: First, perform a spatial range query in the historical database according to the boundary range of the target area (such as the coordinates of the vertices of the minimum bounding rectangle or the geographical fence polygon), and filter out all historical elevation points whose spatial positions are within the target area; Second, sort the filtered data by timestamp, and retain at least three data sets at different historical time points to support trend analysis, for example, select the survey results of three times with time spans of January 2020, June 2021, and March 2023; Subsequently, perform coordinate system normalization processing on the multi-period historical elevation data, and convert it to the same plane coordinate system (such as WGS84 UTM Zone 50N) and elevation datum (such as the EGM2008 geoid) as the current three-dimensional terrain surface data to eliminate systematic deviations caused by datum differences; Finally, construct a spatio-temporal index structure (such as an R-tree or a quadtree), and store the historical elevation data set in layers according to spatial positions and timestamps to support efficient spatial matching and time series analysis.

[0075] Step S512: Perform spatial matching on the elevation point coordinates in the three-dimensional terrain surface data and the coordinates in the historical elevation dataset to determine the successfully matched coordinate pairs.

[0076] Step S512 aims to achieve the precise spatial alignment of the current terrain data and the historical data. Spatial matching refers to, under the same coordinate system, determining whether the elevation points in the three-dimensional terrain surface data and the points in the historical elevation dataset are at the same spatial position through a coordinate tolerance threshold. The specific operation is as follows: First, calculate the distance between the elevation point coordinates (X current , Y current ) in the three-dimensional terrain surface data and the coordinates (X historical , Y historical ) in the historical elevation dataset point by point. The formula is D = √[(X current - X historical ) 2 +(Y current - Y historical ) 2 ; Second, set the spatial tolerance threshold δ (usually 1 / 2 of the data resolution, such as 0.25 meters). If D ≤ δ, it is determined that the match is successful, and a coordinate pair ((X current , Y current ),(X historical , Y historical )) is generated; if D > δ, it is marked as an unmatched point and excluded from subsequent calculations. For datasets with inconsistent resolutions, the low-resolution data needs to be resampled by bilinear interpolation to be consistent with the high-resolution data. For example, when the resolution of the current terrain data is 0.2 meters and a certain historical data is 1 meter, the historical data is encrypted to a 0.2-meter grid through interpolation to ensure the comparability of all coordinate points. The successfully matched coordinate pairs are stored in a list structure, and each entry contains the current elevation value h current , the historical elevation value h historical and the spatial tolerance D.

[0077] Step S513: For each successfully matched coordinate pair, calculate the absolute value of the difference between the current elevation value and the historical elevation value.

[0078] Step S513 is used to quantify the elevation change amplitude of a single point. The formula for the absolute value of the difference is Δh_ i = |h current _ i - h historical _ i |, where h current _ i is the elevation value of the i-th matched point in the three-dimensional terrain surface data, and h historical _ iis the elevation value corresponding to the historical elevation point. During the calculation process, if the historical elevation dataset contains multiple periods of data, it is necessary to select the historical elevation value with the closest time interval to the current data for calculation. For example, if the current data collection time is October 2023, the historical data of March 2023 should be given priority rather than the data of 2020; exclude outliers caused by sensor failures or environmental interferences (such as points where the elevation value exceeds the geographically reasonable range or the mutation exceeds 10 meters).

[0079] Step S514: Statistically calculate the weighted average of the absolute values of the differences of all successfully matched coordinate pairs. The weights for weighting are determined based on the terrain stability coefficients of the regions where the coordinate points are located. Among them, a first weight is set in geologically active regions, and a second weight is set in flat regions. The first weight is less than the second weight.

[0080] Step S514 eliminates the interference of the inherent characteristics of the terrain on the evaluation of the degree of difference through weighted averaging. The terrain stability coefficient is a preset weight value based on the geological activity classification of the target area. Its definition is as follows: Geologically active regions (such as fault zones and landslide-prone areas) are assigned a weight w1 (such as 0.3), flat regions (such as plains and stable sedimentation areas) are assigned a weight w2 (such as 0.7), and transitional regions (such as hills and gentle slopes) are assigned an intermediate weight w3 (such as 0.5). The calculation formula for the weighted average is Δh avg =(Σ(w i ×Δh i )) / Σw i , where w i is the terrain stability weight of the i-th coordinate point, and Δh i is the corresponding absolute value of the difference. The implementation process includes: First, query the classification label of the region where each matching point is located from the preset terrain stability zoning map and map it to the corresponding weight value; Second, traverse all matching points, accumulate the weighted absolute values of the differences and the total weight; Finally, calculate the weighted average and store it as a scalar value.

[0081] Step S515: Multiply the weighted average by the terrain change sensitivity coefficient to obtain the absolute value of the elevation deviation.

[0082] The ultimate goal of Step S515 is to amplify the elevation change signal in high-risk areas through the sensitivity coefficient. The terrain change sensitivity coefficient k sense is a dimensionless parameter, and its numerical range is usually from 1.0 to 3.0, which is obtained through the methods of Steps S5151 to S5156. The calculation formula for the absolute value of the elevation deviation Δh abs is Δh abs =Δh avg ×k sense . For example, if Δh avg =1.24 meters and k sense =1.8, then Δhabs = 1.24 × 1.8 = 2.23 meters. This coefficient can dynamically adjust the difference degree evaluation result according to the regional geological risk, ensuring that the small changes in sensitive areas (such as around volcanoes) are effectively amplified, while the large changes in stable areas are appropriately suppressed.

[0083] As an implementation method, the determination process of the terrain change sensitivity coefficient may include: Step S5151: Obtain the geological structure activity monitoring data of the target area, and the monitoring data includes the surface displacement rate, seismic activity frequency, and groundwater level change parameters.

[0084] The implementation process of Step S5151 involves the collection and preprocessing of multi-source geological monitoring data. The geological structure activity monitoring data includes: 1) Surface displacement rate, obtained through GNSS continuous operation reference station (CORS) or synthetic aperture radar interferometry (InSAR) technology, with the unit of millimeter / year, reflecting the crustal deformation intensity; 2) Seismic activity frequency, extract the historical earthquake event catalog of the target area from the seismic network database, and count the number of earthquakes with a magnitude ≥ 2.0 within a unit time (such as one year); 3) Groundwater level change parameters, calculate the water level change rate through the pressure sensor data of groundwater monitoring wells, with the unit of meter / month. Data preprocessing includes: performing deseasonalization trend processing on the surface displacement rate (eliminating the influence of environmental factors such as rainfall and temperature); performing sliding time window statistics on the seismic activity frequency (such as the number of earthquakes per month); performing Kalman filtering on the groundwater level data to remove noise.

[0085] Step S5152: Compare the surface displacement rate with the reference displacement rate to generate a displacement deviation index.

[0086] Step S5152 is used to quantify the abnormal degree of surface deformation. The reference displacement rate v base is a reference value set according to the regional geological background, usually taken from the long-term (such as 10 years) average deformation rate. The displacement deviation index I disp is calculated by the formula I disp =(v current - v base ) / v base , where v current is the current monitored displacement rate. For example, if in a certain area v base = 5 mm / year and v current = 12.5 mm / year, then I disp =(12.5 - 5) / 5 = 1.5. This index greater than 0 indicates accelerating deformation, and less than 0 indicates decelerating deformation.

[0087] Step S5153: Select a sensitivity adjustment factor according to the grade interval where the seismic activity frequency is located. The grade intervals include low seismic frequency interval, medium seismic frequency interval, and high seismic frequency interval.

[0088] The objective of step S5153 is to adjust the sensitivity according to the crustal activity intensity. The grade intervals are defined as follows: low earthquake frequency interval (annual earthquake times ≤ 2), sensitivity adjustment factor α1 = 1.0; medium earthquake frequency interval (3 ≤ annual earthquake times ≤ 5), α2 = 1.5; high earthquake frequency interval (annual earthquake times ≥ 6), α3 = 2.0.

[0089] Step S5154: Calculate the hydrological influence factor based on the groundwater level change parameter. The hydrological influence factor is positively correlated with the amplitude of water level change.

[0090] Step S5154 is used to quantify the influence of groundwater level change on surface stability. The hydrological influence factor I hydro has the calculation formula of I hydro =1 + 0.2×|Δh water |, where Δh water is the monthly average change amplitude of water level (unit: meter). For example, if Δh water = 0.2 meter, then I hydro = 1 + 0.2×0.2 = 1.04.

[0091] Step S5155: Normalize the displacement deviation index, sensitivity adjustment factor and hydrological influence factor to generate a comprehensive influence coefficient.

[0092] The operation of step S5155 is to eliminate the dimension difference and integrate the influence of multiple factors. The normalization formula is: I disp _ norm =(I disp -I min ) / (I max -I min ), where I min =-1.0, I max = 3.0; α norm =(α - 1.0) / (2.0 - 1.0), which is the normalized α value; I hydro _ norm =(I hydro - 1.0) / (2.0 - 1.0), which is the normalized I hydro value. The comprehensive influence coefficient C comb = 0.5×I disp _ norm + 0.3×α norm + 0.2×I hydro _ norm . For example, I disp = 1.5 → I disp _ norm =(1.5 + 1) / 4 = 0.625; α = 1.5 → α norm= 0.5; I hydro = 1.04 → I hydro _ norm = 0.04; then C comb = 0.5 × 0.625 + 0.3 × 0.5 + 0.2 × 0.04 = 0.3125 + 0.15 + 0.008 = 0.4705.

[0093] Step S5156: Perform a non - linear transformation on the comprehensive influence coefficient through an exponential function to generate a terrain change sensitivity coefficient.

[0094] Step S5156 maps the comprehensive influence coefficient to the sensitivity coefficient interval. The transformation formula is k sense = 1.0 + 2.0×(e (3×Ccomb) - 1) / (e 3 - 1), where the exponential function is used to amplify the changes in the high - influence area. For example, when C comb = 0.4705, k sense = 1.0 + 2.0×(e 1.4115 - 1) / (20.0855 - 1)= 1.0 + 2.0×(4.107 - 1) / 19.0855 ≈ 1.0 + 0.644 = 1.644. This coefficient will be used in the calculation of the absolute value of the elevation deviation in Step S515 to achieve dynamic amplification of terrain change sensitivity.

[0095] As an implementation, in Step S540, each data block is independently verified and interpolated through parallel computing nodes, which may specifically include:

[0096] Step S540: In the real - time streaming update mode, the three - dimensional terrain surface data is segmented into multiple data blocks, and each data block is independently verified and interpolated through parallel computing nodes to generate processed data blocks.

[0097] Step S540 is used to achieve efficient real-time processing of large-scale terrain data through a distributed computing architecture. The data block is a subset of the three-dimensional terrain surface data divided by spatial range, and its division rule is based on the longitude span, latitude span and elevation distribution characteristics of the target area. In specific implementation, firstly, the maximum allowable size of a single data block is determined according to the number of parallel computing nodes and hardware resources (such as GPU memory capacity). For example, each data block is set to cover a geographical range of 0.1°×0.1° (approximately 11 kilometers×11 kilometers) and an elevation range of -100 meters to 9000 meters; secondly, a spatial grid division algorithm is used to cut the original three-dimensional terrain surface data into multiple non-overlapping regular cube data blocks. During the cutting process, ensure that the boundaries of the data blocks are aligned with the geographic coordinate grid lines to avoid cross-blocks. Data redundancy; finally, a unique spatial identification code is generated for each data block. The identification code is composed of a six-tuple of longitude lower limit, longitude upper limit, latitude lower limit, latitude upper limit, elevation lower limit, and elevation upper limit. For example, the identification code "E116.300-E116.400_N39.900-N40.000_H-100-H9000" represents a data block with longitude 116.300° to 116.400°, latitude 39.900° to 40.000°, and elevation -100 meters to 9000 meters. After receiving the data block, the parallel computing node loads the preset terrain data verification rule base, which contains elevation mutation detection rules (such as the difference between the elevation of a single point and the average value of the surrounding 8 neighborhoods exceeds 3 times the standard deviation), contour line continuity rules (such as the spacing between adjacent contour lines must not exceed 1.5 times the preset threshold) and slope consistency rules (such as the absolute value of the deviation between the local slope and the regional average slope is less than 20%). In the verification and interpolation processing flow, the streaming data processing engine (such as Apache Kafka Streams) scans the data block point by point in the order of the timestamps of the elevation points. When a certain elevation point is detected to violate the elevation mutation detection rule, it is immediately marked as an abnormal point and the interpolation compensation mechanism is triggered. For example, if the elevation value of a certain point is 152.3 meters and the average value of its 8 neighborhood points is 148.1 meters with a standard deviation of 1.2 meters, the 3 times standard deviation threshold is 148.1±3.6 meters, and the point is marked as abnormal because it exceeds the upper limit. Interpolation compensation uses the surrounding neighborhood sliding average algorithm. For example, the neighborhood elevation points of a 5×5 window are taken with the abnormal point as the center, and the arithmetic mean of the remaining points is calculated as the compensation value after the highest and lowest 10% of the outliers are removed. After the processing is completed, all abnormal points in the data block are corrected, and the verification log records the abnormal type, location and correction value. The corrected data block is sent to the central node with an additional spatial identification code and processing timestamp.

[0098] Step S541: Assign a unique spatial identification code to each data block, where the spatial identification code includes information on longitude range, latitude range and elevation range.

[0099] Step S541 aims to establish a spatial metadata identification system for data blocks to support distributed processing and data traceability. The spatial identification code is a string code that conforms to the ISO 19115 Geographical Information - Metadata standard. Its generation rule is as follows: Concatenate the lower and upper limits of longitude, lower and upper limits of latitude, and lower and upper limits of elevation of the data block in a fixed format, separated by underscores between each field. The numerical precision is reserved to three decimal places. Longitude and latitude are expressed in decimal degrees, and elevation is in meters. For example, for a data block covering longitude from 116.300° to 116.400°, latitude from 39.900° to 40.000°, and elevation from - 100 meters to 500 meters, its spatial identification code is "E116.300 - E116.400_N39.900 - N40.000_H - 100 - H500". In specific implementation, the allocation operation needs to be executed synchronously with the data segmentation process: When the three - dimensional terrain surface data is cut into cube - shaped data blocks, the boundary coordinates of each block are calculated in real - time and the corresponding spatial identification code is generated. This code is embedded in the metadata header of the data block and simultaneously registered in the spatial index service of the central node. The functions of the spatial identification code include: 1) Quickly locate the geographical range of the data block in parallel computing nodes; 2) Guide the sorting and spatial relationship reconstruction of data blocks during the data splicing stage; 3) Support quickly associating the original data location during abnormal data traceability.

[0100] Step S542: Load a preset terrain data verification rule library in each parallel computing node. The rule library includes elevation mutation detection rules, contour line continuity rules, and slope consistency rules.

[0101] Step S542 is used to deploy a unified data quality control standard in distributed computing nodes. The terrain data verification rule library is a configuration file in XML or JSON format, and its content is defined as follows: 1) The elevation mutation detection rule stipulates that for any elevation point P(x, y), if the absolute difference between its elevation value h(x, y) and the average elevation value μ8 of its surrounding 8 - neighborhood points exceeds 3 times the neighborhood standard deviation σ8, that is, |h(x, y)-μ8|>3σ8, then P is determined to be an abnormal point; 2) The contour line continuity rule requires that the elevation difference Δh contour between adjacent contour lines shall not exceed 1.5 times the contour interval D, that is, Δh contour ≤1.5D. If it is detected that a certain section of the contour line violates this rule, then mark this section as a fracture area; 3) The slope consistency rule is defined as that the absolute value of the deviation between the local slope value S(x, y) and the regional average slope S avg shall not exceed 20%, that is, |S(x, y)-S avg | / S avg≤0.2. The rule library loading process includes: when the parallel computing node starts, downloading the latest version of the rule configuration file from the central storage system (such as HDFS or Amazon S3), parsing it into rule objects in memory, and allocating independent data detection threads for each rule. For example, the elevation mutation detection thread scans the input data stream in real time, calls the neighborhood statistical function to calculate μ8 and σ8, and the slope consistency thread derives the slope raster from the DEM and compares it with the regional average slope.

[0102] Step S543: Use a streaming data processing engine to verify each elevation point in the data block one by one. When any elevation point is detected to violate the elevation mutation detection rule, mark any elevation point as an abnormal point and trigger the interpolation compensation mechanism.

[0103] The specific implementation of step S543 can utilize the high-throughput real-time computing capabilities of a streaming data processing engine (such as Apache Flink or SparkStreaming). The streaming verification process is as follows: the data block is split into an elevation point stream arranged in row-major order, and each point contains longitude, latitude, elevation value, and timestamp; the engine reads each point in sequence, retrieves its surrounding neighborhood points (such as a 5×5 window) according to the spatial coordinates, and calculates the neighborhood statistics (mean, standard deviation); subsequently, substitute the elevation value and statistics of the current point into the elevation mutation detection rule for determination. If the rule is violated, insert an abnormal flag bit (such as setting the abnormal flag to 1) into the metadata of this point and trigger the interpolation compensation event. After the event is triggered, the engine pauses the processing of the current data stream, calls the interpolation compensation module to generate a compensation value, and writes the corrected elevation value back into the data stream.

[0104] Step S544: For the positions marked as abnormal points, use the moving average of the surrounding elevation points for interpolation compensation to generate the compensated elevation points.

[0105] The interpolation compensation algorithm for step S544 can adopt a robust statistical method to reduce the influence of outliers. The specific operation can be, for example: taking the abnormal point as the center, selecting a circular neighborhood with a radius R (such as R = 5 pixels), and extracting all elevation points within this area; sorting the neighborhood points by elevation value, and removing the highest 10% and the lowest 10% of the extreme values; calculating the arithmetic mean of the remaining points as the compensation value, and the formula is h comp =Σh i / (N - 2k), where N is the total number of neighborhood points, k = floor(0.1×N), h comp is the compensated elevation value, h iis the elevation value of each point in the neighborhood. For example, the 5×5 neighborhood of an abnormal point contains 25 points. After removing 2 highest values (such as 155.6 m, 154.9 m) and 2 lowest values (such as 145.2 m, 146.0 m), the average value of the remaining 21 points is 148.3 m, which is used as the compensation value to replace the original abnormal value of 152.3 m. The compensated elevation points retain the original coordinates and timestamps, and record the type of compensation operation (such as "moving average compensation_R5") and the number of neighborhood points involved in the calculation in the metadata.

[0106] Step S545: Send the verified and compensated data blocks to the central node so that the central node can sort and splice the data blocks according to the spatial identification code.

[0107] The data transmission protocol in Step S545 needs to ensure the integrity and temporal consistency of the processed data blocks. The specific implementation includes: after each parallel computing node finishes processing the data block, it encapsulates the data block into a data packet according to the lexicographical order of the spatial identification code (ascending longitude → ascending latitude → ascending elevation) and sends it to the central node through the TCP protocol; after receiving the data packet, the central node parses the spatial identification code and registers it in the global spatial index table. The index table is organized in a quadtree structure to support fast spatial queries. During the sorting process, the central node arranges the data blocks into a two-dimensional grid sequence according to the longitude and latitude ranges of the spatial identification code. For example, the identification code "E116.300-E116.400_N39.900-N40.000" is arranged after "E116.200-E116.300_N39.900-N40.000", forming a splicing order from west to east and from south to north. The data blocks are cached in memory in this order, waiting for the splicing instruction.

[0108] Step S546: During the splicing process, detect the elevation difference at the boundary of adjacent data blocks. If the elevation difference exceeds the preset splicing threshold, start the boundary smoothing algorithm to re-interpolate the boundary area.

[0109] Step S546 is used to eliminate the block indirect seam problem caused by distributed processing. The preset splicing threshold γ is set according to the data accuracy requirements. For example, γ = 0.1 m means that the elevation difference between corresponding points at the boundary of adjacent blocks shall not exceed 0.1 m. The detection process is as follows: extract the overlapping area with a width of 2 pixels at the boundary of adjacent data blocks (such as the east boundary of data block A and the west boundary of data block B), and calculate the absolute value of the elevation difference Δh = |h A(x,y) -h B(x,y) |; count the proportion of the number of points where Δh exceeds γ in the total number of points. If the proportion exceeds 5%, it is determined as a serious seam problem and the boundary smoothing algorithm is triggered. For example, the overlapping area of adjacent blocks A and B contains 100 points, and among them, 8 points have Δh > 0.1 m, with a proportion of 8% > 5%, so smoothing processing needs to be performed.

[0110] As an implementation, in step S546, starting the boundary smoothing algorithm to re- interpolate the boundary region may include:

[0111] Step S5461: Extract the elevation points within a preset width range at the boundary of adjacent data blocks to generate a boundary elevation point set.

[0112] The preset width in step S5461 can be set to, for example, 2 to 5 pixel widths of the data block boundary. For example, for 0.5-meter resolution data, the preset width is 2 pixels, which is 1.0 meter. The extraction operation includes: Extracting all elevation points in the column where the maximum longitude value is located and the previous column (a total of 2 columns) from the east boundary of data block A, and at the same time extracting all points in the column where the minimum longitude value is located and the next column (a total of 2 columns) from the west boundary of data block B, and combining them to generate a boundary elevation point set. For example, if the longitude of the east boundary of data block A is 116.400°, then extract all points with longitudes from 116.395° to 116.400° (2 pixel widths); if the longitude of the west boundary of data block B is 116.400°, extract points from 116.400° to 116.405°. Each point in the set records the source block identifier (A or B) and the original elevation value.

[0113] Step S5462: Calculate the average elevation value and elevation standard deviation of all points in the boundary elevation point set.

[0114] The statistical calculation in step S5462 aims to quantify the elevation consistency of the boundary region. For example, the average elevation value μ boundary =Σh i / N, where N is the total number of boundary points; the elevation standard deviation σ boundary =√[Σ(h i -μ boundary ) 2 / (N - 1)], where h i refers to the elevation value of each point in the boundary elevation point set. For example, the boundary set contains 200 points, the average elevation μ = 150.2 meters, and the standard deviation σ = 0.3 meters. These statistical values are used for subsequent Gaussian model construction and confidence weight calculation.

[0115] Step S5463: Based on the average elevation value, construct a Gaussian distribution model to generate the elevation confidence weight of each elevation point.

[0116] For example, the Gaussian distribution model is defined as the probability density function P(h)=exp[-(h - μ) 2 / (2σ 2 )] / (σ√(2π)), and the confidence weight w i =P(h i ) / P max , where P maxis the distribution peak (P(μ)=1 / (σ√(2π))), h is the elevation value, and σ is the standard deviation. The weights are normalized to the range of 0 to 1, and the points closer to μ have higher weights.

[0117] Step S5464: Weightedly fuse the elevation points in the boundary area according to the elevation confidence weights to generate the fused boundary elevation value.

[0118] For example, the weighted fusion formula is h fused =Σ(w Ai ×h Ai +w Bi ×h Bi ) / Σ(w Ai +w Bi ), where w Ai and w Bi are the corresponding point weights from data blocks A and B respectively, h Ai is the elevation value of data block A, and h Bi is the elevation value of data block B.

[0119] Step S5465: Use bilinear interpolation to perform a smooth transition on the fused boundary elevation value to eliminate the stitching gap.

[0120] The bilinear interpolation in step S5465 constructs a continuous elevation transition surface within the boundary area. The specific operation is as follows: Take the fused boundary elevation points as control points, generate an interpolation grid in the east-west direction (longitude) and north-south direction (latitude). For any point to be interpolated (x, y), its elevation h(x,y)=a×x + b×y + c×x×y + d, and the coefficients a, b, c, d are solved by least squares fitting of the control points. For example, generate an interpolation grid with a resolution of 0.2 meters within a 1.0-meter-wide boundary area to smoothly transition the elevation value of data block A from 150.4 meters to 150.1 meters of data block B.

[0121] Step S5466: Update the processed boundary elevation value to the adjacent data blocks and recalculate the elevation difference in the stitching area.

[0122] For example, write the interpolated boundary elevation values back to the corresponding positions of data blocks A and B respectively, overwriting the original boundary data; then, re-extract the elevation points in the boundary area and calculate Δh new =|h_A new (x,y)-h_B new (x,y)|, where h_A new (x,y) and h_B new (x,y) respectively refer to the elevation values of the updated data blocks A and B at the coordinate (x, y), and verify whether it meets the γ threshold. For example, the original Δh max= 0.3 m, after smoothing, Δh max = 0.05 m < γ = 0.1 m, the splicing gap is eliminated. The updated data block is resent to the central node to participate in the construction of the global terrain model.

[0123] As an implementation manner, the generation process of the spatio-temporal consistency mark may specifically include:

[0124] Step S551: Attach a data acquisition timestamp and a spatial coordinate range identifier to each verified data block.

[0125] Step S551 is used to assign accurate spatio-temporal attribute identifiers to the terrain data blocks generated by distributed processing. The data acquisition timestamp is a Coordinated Universal Time (UTC) time code following the ISO 8601 standard, and the accuracy is not limited. For example, it is at the millisecond level, and the format can be set to "YYYYMMDDThhmmss.sssZ". For example, "20231015T083045.123Z" represents 08:30:45.123 seconds on October 15, 2023. The spatial coordinate range identifier is generated according to the geographical coverage range of the data block, including a six-tuple of lower longitude, upper longitude, lower latitude, upper latitude, lower elevation, and upper elevation. The numerical accuracy is reserved to three decimal places. Longitude and latitude are represented in decimal degrees, and elevation is in meters. For example, the identifier "E116.300 - E116.400_N39.900 - N40.000_H - 100 - H500" represents a data block with a longitude range of 116.300° to 116.400°, a latitude range of 39.900° to 40.000°, and an elevation range of -100 m to 500 m. In a specific implementation, the timestamp is synchronously generated by the Global Positioning System (GPS) timing module of the data acquisition device, and the spatial coordinate range identifier is automatically calculated by a spatial grid division algorithm during the data block segmentation stage. The attachment operation is implemented by modifying the metadata segment of the data block: the timestamp is written into the "AcquisitionTime" field, and the spatial coordinate range identifier is written into the "Spatial Extent" field, both of which use binary encoding to improve storage efficiency.

[0126] Step S552: Verify the topological relationship of the spatial coordinate range identifiers of all data blocks with the same timestamp until there are no overlapping or missing areas.

[0127] Step S552 aims to ensure the spatial integrity and seamlessness of data blocks at the same acquisition moment. For example, the topological relationship verification includes two checks: 1) no omission in spatial coverage, that is, the union of the spatial coordinate range identifiers of all data blocks needs to completely cover the target area; 2) no spatial overlap, that is, the intersection of the spatial coordinate range identifiers of any two data blocks is empty. During implementation, first, extract the spatial coordinate range identifiers of all data blocks with the same timestamp from the spatio-temporal index library of the central node and arrange them in ascending order of longitude and ascending order of latitude; subsequently, use a computational geometry library (such as GEOS) to perform topological relationship analysis on the boundaries of adjacent data blocks: by calculating whether the upper longitude limit of data block A is equal to the lower longitude limit of data block B and whether the latitude ranges are continuous, determine whether there are gaps or overlaps. If a gap is detected (for example, the upper longitude limit of data block A is 116.400° while the lower longitude limit of data block B is 116.405°), generate a regional omission warning and trigger the data re-acquisition process; if an overlap is detected (for example, the upper longitude limit of data block A is 116.400° and the lower longitude limit of data block B is 116.395°), mark it as a conflict area and initiate the re-segmentation of the data block.

[0128] Step S553: Generate a processing pipeline number according to the processing order of the data block, and the number includes the node identifier and the data processing stage code.

[0129] The goal of Step S553 is to record the processing path of data blocks in distributed computing to support end-to-end traceability. The encoding rule of the processing pipeline number is: the node identifier uses the globally unique node ID assigned by the distributed computing framework (such as "Node07"), and the data processing stage code consists of the English abbreviation of the processing link and the serial number (such as "VAL1" represents the first verification and "INT2" represents the second interpolation). The number generation logic is: when a data block enters a certain computing node, the node scheduler appends the current processing stage code to its metadata and records the node ID and timestamp. For example, if a data block passes through the verification stage (VAL1) of node Node03 and the interpolation stage (INT2) of node Node12 in sequence, the pipeline number sequence is "Node03_VAL1_20231015T083045.123Z → Node12_INT2_20231015T083102.456Z". The number is stored in the "PipelineID" field of the data block metadata and uses the JSON array format to record the multi-stage processing history. This numbering mechanism enables the location of the problem link by backtracking the pipeline number in case of data anomalies. For example, if a data block is marked as abnormal at the VAL1 stage of node Node05, the logs of this node can be quickly retrieved to analyze the root cause.

[0130] Step S554: Combine and encode the timestamp, spatial coordinate range identifier, and processing pipeline number to generate a globally unique spatio-temporal consistency tag.

[0131] The implementation process of Step S554 ensures the global uniqueness and information integrity of the tag through structured encoding. The combination encoding rule is as follows: the timestamp, spatial coordinate range identifier, and processing pipeline number are concatenated in a fixed order, separated by "#" between each field. For example, "20231015T083045.123Z#E116.300-E116.400_N39.900-N40.000_H-100-H500#Node03_VAL1→Node12_INT2". The encoded string is converted to a compact format through the Base64 algorithm. For example, the original string is converted to "MjAyMzEwMTVUMDgzMDQ1LjEyM1ojRTExNi4zMDAtRTExNi40MDBfTjM5LjkwMC1OMzkuOTAwX0gtMTAwLUg1MDAjTm9kZTAzX1ZBTDEtPk5vZGUxMl9JTlQy". The global uniqueness is guaranteed by the following mechanisms: 1) The timestamp is accurate to the millisecond level to avoid conflicts between different batches of data; 2) The spatial coordinate range identifier covers a unique geographical area; 3) The processing pipeline number includes the node ID and timestamp to ensure differences in the processing path. For example, data blocks collected in the same geographical area at different times will generate different tags due to different timestamps; adjacent data blocks at the same time will also have different tags due to different spatial ranges.

[0132] Step S555: During the process of outputting the real-time terrain mapping result, embed the spatio-temporal consistency tag into the metadata segment of the data file and perform a hash binding with the three-dimensional terrain surface data.

[0133] Step S555 is used to ensure data integrity and anti-tampering through digital signature technology. The embedding operation includes: writing the spatio-temporal consistency tag into the metadata segment of the standard geographic data format (such as GeoTIFF or LAS), specifically at the "CustomTags" field in the file header. The hash binding process is as follows: First, use the SHA-256 algorithm to calculate the hash of all elevation points of the three-dimensional terrain surface data, generating a 64-bit hexadecimal digest value (such as "a1b2c3d4e5f6..."); Subsequently, encrypt the digest value and merge it with the spatio-temporal consistency tag to generate a digital signature string; Finally, write the signature into the "DigitalSignature" field of the data file. For example, the metadata of a certain data file contains the fields "CustomTags=Base64_MjAyMzEwMTVUMDgzMDQ1LjEyM1ojRTExNi4zMDAtRTExNi40MDBfTjM5LjkwMC1OMzkuOTAwX0gtMTAwLUg1MDAjTm9kZTAzX1ZBTDEtPk5vZGUxMl9JTlQy" and "DigitalSignature=sha256:a1b2c3d4e5f6...". During data verification, the recipient can confirm that the data has not been tampered with by recalculating the hash value and comparing it with the signature. If a certain data block is maliciously modified during transmission (such as changing the elevation value of a certain point from 150.3 meters to 160.3 meters), its hash value will change significantly, resulting in a failed signature verification and triggering a security alert.

[0134] In an optional derivative implementation manner, after outputting the real-time terrain mapping result in step S500, the method provided by the embodiments of the present invention may further include:

[0135] Step S600: Detect abnormal terrain areas in the real-time terrain mapping result, and generate boundary coordinates of the abnormal area and an abnormal type identifier.

[0136] Step S600 is used to identify abnormal areas in the real-time terrain data that do not conform to the natural landform evolution law. The real-time terrain mapping result is the three-dimensional terrain surface data with spatio-temporal consistency tags output in step S500, and its format is a gridded elevation point set or a triangulation network model. The detection of abnormal terrain areas is achieved through multi-dimensional terrain feature analysis, specifically including the calculation of the elevation gradient change rate, the assessment of historical terrain stability, and the modeling of abnormal probability. During implementation, first extract the elevation gradient change rate from the real-time terrain data. This parameter is defined as the change in elevation value within a unit horizontal distance, and the calculation formula is G = √[(∂h / ∂x) 2 +(∂h / ∂y) 2 , where ∂h / ∂x and ∂h / ∂y are the elevation gradient components in the east-west and north-south directions respectively. The preset threshold Gthreshold Dynamically set according to the regional terrain characteristics. For example, it is set to 0.5 m / m in plain areas and 2.0 m / m in mountainous areas. When the elevation gradient change rate of a certain continuous area exceeds G threshold , it is marked as a candidate abnormal area. For example, in a certain area, G = 2.5 m / m is detected within a 5×5 pixel range, exceeding the threshold for mountainous areas and being listed as a candidate anomaly. Subsequently, obtain the historical terrain stability index S index of the candidate area, and its calculation formula is S index = 1 / (σ 2 × e (λΔt) ), where σ 2 is the variance value of historical elevation changes, Δt is the interval between the current time and the historical data collection time, and λ is the time decay coefficient (e.g., λ = 0.01 / day). The weighted comparison process combines the elevation gradient change rate G and S index to generate the regional anomaly probability value P anomaly = α × G / G max + (1 - α) × (1 - S index ), where α is the weight factor (e.g., α = 0.6), and G max is the maximum gradient value of the area. When P anomaly exceeds the dynamic threshold P threshold (e.g., 0.7), it is determined as a valid abnormal area, and an anomaly type identifier is generated by comparing the gradient change pattern with historical data: natural erosion identifier (gradual gradient change and matching the historical erosion pattern), artificial construction identifier (sudden gradient change and matching the building outline), and geological collapse identifier (gradient distributed in concentric circles and accompanied by negative elevation changes).

[0137] Step S700: Obtain the multi-spectral remote sensing image data corresponding to the boundary coordinates of the abnormal area, and extract the surface reflectance characteristics and vegetation cover density parameters within the abnormal area.

[0138] Step S700 is used to enhance the surface property analysis of the abnormal area through multi-spectral remote sensing data. The multi-spectral remote sensing image data is a high-resolution image containing bands such as visible light, near-infrared, and short-wave infrared (e.g., Sentinel-2 or Landsat-8 data), and its spatial resolution needs to match the boundary coordinates of the abnormal area (e.g., 10 m). The extraction of surface reflectance characteristics includes: 1) Calculate the normalized difference vegetation index (NDVI) = (NIR - Red) / (NIR + Red), where NIR is the reflectance of the near-infrared band and Red is the reflectance of the red band; 2) Calculate the modified soil-adjusted vegetation index (MSAVI) = (2×NIR + 1 - √((2×NIR + 1) 2-8×(NIR - Red)) / 2, which is used to reduce soil background interference; 3) Calculate the Normalized Difference Water Index (NDWI) = (Green - NIR) / (Green + NIR) to identify the water body distribution. The vegetation coverage density parameter is calculated by the pixel decomposition model, which decomposes the mixed pixel into the proportions of endmembers such as vegetation, bare soil, and water body. The vegetation coverage density V cover = Proportion of vegetation endmember × 100%.

[0139] Step S800: Match the surface reflectance characteristics with the preset geological disaster characteristic library. When the match is successful, activate the terrain review mode, send a high - density sampling instruction to the GNSS - RTK receiver, and trigger the collaborative mapping of the airborne lidar.

[0140] Step S800 aims to identify geological disaster risks through multi - spectral feature matching and initiate the review mechanism. The geological disaster characteristic library is a database containing the spectral characteristics of typical disasters (landslide, debris flow, collapse), and its storage form is a collection of multi - dimensional vectors. Each vector contains historical statistical values of parameters such as NDVI, MSAVI, NDWI, V cover and so on. The matching algorithm uses cosine similarity calculation: Sim = Σ(F i × F' i ) / √(ΣF i 2 × ΣF' i 2 ), where F i is the current surface reflectance characteristic, and F' i is the reference vector in the characteristic library. When Sim ≥ 0.85, it is determined that the match is successful and the terrain review mode is activated. For example, the similarity Sim between the current feature vector and the "landslide" feature in the library is 0.89, triggering the review process. The high - density sampling instruction increases the data acquisition frequency of the GNSS - RTK receiver from 1 Hz to 10 Hz. At the same time, the airborne lidar starts scanning, and the scanning parameters are set as: pulse frequency 500 kHz, scanning angle ±30°, and point density ≥ 50 points per square meter.

[0141] Step S900: Based on the spatial registration result of the high - density sampling data and the lidar point cloud data, correct the position of the abnormal area boundary coordinates, generate the corrected terrain mapping result, and update it to the three - dimensional terrain surface data.

[0142] Step S900 is used to improve the positioning accuracy of the abnormal boundary through multi-source data fusion. For spatial registration, the Iterative Closest Point (ICP) algorithm is adopted to perform rigid transformation (translation and rotation) on the GNSS-RTK high-density sampled point cloud and the lidar point cloud to minimize the point distance error. After registration, a fused positioning reference plane is generated, and its coordinate system is unified to WGS84 UTM Zone 50N, and the elevation datum is EGM2008. The double-constraint adjustment calculation introduces the positioning accuracy model of GNSS-RTK (horizontal ±1 cm, elevation ±2 cm) and the scanning error model of lidar (angle error ±0.01°) to construct a joint adjustment equation: Σ(w gnss ×(x gnss -x adj ) 2 +w lidar ×(x lidar -T(x adj )) 2 )→min, where T is the coordinate transformation matrix, w gnss =1 / σ gnss 2 , w lidar =1 / σ lidar 2 , x gnss refers to the coordinate value measured by GNSS (Global Navigation Satellite System), x adj is the adjusted coordinate value, w gnss , w lidar are the weights calculated according to their respective error models. For example, the horizontal error of a certain boundary point is reduced from ±5 cm to ±1.5 cm after adjustment. For morphological closing operation, a 3×3 circular structural element is used to perform dilation-erosion operation on the boundary point set to fill pores with a diameter less than 2 meters. For slope continuity transition processing, anisotropic diffusion filtering is applied to smooth the gradient along the boundary normal direction to ensure that the corrected boundary line is naturally connected to the surrounding terrain. Finally, the updated three-dimensional terrain surface data replaces the elevation values of the original abnormal area in an incremental manner, and a terrain classification label (such as "Geological Collapse Area_L2") is added.

[0143] As an implementation manner, step S600, for detecting the terrain abnormal area from the real-time terrain mapping result and generating the boundary coordinates and abnormal type identification of the abnormal area, may include:

[0144] Step S610: Extract continuous areas with the elevation gradient change rate exceeding a preset threshold from the real-time terrain mapping result and mark them as candidate abnormal areas.

[0145] For example, perform Sobel operator convolution calculation on the real-time terrain data to calculate the elevation gradient, and set the dynamic threshold G threshold =μ G +3σ G, where μ G is the regional average gradient, and σ G is the standard deviation. For example, in a certain region, μ G = 0.8 m / m, and σ G = 0.2 m / m, then G threshold = 0.8 + 3×0.2 = 1.4 m / m. A continuous region is defined as a connected domain with at least 5×5 pixels and a gradient exceeding the limit. The candidate abnormal regions are marked through the region growing algorithm. The coordinates of the candidate regions are stored in the polygon vector format, and the attribute table records the gradient mean, area, and centroid position.

[0146] Step S620: Obtain the historical terrain stability index of the candidate abnormal regions. The stability index is obtained by calculating the variance value of the elevation change in the historical surveying and mapping data and the time decay coefficient.

[0147] In step S620, the historical elevation change variance σ 2 is obtained by calculating the difference between multiple periods of DEMs. For example, select 3 periods of DEM data in the recent 5 years, calculate the elevation change amount Δh t , and the variance σ 2 = Σ(Δh t - μ Δh ) 2 / (n - 1), where μ Δh refers to the average value of the historical elevation change amounts. The time decay coefficient λ = 0.01 / day. For example, in a certain region, σ 2 = 0.25 m 2 , and Δt = 365 days, then S index = 1 / (0.25×e (0.01×365) ) = 1 / (0.25×38.5) = 0.104, indicating low stability.

[0148] Step S630: Compare the elevation gradient change rate and the historical terrain stability index with weights to generate a regional abnormal probability value.

[0149] For example, the weighting formula can be P anomaly = 0.6×(G / G max ) + 0.4×(1 - S index ), where G max = 5.0 m / m. For example, G = 2.5 m / m, and S index = 0.104, then P anomaly = 0.6×(2.5 / 5.0) + 0.4×(1 - 0.104) = 0.3 + 0.358 = 0.658. When it is lower than the threshold of 0.7, no alarm is triggered.

[0150] Step S640: When the regional abnormal probability value exceeds the dynamic threshold, it is determined as a valid abnormal region and an abnormal type identifier is generated.

[0151] Dynamic threshold P threshold = 0.7. When P anomaly ≥ 0.7, generate an identifier by combining the gradient direction distribution and the spectral feature matching result. For example, for a certain area where P anomaly = 0.82, the gradient is radially distributed, and the NDVI decreases by 0.2, it is determined as a geological collapse identifier.

[0152] Step S650: Determine the corresponding verification strategy according to the abnormal type identifier.

[0153] For example, the verification strategy may include: Invoking historical multi - period satellite images (such as at 5 - year intervals) for change detection in the natural erosion identifier area; Querying the local government construction permit database in the artificial construction identifier area to match the coordinates with the permit scope; Starting a ground - penetrating radar (GPR) or resistivity imaging equipment in the geological collapse identifier area to scan for underground cavities.

[0154] As an implementation, in step S900, based on the spatial registration result of the high - density sampling data and the lidar point cloud data, correct the position of the abnormal area boundary coordinates to generate a corrected topographic mapping result, which may include:

[0155] Step S910: Unify the spatial coordinate systems of the GNSS - RTK elevation points in the high - density sampling data and the lidar point cloud data to generate a fused positioning reference plane.

[0156] For example, the coordinate system unification uses the seven - parameter Helmert transformation to convert the lidar point cloud from the scanning coordinate system to the WGS84 frame of GNSS - RTK, with translation parameters ΔX = 1.2m, ΔY = - 0.8m, ΔZ = 0.5m, rotation angles ω = 0.01°, φ = 0.005°, κ = 0.003°, and scale factor s = 1.000015.

[0157] Step S920: Perform double - constraint adjustment calculations on the abnormal area boundary coordinates on the fused positioning reference plane to eliminate the GNSS - RTK signal jitter error and the lidar scanning angle error.

[0158] For example, the adjustment model can be: minΣ( (x gnss -x) 2 / σ gnss 2 +(x lidar -T(x)) 2 / σ lidar 2 ), where σ gnss = 0.01m, σ lidar= 0.02 m, where T is the coordinate transformation matrix. After adjustment, the horizontal accuracy of the boundary points is improved to ±0.015 m. x gnss is the coordinate value observed by the GNSS-RTK receiver (plane coordinates x, y or three-dimensional coordinates x, y, z); x is the "true" coordinate value to be solved (the optimal estimated value after adjustment); x lidar is the original point cloud coordinate value obtained by lidar scanning, and T(x) is the coordinate transformation function that transforms the lidar point cloud from the local coordinate system to the global coordinate system of GNSS-RTK; σ gnss and σ lidar are the standard deviations of the observations of GNSS-RTK and lidar, representing the data accuracy.

[0159] Step S930: Extract the adjusted boundary point set and perform morphological closing operation to fill the boundary pores caused by data loss.

[0160] For example, the morphological closing operation uses a 3×3 circular structuring element, first dilating (filling pores) and then eroding (restoring the boundary shape). For example, a boundary pore with a diameter of 1.5 meters is completely filled after the closing operation.

[0161] Step S940: Perform slope continuity transition processing on the processed boundary point set and the terrain data of the surrounding normal areas to generate a smoothly transitioned corrected boundary line.

[0162] For example, the slope transition algorithm is as follows: within a 10-meter buffer zone on both sides of the boundary line, the slope values are interpolated by distance weighting, and the formula is set as S(x,y) = w b × S b + w n × S n where w b = 1 - d / D, w n = d / D, d is the distance from the boundary, D = 10 m. S(x,y) is the mixed slope value; S b is the slope value of the boundary area; S n is the slope value of the normal area; w b and w n respectively represent the mixing ratios of the boundary slope (Sb) and the normal area slope (Sn); d is the vertical distance from the current point (x,y) to the boundary line; D is the total width of the buffer zone.

[0163] Step S950: Re-divide the terrain partition based on the corrected boundary line, and update the elevation values and terrain classification labels of the corresponding areas in the three-dimensional terrain surface data.

[0164] For example, the partition update operation includes: converting the correction boundary line into a raster mask, replacing the elevation values of the original abnormal area, and adding classification labels to the attribute table. For example, the original label "unclassified" is updated to "geological collapse area_L2", and the elevation value correction amount Δh = +0.3 m.

[0165] Figure 2 The following is a schematic diagram of the hardware entity of a topographic surveying and mapping system provided by an embodiment of the present invention. As Figure 2 shown, the hardware entity of the topographic surveying and mapping system 1000 includes: a processor 1001 and a memory 1002. Among them, the memory 1002 stores a computer program that can run on the processor 1001, and when the processor 1001 executes the program, it implements the steps in the method of any of the above embodiments.

Claims

1. A topographic surveying method based on GNSS-RTK, characterized in that, The method includes: Obtaining multi-source positioning signal data within a target area, where the multi-source positioning signal data includes satellite raw observations, receiver antenna phase center deviation parameters, and real-time kinematic differential correction information; Performing multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression; Inputting the phase observation sequence after interference suppression into a terrain feature calculation model to extract terrain elevation change feature parameters of the target area; Performing spatial fusion processing based on the terrain elevation change feature parameters and reference elevation parameters in a preset terrain database to generate three-dimensional terrain surface data of the target area; Dynamically adjusting the topographic map update frequency according to the difference degree between the three-dimensional terrain surface data and historical topographic survey results, and outputting real-time topographic survey results.

2. The method according to claim 1, wherein The performing multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression includes: Obtaining the carrier phase measurement value and pseudorange measurement value in the satellite raw observations, and respectively constructing a carrier phase fluctuation sequence and a pseudorange noise distribution sequence; Performing antenna phase deviation compensation on the carrier phase fluctuation sequence based on the receiver antenna phase center deviation parameters to generate a compensated carrier phase sequence; Establishing an association model between a multipath effect influence factor and the pseudorange noise distribution sequence according to the horizontal direction positioning error parameter in the real-time kinematic differential correction information; Performing dynamic threshold segmentation on the association model through a sliding time window to identify mutation noise points in the pseudorange noise distribution sequence; Performing combined adjustment calculation on the compensated carrier phase sequence and the pseudorange measurement value after removing mutation noise points to generate the phase observation sequence after interference suppression.

3. The method according to claim 2, characterized in that The method further includes the training process of the terrain feature calculation model, including: Collecting multiple groups of training data in historical topographic survey tasks, where each group of training data includes a phase observation sequence sample after multipath interference suppression processing and lidar elevation verification data corresponding to the phase observation sequence sample; Constructing a deep neural network model, where the deep neural network model includes a signal feature extraction layer, an elevation mapping layer, and a terrain surface generation layer; Inputting the phase observation sequence sample into the signal feature extraction layer to extract phase fluctuation features, satellite elevation distribution features, and receiver motion state features; Performing feature fusion on the phase fluctuation features and the satellite elevation distribution features to generate a spatially correlated feature vector, and inputting the spatially correlated feature vector into the elevation mapping layer for nonlinear transformation to obtain a nonlinear transformation result; Performing weighted superposition on the receiver motion state features and the nonlinear transformation result to generate preliminary terrain elevation parameters; Performing surface smoothing processing on the preliminary terrain elevation parameters through the terrain surface generation layer to output predicted terrain surface data; Calculating the elevation difference loss value between the predicted terrain surface data and the lidar elevation verification data, and optimizing the parameters of the deep neural network model using the backpropagation algorithm until the elevation difference loss value is less than a preset threshold.

4. The method according to claim 3, characterized in that, Performing spatial fusion processing based on the terrain elevation change characteristic parameters and the reference elevation parameters in the preset terrain database to generate three-dimensional terrain surface data for the target area, including: Extracting the reference elevation parameters corresponding to the target area from the preset terrain database, where the reference elevation parameters include contour line distribution data, slope parameters, and surface cover type parameters; Performing spatial superposition of the elevation increment data in the terrain elevation change characteristic parameters and the contour line distribution data to generate a superimposed elevation distribution map; Performing terrain continuity correction on the superimposed elevation distribution map according to the slope parameters to eliminate elevation mutation areas and obtain a corrected elevation distribution map; Performing elevation compensation for the vegetation-covered areas on the corrected elevation distribution map based on the surface cover type parameters to generate compensated elevation distribution data; Converting the compensated elevation distribution data into a gridded elevation point set and generating the three-dimensional terrain surface data through a cubic spline interpolation algorithm.

5. The method according to claim 4, characterized in that Dynamically adjusting the topographic map update frequency according to the difference degree between the three-dimensional terrain surface data and the historical topographic survey results and outputting real-time topographic survey results, including: Calculating the absolute value of the elevation deviation at the corresponding positions in the three-dimensional terrain surface data and the historical topographic survey results; Selecting a topographic map update mode according to the preset range in which the absolute value of the elevation deviation is located. Specifically, when the absolute value of the elevation deviation is less than the first threshold, a timed polling update mode is adopted; when the absolute value of the elevation deviation is between the first threshold and the second threshold, an event-triggered update mode is adopted; when the absolute value of the elevation deviation is greater than the second threshold, a real-time streaming update mode is started; In the event-triggered update mode, dynamically adjusting the data sampling interval according to the elevation deviation change rate. Specifically, when the elevation deviation change rate exceeds the preset rate threshold, shortening the data sampling interval and increasing the data output frequency of the GNSS-RTK receiver; In the real-time streaming update mode, splitting the three-dimensional terrain surface data into multiple data blocks, and performing independent verification and interpolation processing on each data block through parallel computing nodes to generate processed data blocks; Performing spatial splicing and timestamp alignment on the processed data blocks and outputting real-time topographic survey results with spatio-temporal consistency marks.

6. The method according to claim 5, characterized in that, The calculating the absolute value of the elevation deviation at the corresponding positions in the three-dimensional terrain surface data and the historical topographic survey results includes: Extracting a historical elevation data set corresponding to the target area from the historical topographic survey results, where the historical elevation data set contains elevation point coordinates and corresponding elevation values collected at multiple historical time points; Performing spatial matching of the elevation point coordinates in the three-dimensional terrain surface data and the coordinates in the historical elevation data set to determine the successfully matched coordinate pairs; For each successfully matched coordinate pair, calculating the absolute value of the difference between the current elevation value and the historical elevation value; Statistically calculate the weighted average of the absolute values of the differences of all the successfully matched coordinate pairs, where the weights for weighting are determined according to the terrain stability coefficients of the regions where the coordinate points are located. Specifically, a first weight is set in geologically active regions, and a second weight is set in flat regions, and the first weight is less than the second weight; Multiply the weighted average by the terrain change sensitivity coefficient to obtain the absolute value of the elevation deviation; Among them, the determination process of the terrain change sensitivity coefficient includes: Obtain the geological structure activity monitoring data of the target region, and the monitoring data includes surface displacement rate, seismic activity frequency, and underground water level change parameters; Compare the surface displacement rate with the reference displacement rate to generate a displacement deviation index; Select a sensitivity adjustment factor according to the grade interval where the seismic activity frequency is located, and the grade intervals include a low seismic frequency interval, a medium seismic frequency interval, and a high seismic frequency interval; Calculate a hydrological influence factor based on the underground water level change parameters, and the hydrological influence factor is positively correlated with the amplitude of water level change; Normalize the displacement deviation index, the sensitivity adjustment factor, and the hydrological influence factor to generate a comprehensive influence coefficient; Perform a non-linear transformation on the comprehensive influence coefficient through an exponential function to generate the terrain change sensitivity coefficient.

7. The method according to claim 6, characterized in that, The independent verification and interpolation processing of each data block by the parallel computing nodes includes: Assign a unique spatial identification code to each data block, and the spatial identification code includes longitude range, latitude range, and elevation range information; Load a preset terrain data verification rule library in each parallel computing node, and the rule library includes elevation mutation detection rules, contour continuity rules, and slope consistency rules; Use a streaming data processing engine to perform point-by-point verification on the elevation points in the data block. Among them, when it is detected that any elevation point violates the elevation mutation detection rule, mark the any elevation point as an abnormal point and trigger an interpolation compensation mechanism; For the positions marked as abnormal points, use the sliding average of the surrounding elevation points for interpolation compensation to generate compensated elevation points; Send the verified and compensated data blocks to the central node so that the central node can sort and splice the data blocks according to the spatial identification code; During the splicing process, detect the elevation difference at the boundary of adjacent data blocks. If the elevation difference exceeds a preset splicing threshold, start a boundary smoothing algorithm to re-interpolate the boundary area.

8. The method according to claim 7, wherein The start of the boundary smoothing algorithm to re-interpolate the boundary area includes: Extract the elevation points within a preset width range at the boundary of adjacent data blocks to generate a boundary elevation point set; Calculate the average elevation value and elevation standard deviation of all points in the boundary elevation point set; Based on the average elevation value, construct a Gaussian distribution model to generate the elevation confidence weight of each elevation point; Perform weighted fusion on the elevation points in the boundary area according to the elevation confidence weight to generate a fused boundary elevation value; Use the bilinear interpolation method to perform a smoothing transition process on the fused boundary elevation value to eliminate the splicing gap; Update the processed boundary elevation values to adjacent data blocks and recalculate the elevation differences in the stitching area.

9. The method according to claim 8, wherein The generation process of the spatio-temporal consistency tag includes: Attach a data acquisition timestamp and a spatial coordinate range identifier to each verified data block; Verify the topological relationship of the spatial coordinate range identifiers of all data blocks under the same timestamp until there are no overlapping or missing areas; Generate a processing pipeline number according to the processing order of the data blocks, and the number includes a node identifier and a data processing stage code; Combine and encode the timestamp, the spatial coordinate range identifier, and the processing pipeline number to generate a globally unique spatio-temporal consistency tag; During the process of outputting the real-time topographic mapping result, embed the spatio-temporal consistency tag into the metadata segment of the data file and perform hash binding with the three-dimensional terrain surface data.

10. A topographic mapping system, comprising a memory and a processor, the memory storing a computer program that can run on the processor, characterized in that, When the processor executes the program, it implements the steps in the method according to any one of claims 1 to 9.

Images