A terrain surveying and mapping method and system based on GNSS-RTK

Abstract

This invention provides a GNSS-RTK-based topographic mapping method and system. The method includes: acquiring multi-source positioning signal data within a target area, including raw satellite observations, receiver antenna phase center deviation parameters, and real-time dynamic 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 interference-suppressed phase observation sequence into a topographic feature calculation model to extract topographic elevation change feature parameters of the target area; performing spatial fusion processing based on the topographic elevation change feature parameters and benchmark elevation parameters in a preset topographic database to generate three-dimensional topographic surface data of the target area; dynamically adjusting the topographic map update frequency according to the difference between the three-dimensional topographic surface data and historical topographic mapping results, and outputting real-time topographic mapping results. This invention can balance accuracy, real-time performance, and resource efficiency in complex terrain environments.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a GNSS-RTK-based topographic mapping method and system. Background Technology

[0002] In existing technologies, GNSS-RTK-based topographic mapping systems typically employ a single data source processing mode to acquire topographic information. For example, they rely on carrier phase calculations from satellite observations to generate elevation data, or post-process the positioning results using static differential correction parameters. While these methods can achieve decimeter-level positioning accuracy in open, flat areas, in environments with significant multipath effects, such as densely built-up urban areas and vegetated areas, the accumulation of pseudorange noise caused by satellite signal reflection and diffraction leads to centimeter- to decimeter-level fluctuations in elevation calculation results. Some improved solutions attempt to suppress multipath interference by increasing the number of observed satellites or extending the observation time, but these methods struggle to meet the timeliness requirements of real-time topographic monitoring. Another approach loosely couples GNSS-RTK data with inertial navigation system data, which can temporarily improve positioning stability in dynamic environments, but the time-varying nature of sensor errors makes it impossible to guarantee the reliability of long-term continuous mapping. At the data fusion level, traditional methods linearly overlay GNSS-RTK elevation data with historical terrain databases using fixed weights, failing to consider the impact of terrain slope variations and land cover types on the data fusion rules. This results in elevation jumps at the boundaries between vegetated and bare land areas in the fused terrain surface. Furthermore, existing topographic map update mechanisms generally employ a global update strategy with preset time intervals, failing to differentiate the data change sensitivity between geologically stable and active areas, leading to both wasted computational resources and delays in updating key areas. These technical shortcomings make it difficult for traditional GNSS-RTK mapping methods to balance accuracy, real-time performance, and resource efficiency in complex terrain environments, limiting their widespread application in highly dynamic scenarios such as disaster early warning and engineering monitoring. Summary of the Invention

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

[0004] On one hand, embodiments of the present invention provide a GNSS-RTK-based topographic mapping method, the method comprising: acquiring multi-source positioning signal data within a target area, the multi-source positioning signal data including original satellite observations, receiver antenna phase center deviation parameters, and real-time dynamic 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 topographic feature calculation model to extract topographic elevation change feature parameters of the target area; performing spatial fusion processing based on the topographic elevation change feature parameters and benchmark elevation parameters in a preset topographic database to generate three-dimensional topographic surface data of the target area; dynamically adjusting the topographic map update frequency according to the difference between the three-dimensional topographic surface data and historical topographic mapping results, and outputting real-time topographic mapping results.

[0005] On the other hand, the present invention provides a topographic mapping system, including a memory and a processor, wherein the memory stores a computer program that can run on the processor, and the processor executes the program to implement the steps in the above-described method.

[0006] The GNSS-RTK-based terrain mapping method provided by this invention acquires multi-source positioning signal data within a target area. This multi-source positioning signal data includes raw satellite observations, receiver antenna phase center deviation parameters, and real-time dynamic differential correction information. It comprehensively utilizes the spatial distribution characteristics of satellite signals, receiver hardware error parameters, and the correction capabilities of real-time differential data. The method performs multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence with suppressed interference, effectively eliminating multipath signal aliasing caused by urban building reflections and vegetation obstruction. The phase observation sequence is then input into a terrain feature calculation model to extract terrain elevation changes. The system utilizes deep learning models to adaptively analyze the correlation between signal propagation patterns and elevation changes in complex terrain. Based on the terrain elevation change feature parameters and the baseline elevation parameters of a preset terrain database, it performs spatial fusion processing to generate three-dimensional terrain surface data. Through contour line distribution correction, slope continuity optimization, and vegetation cover compensation, it achieves spatial consistency fusion of multi-source terrain data. The system dynamically adjusts the update frequency according to the difference between the three-dimensional terrain surface data and historical results and outputs real-time terrain mapping results. It can intelligently switch between timed polling, event triggering, and real-time streaming update modes according to the actual rate of terrain change, optimizing system resource allocation while ensuring mapping accuracy. In this way, the raw satellite observations can provide high-frequency positioning reference information, the receiver phase center deviation parameter can compensate for inherent hardware errors, and the real-time dynamic differential correction information can eliminate environmental interference such as atmospheric delay. The three work together to form the 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 micro-change trend of surface elevation from phase fluctuation features through strong correlation learning with lidar verification data during the training phase. Spatial fusion processing solves the inherent defect of insufficient GNSS signal penetration through the vegetation cover area compensation mechanism, ensuring the continuity of elevation data in bare surfaces and vegetated areas. The difference-driven dynamic update strategy integrates the terrain stability index and change sensitivity coefficient into the threshold judgment 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, it achieves a comprehensive improvement in mapping accuracy, real-time performance, and resource efficiency in complex terrain environments, expanding the practical boundaries of GNSS-RTK technology in the field of dynamic terrain monitoring. Attached Figure Description

[0007] Figure 1 This is a schematic diagram illustrating the implementation process of a GNSS-RTK-based topographic mapping method provided in an embodiment of the present invention.

[0008] Figure 2 This is a schematic diagram of the hardware entity of a topographic mapping system provided in an embodiment of the present invention. Detailed Implementation

[0009] This invention provides a GNSS-RTK-based topographic mapping method, which can be executed by the processor of a topographic mapping system. The topographic mapping system can refer to a server, laptop, tablet, or desktop computer with data processing capabilities.

[0010] Figure 1 This is a schematic diagram illustrating the implementation process of a GNSS-RTK-based topographic mapping method provided in an embodiment of the present invention, as shown below. Figure 1 As shown, the method includes the following steps:

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

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

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

[0014] Step S200 is implemented to eliminate the interference of multipath effects caused by signal reflection on positioning data. Multipath interference refers to the superposition of multiple signals formed after GNSS signals are reflected by obstacles such as the ground surface and buildings during propagation. This causes 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, specifically including the following core operations: First, time-frequency analysis is performed on the carrier phase measurement sequence in the original satellite observations to identify the periodic fluctuation patterns caused by multipath effects; second, a multipath noise model is constructed based on the receiver motion state parameters, and the real signal and reflected signal components are separated through an adaptive filtering algorithm; finally, the spatial correlation correction of residual noise is performed by combining the horizontal positioning error parameters in the real-time dynamic differential correction information.

[0015] In practical implementation, generating the phase observation sequence after interference suppression requires processing multi-source positioning signal data in stages. First, antenna phase center deviation compensation is performed on the carrier phase measurements. The original phase observations are corrected in three-dimensional space using the receiver antenna phase center deviation parameter to eliminate systematic deviations caused by antenna physical characteristics. Second, for multipath noise in pseudorange measurements, a 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 within each segment, abrupt noise points exceeding the dynamic threshold are identified and marked as outliers for removal. Subsequently, the compensated carrier phase sequence and the filtered pseudorange measurements are jointly adjusted. The adjustment model incorporates receiver motion state parameters (such as velocity and acceleration) as constraints, and the optimal phase observation solution is solved using a least squares optimization algorithm. During this process, the horizontal direction error parameter in the real-time dynamic differential correction information is used to construct the correlation matrix between multipath effects and pseudorange noise. Matrix factorization techniques are used to extract noise components related to multipath and separate them from the original observation sequence. The final output phase observation sequence after interference suppression is a high-precision carrier phase dataset that has undergone bias compensation, noise removal and adjustment optimization. Its data format is a continuous phase observation sequence with timestamp alignment.

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

[0017] Step S300 converts the phase observation data into elevation feature parameters characterizing terrain undulations using a mathematical model. The terrain feature calculation model is a deep learning-based multimodal data fusion architecture. Its input is the phase observation sequence after interference suppression, and its output is the elevation change gradient, surface curvature, and slope distribution parameters within the target area. The model first analyzes the spatial correlation features in the phase observation sequence through a signal feature extraction layer, including the spatial distribution characteristics of carrier phase fluctuation amplitude, satellite elevation angle distribution pattern, and receiver motion trajectory. Then, it models the nonlinear relationship between signal features and terrain elevation through an elevation mapping layer, generating preliminary elevation estimates. Finally, it smooths and optimizes the preliminary estimates through a terrain surface generation layer, eliminating local noise and enhancing the expression of terrain continuity. In specific implementation, the extraction of terrain elevation change feature parameters requires phase observation sequence processing in stages. First, the signal feature extraction layer of the terrain feature calculation model uses a convolutional neural network (CNN) structure to extract temporal and spatial features from the input phase observation sequence. Temporal features, including short-term fluctuation trends and long-term variation cycles of phase observations, are calculated using a one-dimensional convolutional kernel sliding operation. Spatial features are extracted by analyzing the spatial geometric distribution of different satellite signals (e.g., the combination of satellite elevation and azimuth angles), specifically using a Graph Convolutional Network (GCN) to model the spatial topology of the satellite-receiver system. Next, the elevation mapping layer fuses the extracted signal features with pre-defined terrain prior knowledge (e.g., the elevation distribution patterns of typical landforms), generating preliminary elevation change parameters, including elevation increment, elevation change rate, and local curvature, through a fully connected neural network. Finally, the terrain surface generation layer employs a physically constrained surface optimization algorithm to jointly optimize the preliminary elevation parameters with terrain continuity priors (e.g., the consistency of elevation gradients in adjacent areas), outputting terrain elevation change feature parameters. These feature parameters are stored in a rasterized data format, with each raster cell containing an elevation value, elevation change direction, and confidence index, used for subsequent 3D terrain reconstruction.

[0018] Step S400: Spatial fusion processing is performed based on the terrain elevation change characteristic parameters and the benchmark elevation parameters in the preset terrain database to generate three-dimensional terrain surface data of the target area.

[0019] Step S400 is used to fuse the real-time calculated elevation change features with historical benchmark topographic data to construct a high-precision 3D topographic model. The preset topographic database is a standardized dataset containing historical surveying results of the target area, and its benchmark elevation parameters include contour line distribution data, slope parameters, and land cover type parameters. The contour line distribution data is a contour line vector map obtained through aerial photogrammetry or lidar scanning, the slope parameters are slope raster data calculated based on the Digital Elevation Model (DEM), and the land cover type parameters are land cover category labels (such as vegetation, water, and bare soil) obtained through remote sensing image classification. The core technologies of spatial fusion processing include elevation data overlay, topographic continuity correction, and land cover compensation.

[0020] In practice, generating 3D terrain surface data requires multi-level spatial data fusion. First, the elevation increment data from the terrain elevation change characteristic parameters is spatially overlaid with contour line distribution data from a pre-defined terrain database. The overlay process employs raster-to-vector conversion technology, converting the contour line vector data into a raster format with the same resolution as the elevation increment data, and generating the overlaid elevation distribution map through pixel-by-pixel addition. Second, terrain continuity correction is performed on the overlaid elevation distribution map based on slope parameters: by calculating the elevation gradient between adjacent raster cells, areas of abrupt elevation changes that do not conform to the natural terrain evolution (such as the illusion of steep cliffs caused by data noise) are identified, and an anisotropic diffusion algorithm is used to smooth these abrupt changes. Subsequently, elevation compensation for vegetation cover areas is performed on the corrected elevation distribution map based on land cover type parameters: for densely vegetated areas such as forests and shrubs, a pre-defined vegetation layer thickness compensation value (e.g., 0.5 to 3 meters) is added to the original elevation value. This compensation value is dynamically adjusted using an empirical database of vegetation type and height relationships. Finally, the compensated elevation distribution data is converted into a gridded elevation point set, and continuous three-dimensional terrain surface data is generated through a cubic spline interpolation algorithm.

[0021] Step S500: Dynamically adjust the topographic map update frequency based on the difference between the 3D topographic surface data and the historical topographic mapping results, and output the real-time topographic mapping results.

[0022] The difference degree refers to the quantitative index of the elevation deviation between the current 3D terrain surface data and historical terrain mapping results at the same spatial location. Its calculation process includes statistical analysis of the absolute value of elevation deviation, terrain stability weighting, and sensitivity coefficient correction. The strategy of dynamically adjusting the topographic map update frequency is based on a difference degree threshold grading mechanism: when the difference degree is below the first threshold, a timed polling update mode is adopted, updating data at fixed time intervals (such as 24 hours); when the difference degree is between the first and second thresholds, the system switches to an event-triggered update mode, initiating local updates only when significant terrain changes are detected; when the difference degree exceeds the second threshold, a real-time streaming update mode is enabled, using parallel computing nodes to process and verify the data for the entire region in real time.

[0023] In practice, outputting real-time topographic mapping results requires multi-level data processing and scheduling operations. First, when calculating the difference, the 3D topographic surface data and historical topographic mapping results need to be spatially matched: coordinate transformation unifies the elevation point sets of both to the same coordinate system, and nearest neighbor interpolation is used to align spatial resolution. Second, for each matched elevation point pair, the absolute deviation between the current elevation value and the historical value is calculated point by point, and a weighted average is performed based on the topographic stability coefficient of the area where the point is located (e.g., 0.3 for geologically active areas and 0.7 for flat areas) to obtain a preliminary difference index. Subsequently, the difference is corrected using a topographic change sensitivity coefficient: this sensitivity coefficient is calculated from geological tectonic activity monitoring data of the target area (e.g., surface displacement rate, earthquake frequency), and its function is to amplify the difference weight in tectonically active areas. Based on the corrected difference value, the system automatically selects the update mode and triggers the corresponding processing flow. In real-time streaming update mode, 3D terrain surface data is divided into multiple spatial data blocks. Each data block is assigned to an independent parallel computing node for outlier detection and interpolation compensation. The processed data blocks are then spatially stitched and timestamped to generate a final result with spatiotemporal consistency markers. The spatiotemporal consistency markers include data acquisition time, spatial range, and processing pipeline number, and are bound to the terrain data using a hash algorithm to ensure data integrity and traceability.

[0024] As one implementation method, step S200 involves performing multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression, which may specifically include:

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

[0026] Raw satellite observations are the unprocessed, raw signal measurement results output by the Global Navigation Satellite System (GNSS) receiver, specifically including two core data types: carrier phase measurements and pseudorange measurements. Carrier phase measurements refer to the phase difference along the satellite carrier signal propagation path recorded by the receiver. Their values ​​consist of integer and fractional periods, reflecting the precise length variation of the signal propagation path, possessing millimeter-level accuracy but suffering from integer ambiguity. Pseudorange measurements refer to the geometric distance between the satellite and the receiver calculated based on the signal propagation time. Their values ​​are obtained through modulation code phase measurements, exhibiting unambiguous characteristics but being susceptible to multipath effects and ionospheric delay. The process of constructing a carrier phase fluctuation sequence involves extracting the carrier phase measurements of each satellite from raw satellite observations over a continuous time period, arranging them chronologically to form a time series, and then eliminating the common error between the receiver clock bias and the satellite clock bias through differential operations to generate sequence data reflecting short-term carrier phase fluctuations. The process of constructing the pseudorange noise distribution sequence is as follows: pseudorange measurements within the same time period are classified by satellite number; the pseudorange measurements for each satellite are time-series processed; and low-frequency trend terms and high-frequency noise terms are separated by moving average filtering. The high-frequency noise term is then extracted to construct the pseudorange noise distribution sequence. During this process, it is crucial 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 a compensated carrier phase sequence.

[0028] Step S220 is used to eliminate the systematic carrier phase deviation introduced by the receiver antenna hardware characteristics. 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 documents provided by the manufacturer. The operation process of antenna phase deviation compensation is as follows: First, according to the receiver antenna model, the corresponding phase center deviation parameter is retrieved from the pre-stored antenna parameter database. This parameter is usually stored in the form of three-dimensional offset in the Northeast-Northeast-Upper-Side (ENU) coordinate system, represented as ΔX, ΔY, and ΔZ. Second, each measurement value in the carrier phase fluctuation sequence is combined with its corresponding satellite azimuth and elevation information. The phase center deviation parameter is projected onto the satellite signal propagation direction through a coordinate transformation model to calculate the direction-related phase deviation correction. Finally, the calculated phase deviation correction is subtracted from the original carrier phase fluctuation sequence one by one to generate the compensated carrier phase sequence. For example, the phase deviation correction ΔL1(ti) of satellite G01 at time ti can be calculated using 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, hardware-related errors in the carrier phase fluctuation sequence are effectively suppressed, and the compensated carrier phase sequence more accurately reflects the physical changes in the signal propagation path.

[0029] Step S230: Based on the horizontal positioning error parameters in the real-time dynamic differential correction information, establish a correlation model between the multipath effect influencing factor and the pseudorange noise distribution sequence.

[0030] The horizontal positioning error parameter in real-time dynamic differential correction information refers to the horizontal (eastward and northward) positioning residuals calculated by the base station network using real-time differential technology. This parameter reflects the remaining positioning deviation of the base station receiver after eliminating common errors, and its value is affected by both multipath effects and atmospheric delay residual errors. 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 path length of the reflected signal and the material of the reflector. The process of establishing the correlation model includes the following steps: First, extract the time series of the horizontal positioning error parameter from the real-time dynamic differential correction information; then, synthesize the eastward error sequence E(t) and the northward error sequence N(t) to generate the horizontal error magnitude sequence M(t) = √(E(t)). 2 +N(t) 2Secondly, the pseudorange noise distribution sequence P(t) and the horizontal error magnitude sequence M(t) are synchronized in time. Finally, a correlation model between P(t) and M(t) is constructed using multiple regression analysis, specifically in the form P(t) = k·M(t) + C + ε, where k is the proportionality coefficient, C is a constant term, and ε is random noise. Through this model, the multipath effect factor k can quantitatively characterize the contribution of horizontal positioning error to pseudorange noise.

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

[0032] Step S240 is used to detect anomalous abrupt changes in the pseudorange noise distribution sequence to support noise removal. A sliding time window refers to a data segmentation mechanism that slides along the time axis with a fixed step size. Its window length is dynamically set based on the time-varying characteristics of the signal sampling rate and multipath interference; for example, a 5-second window length can be used at a 1Hz sampling rate. The dynamic threshold segmentation process is as follows: First, the pseudorange noise distribution sequence is divided into multiple subsequences according to the sliding time window. Statistical characteristic analysis is performed on the noise data within each subsequence, calculating its mean μ and standard deviation σ. Second, a dynamic threshold is set for each subsequence based on the statistical characteristics; for example, a threshold of μ ± 3σ is set. Finally, noise points exceeding the dynamic threshold in the subsequences are detected and marked as abrupt noise points. During this process, the threshold setting strategy needs to be optimized based on the output results of the correlation model: for subsequences with a high multipath effect influence factor k, the threshold is appropriately lowered to enhance the sensitivity of abrupt change detection; conversely, for subsequences with a low k value, the threshold is increased to avoid false detections.

[0033] Step S250: Perform joint adjustment calculations on the compensated carrier phase sequence and the pseudorange measurements after removing abrupt noise points to generate a phase observation sequence after suppressing interference.

[0034] Step S250 aims to fuse the compensated high-precision carrier phase data with the filtered pseudorange observation data, and generate an anti-interference optimized phase observation sequence through a joint adjustment algorithm. Joint adjustment calculation refers to incorporating measurements from different observation types into a unified mathematical model for overall solution. Its advantage lies in leveraging the high precision of the carrier phase and the unambiguous nature of the pseudorange to complement each other and improve solution accuracy. The specific implementation process includes the following steps: First, establish a joint adjustment model for the carrier phase observation equation and the pseudorange observation equation, where the carrier phase equation is expressed as λ·Φ=ρ+c·(dT-dt)+Trop+Iono+ε Φ The pseudorange equation is expressed as P = ρ + c·(dT - dt) + Trop + Iono + ε PWhere λ is the carrier wavelength, Φ is the carrier phase measurement, ρ is the geometric distance, dT and dt are the receiver clock bias and satellite clock bias, respectively, Trop is the tropospheric delay, Iono is the ionospheric delay, and ε is the ionospheric delay. Φ With ε P The first step involves considering the observation noise of the carrier phase and pseudorange, respectively. The second step involves substituting the compensated carrier phase sequence and the pseudorange measurements after removing abrupt noise points into the joint adjustment model. Receiver motion state parameters (such as velocity and acceleration) are introduced as constraints, and the least squares estimation algorithm is used to solve for the optimal integer ambiguity combination and position correction. Finally, the calculated fixed integer ambiguity values ​​are substituted back into the carrier phase observation equation to generate the phase observation sequence after interference suppression.

[0035] As one implementation method, the method also includes a training process for the terrain feature calculation model, which may specifically include the following steps:

[0036] Step S301: Collect multiple sets of training data from historical topographic mapping tasks. Each set of training data includes phase observation sequence samples processed by multipath interference suppression and lidar elevation verification data corresponding to the phase observation sequence samples.

[0037] Step S301 is used to construct a high-quality training dataset, the core of which is to obtain multi-source terrain observation data pairs with spatiotemporal consistency. Historical terrain mapping tasks refer to previously completed and validated mapping projects. The data includes complete raw observation data streams from the Global Navigation Satellite System, along with high-precision lidar elevation verification data. Simultaneously, multipath interference suppression processing has been completed, generating standardized phase observation sequences. The phase observation sequence samples processed by multipath interference suppression refer to carrier phase time-series data processed by the method in step S200. The data format is a phase observation value matrix arranged by timestamps, with matrix dimensions including satellite number, observation epoch, and phase compensation value. Lidar elevation verification data refers to a Digital Elevation Model (DEM) generated from three-dimensional point cloud data acquired through airborne or ground-based lidar systems. Its spatial resolution must reach sub-meter level (e.g., 0.5 meters), elevation accuracy better than ±5 centimeters, and it must be strictly aligned with the phase observation sequence samples in both time and space. In practice, data acquisition requires the following steps: First, select task records from the historical surveying project database that meet the time span requirements, ensuring that the selected tasks cover diverse terrain types (such as mountains, plains, and urban built-up areas). Second, perform multipath interference suppression processing step S200 on the raw observation data of the Global Navigation Satellite System in each task to generate phase observation sequence samples. Simultaneously, 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. Finally, perform spatiotemporal alignment processing on the phase observation sequence samples and 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, which 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-tiered structure. The signal feature extraction layer is responsible for parsing multi-level spatiotemporal features from the phase observation sequence. Its input is a phase observation matrix from multiple satellites and multiple epochs, and its output is a feature vector containing phase fluctuation patterns, satellite geometric distribution characteristics, and receiver motion state. The elevation mapping layer is used to establish a nonlinear mapping relationship between signal features and terrain elevation. It realizes the conversion from feature space to elevation parameters through a fully connected network and attention mechanism. The terrain surface generation layer focuses on reconstructing discrete elevation estimates into continuous terrain surfaces, and uses a physical constraint-based interpolation algorithm and a smoothing filter to eliminate local errors. In practical implementation, the model construction needs to be implemented in modules. The signal feature extraction layer can adopt a hybrid architecture of convolutional neural network and long short-term memory network. One-dimensional convolutional kernel (size 3×1) is used to extract local fluctuation features of phase observation sequence, and LSTM unit is used to capture the dynamic pattern of satellite elevation angle distribution changing over time. Secondly, the elevation mapping layer can be composed of four fully connected network layers, with the number of neurons in each layer decreasing sequentially (e.g., 512→256→128→64). The activation function adopts the rectified linear unit (ReLU), and residual connection is introduced in the output layer to accelerate convergence. Finally, the terrain surface generation layer integrates deconvolutional network and thin plate spline (TPS) interpolation algorithm. The deconvolutional network is responsible for upsampling the feature vector to the target spatial resolution, and the TPS algorithm generates a smooth surface based on 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 angle distribution features and receiver motion state features.

[0041] Phase fluctuation characteristics refer to the short-term fluctuation patterns of carrier phase observations over time. Their values ​​reflect subtle changes in surface reflection characteristics and signal propagation paths, and are obtained through sliding calculations using a one-dimensional convolutional kernel across the time dimension. Satellite elevation angle distribution characteristics refer to the statistical characteristics of the elevation angles of the participating satellites in the sky over time. Their values ​​affect the intensity of multipath effects and the signal-to-noise ratio (SNR). Spatial correlation features are extracted by constructing a satellite-epoch elevation angle matrix and applying a Graph Convolutional Network (GCN). Receiver motion state characteristics include receiver velocity, acceleration, and trajectory curvature parameters. These values ​​are derived from the dynamic positioning results of the Global Navigation Satellite System (GNSS), and the temporal evolution of the motion state is captured using an LSTM network. In practical implementation, feature extraction requires channel-specific processing: First, the phase observation sequence samples are reshaped into a three-dimensional tensor [number of satellites × number of epochs × 1 channel], input to the CNN module of the signal feature extraction layer, and multi-scale phase fluctuation features are extracted through three convolutional layers (with filters of 16, 32, and 64 respectively), outputting a feature map with dimensions of [number of satellites × number of epochs × 64 channels]; Second, the satellite elevation angle distribution features are generated by constructing an elevation angle-time matrix (with dimensions 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 angle distribution feature vector; Simultaneously, the receiver motion state features (including eastward velocity, northward velocity, elevation velocity, and three-axis acceleration) are arranged in epoch time order as a time series, input to the LSTM network (with 64 hidden layer units), and outputting a 128-dimensional motion state feature vector.

[0042] Step S304: Perform feature fusion between phase fluctuation features and satellite elevation angle distribution features to generate spatially correlated feature vectors, and input the spatially correlated feature vectors into the elevation mapping layer for nonlinear transformation to obtain the nonlinear transformation result.

[0043] Step S304 integrates the temporal characteristics of phase fluctuations with the spatial distribution characteristics of satellite elevation angles to construct spatially correlated features reflecting changes in terrain elevation. Feature fusion employs a strategy combining cascading and attention mechanisms: First, the phase fluctuation feature matrix (size [number of satellites × number of epochs × 64 channels]) is max-pooled along the satellite dimension, compressed into a temporal feature sequence of [number of epochs × 64 channels]. Second, the satellite elevation angle distribution feature vector is expanded to a sequence of [number of epochs × 32 channels], concatenated with the phase fluctuation temporal feature sequence along the channel dimension to generate a preliminary fused feature of [number of epochs × 96 channels]. Finally, a self-attention mechanism is introduced to calculate the correlation weights between features of different epochs, and a spatially correlated feature vector of [number of epochs × 96 channels] is generated through weighted summation. The nonlinear transformation process of the elevation mapping layer is as follows: spatially relevant feature vectors are input into a four-layer fully connected network in epochal order. The first layer maps 96-channel features 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 nonlinear transformation result. During this process, batch normalization and ReLU activation functions are applied after each fully connected network layer, and residual connections are added between the third and fourth layers to suppress gradient vanishing.

[0044] Step S305: Weight the receiver motion state characteristics and the nonlinear 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 weight superposition adopts a combination of feature concatenation and adaptive weighting: First, the receiver motion state feature vector (size

[128] ) is copied and extended along the epoch dimension to the same size as the nonlinear transformation result ([100×128]); Second, the extended motion state features and the nonlinear transformation result (size [100×64]) are concatenated along the feature dimension to generate a fusion feature matrix of [100×192]; Subsequently, the fusion features are channel-weighted by a learnable weight matrix. The weight matrix is ​​dynamically generated by a fully connected network. Its input is the motion state features and the nonlinear transformation result of the current epoch, and its output is a 192-dimensional weight vector; Finally, the weighted fusion features are input to the fully connected layer (output dimension is 1) to generate the elevation change parameters corresponding to each epoch. After aggregation according to spatial location, a preliminary terrain elevation parameter matrix is ​​formed.

[0046] Step S306: The preliminary terrain elevation parameters are smoothed by the terrain surface generation layer, and the predicted terrain surface data is output.

[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 employs a combined algorithm of Thin Plate Spline (TPS) interpolation and Anisotropic Diffusion Filter: First, each grid cell in the initial terrain elevation parameter matrix is ​​treated as a control point, and the TPS algorithm is applied to generate a continuous surface, whose energy function minimizes bending energy to achieve global smoothness. Second, an anisotropic diffusion coefficient matrix is ​​constructed based on land cover type data (such as vegetation, bare soil, and water bodies), and edge-preserving smoothing is performed on the TPS surface, i.e., the smoothing intensity is enhanced in homogeneous areas (such as flat farmland), while detailed features are preserved in heterogeneous areas (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, upsampling the initial 0.5-meter resolution elevation matrix to 0.2 meters, and then outputting high-resolution terrain surface data 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 model parameters through supervised learning to minimize the solution error. The elevation difference loss 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 location, and its calculation formula is: Where N is the number of effective grid points, is the predicted value, These are the verification values ​​for the LiDAR. The backpropagation algorithm uses the Adaptive Moment Estimator (Adam) optimizer with an initial learning rate of 1e-4 and a learning rate decay strategy (decay factor 0.5 every 10 rounds). During training, the following steps are required: First, input the batch data (e.g., 32 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 gradient of the parameters of each layer through backpropagation; finally, update the model parameters and repeat the iteration until the loss value converges to a preset threshold (e.g., 1e-3 meters). 2 ).

[0050] As one implementation method, step S400 involves spatially fusing terrain elevation change characteristic parameters with benchmark elevation parameters in a preset terrain database to generate three-dimensional terrain surface data of the target area. Specifically, this may include:

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

[0052] Step S410 is used to call basic geographic information data from a preset terrain database to support terrain fusion calculation. The preset terrain database is a standardized geospatial database containing historical surveying results of the target area. Its baseline elevation parameters are generated from multi-source geographic data through standardization processing, specifically including contour line distribution data, slope parameters, and land cover type parameters. The contour line distribution data is a contour line vector dataset obtained through aerial photogrammetry or lidar scanning technology. Each contour line represents a continuous spatial distribution of the same elevation value. Its data attributes include elevation values, coordinate point sequences, and accuracy level information. The storage format is a Shapefile or GeoJSON file conforming to the Open Geospatial Consortium (OGC) standard. The slope parameters are a raster dataset generated based on a Digital Elevation Model (DEM) using a slope calculation algorithm. Each raster cell stores a slope value (in degrees or percentages), reflecting the degree of surface inclination. The calculation uses the third-order inverse-squared difference method to balance computational efficiency and accuracy. Land cover type parameters are raster data of land cover categories obtained through multispectral remote sensing image classification. The classification system follows international geographic information standards (such as ISO 19144) and includes categories such as vegetation, water bodies, bare soil, and buildings. Each raster cell stores the category code and classification confidence level. In practice, the extraction operation requires the following steps: First, a spatial query is performed from a pre-defined terrain database based on the boundary coordinates of the target area (such as a minimum bounding rectangle or a geofence polygon) to filter out data blocks of benchmark elevation parameters that fully contain or partially overlap the target area. Second, the filtered data blocks undergo format conversion and coordinate system unification processing. The contour line distribution data is converted from vector format to the same raster format (such as GeoTIFF) as the slope parameters and land cover type parameters, and all data are converted to the same plane coordinate system (such as WGS84 UTM) and elevation benchmark (such as the EGM2008 geoid). Finally, the converted data is resolution aligned, and the spatial resolution of the slope parameters and land cover type parameters is adjusted to be consistent with the contour line distribution data (such as 0.5 meters) using bilinear interpolation.

[0053] Step S420: Spatially overlay the elevation increment data in the topographic elevation change characteristic parameters with the contour line distribution data to generate an overlaid elevation distribution map.

[0054] Step S420 integrates the real-time calculated elevation change data with historical benchmark topographic data to generate an updated elevation distribution map. The elevation increment data in the topographic elevation change feature parameters refers to the rasterized elevation change output by the topographic feature calculation model in step S300. Its value represents the elevation change of each location in the target area relative to the historical benchmark within the current calculation cycle (positive values ​​indicate uplift, negative values ​​indicate subsidence). The data format is a floating-point raster with the same resolution as the benchmark elevation parameters. Spatial overlay operation uses a raster algebra method. The specific process is as follows: First, the contour line distribution data is converted from a vector format to the same raster format as the elevation increment data. During the conversion, each raster cell is assigned the corresponding contour line elevation value, and areas not directly covered by contour lines are filled using linear interpolation. Second, the elevation increment data and the rasterized contour line data are added pixel-by-pixel 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) The baseline elevation value is ΔH(x,y), and the elevation increment value is ΔH(x,y). Finally, the data validity of the overlay results is verified, and outliers caused by coordinate offset or resolution mismatch (such as abrupt values ​​that exceed the reasonable elevation range) are removed.

[0055] Step S430: Based on the slope parameters, perform terrain continuity correction on the superimposed elevation distribution map to eliminate areas of abrupt elevation changes, and obtain the corrected elevation distribution map.

[0056] Step S430 is used to eliminate abrupt elevation changes caused by data noise or fusion errors, ensuring that the terrain surface conforms to the continuity of natural landforms. The terrain continuity correction adopts a slope-constrained anisotropic diffusion filter algorithm. Its principle is to dynamically adjust the filtering intensity according to the slope parameter, retain detailed features in steep slope areas, and enhance the smoothing effect in gentle slope areas. The specific process includes: First, calculating 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, calculating the anisotropic diffusion coefficient matrix based on the slope parameter, where 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 areas with a slope greater than k (such as cliffs) to suppress smoothing, while the diffusion coefficient increases in areas with a slope less than k to promote smoothing; Finally, applying an explicit iterative method to perform diffusion filtering on the elevation distribution map, with the number of iterations dynamically set according to the severity of elevation abrupt changes (e.g., 5 to 20 times), and 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 vegetation cover areas on the corrected elevation distribution map based on the land cover type parameter to generate compensated elevation distribution data.

[0058] Step S440 is used to correct the elevation underestimation error caused by signal penetration effect in vegetation-covered areas. Elevation compensation for vegetation-covered areas adopts a category-related compensation strategy, with the compensation value dynamically determined based on an empirical database of vegetation type and elevation. The specific implementation steps are as follows: First, vegetation category subclasses (such as coniferous forest, broadleaf forest, and shrubs) are extracted from the land cover type parameters, and a preset vegetation layer thickness compensation value is configured for each subclass (e.g., 0.8 meters for coniferous forest, 1.2 meters for broadleaf forest, and 0.3 meters for shrubs). Second, the corrected elevation distribution map is spatially aligned with the vegetation cover type raster, and elevation compensation calculation is performed on each raster cell in the vegetation-covered area. Finally, edge transition processing is performed on the compensated elevation data, using Gaussian filtering to smoothly transition the elevation values ​​at the boundary between vegetation and non-vegetated areas, avoiding abrupt changes in compensation values.

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

[0060] Step S450 converts raster-formatted elevation data into a continuous surface model to support 3D visualization and analysis. A gridded elevation point set refers to a set of discrete elevation points arranged at regular intervals. Its data format is a sequence of triplets containing longitude, latitude, and elevation values, with the grid spacing consistent with the resolution of the compensated elevation distribution data (e.g., 0.5 meters). The conversion process includes the following operations: First, the compensated elevation distribution raster data is traversed in row and column order, generating a corresponding elevation point for each valid raster cell. The coordinates are calculated based on the raster origin coordinates, pixel size, and row and column indices. Second, a topological check is performed on the elevation point set to remove isolated points caused by missing or invalid data, and missing areas are filled using nearest neighbor interpolation. 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) must be set during the interpolation process to ensure smooth transitions at the surface edges. The final output 3D terrain surface data is stored in Triangulated Irregular Network (TIN) format. Each triangle vertex contains precise coordinates and elevation values, along with interpolation accuracy metadata, which can be directly imported into a Geographic Information System (GIS) platform for 3D rendering or terrain analysis.

[0061] As one implementation method, step S500 involves dynamically adjusting the topographic map update frequency based on the difference between the three-dimensional topographic surface data and historical topographic mapping results, and outputting real-time topographic mapping results. Specifically, this may include:

[0062] Step S510: Calculate the absolute value of the elevation deviation between the three-dimensional terrain surface data and the corresponding location in the historical terrain mapping results.

[0063] Step S510 is used to quantify the spatial difference between the current terrain data and the historical benchmark. Its 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 of the same spatial location in the 3D terrain surface data and the elevation value in the historical terrain mapping results. Its calculation must meet the spatiotemporal consistency constraints: 1) the spatial coordinate system and the elevation benchmark must be strictly unified; 2) the historical mapping results corresponding to the timestamp must be comparable to the collection period of the current terrain data. The specific implementation process includes: First, extracting the historical elevation dataset of the target area from the historical terrain mapping results. This dataset contains the coordinates of elevation points collected at multiple historical time points and their corresponding elevation values, stored in a spatiotemporal four-dimensional array (longitude, latitude, elevation, time); Second, spatially matching the elevation point coordinates in the 3D terrain surface data with the coordinates in the historical elevation dataset. The matching method uses bilinear interpolation to resample the historical data to the resolution grid of the current terrain data; Subsequently, calculating the elevation difference for each successfully matched coordinate pair, using the formula Δh=|h current(x,y) -h historical(x,y) |, where h current(x,y) h is the current elevation value. historical(x,y) The historical elevation values ​​are used. Finally, a weighted average is calculated based on the absolute values ​​of the differences using a topographic stability coefficient. This coefficient is set according to the geological activity level; for example, a lower weight of 0.3 is assigned to landslide-prone areas (geologically active areas) to weaken the impact of instantaneous changes, while a higher weight of 0.7 is assigned to plain areas (topographically stable areas) to enhance the significance of long-term trends. The weighted average formula is Δh. avg =Σ(w i ·Δh i ) / Σw i , where w i Let Δh be the terrain stability coefficient of the area where coordinate point i is located. i This represents the elevation difference. Furthermore, by introducing a terrain change sensitivity coefficient (calculated from surface displacement rate, seismic activity frequency, and groundwater level change parameters) to nonlinearly amplify the weighted average, the final expression for the absolute value of the elevation deviation is Δh. abs =Δh avg ×k sense , where k sense This is the sensitivity coefficient.

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

[0065] Step S520 aims to dynamically switch topographic map update strategies to balance data freshness and computational resource consumption. The preset range is defined by a first threshold α1 and a second threshold α2, whose values ​​are set according to the sensitivity of the target area to topographic changes. For example, α1=5 cm and α2=15 cm are set in ordinary areas, while α1=2 cm and α2=8 cm are set in disaster-sensitive areas. The timed polling update mode updates the entire topographic data at fixed time intervals (e.g., 24 hours), suitable for scenarios with slow or stable topographic changes; the event-triggered update mode initiates incremental updates when a local elevation deviation exceeds α1 but is below α2, processing only data in the changed area; the real-time streaming update mode performs millisecond-level continuous processing and publishing of the entire area's data, suitable for scenarios such as geological disaster emergency response. In specific implementation, the update mode selection logic is as follows: First, the calculated absolute value of the elevation deviation Δh... abs Compared with α1 and α2; if Δh abs If α1 < Δh, the timed polling update mode is activated, and the system calls a preset timed task scheduler (such as a Cron job) to trigger the data update process at a period of T (e.g., T = 86400 seconds); if α1 ≤ Δh abs <α2, switch to event-triggered update mode, the system registers an elevation change event listener, when the Δh of a specific grid changes... abs If α1 is exceeded three times consecutively, an event message is generated and local data resampling and fusion are triggered; if Δh abs ≥α2, immediately start the real-time streaming update mode, enable high-frequency data output of the GNSS-RTK receiver (e.g., increase from 1Hz to 10Hz), and allocate parallel computing resources to perform streaming processing.

[0066] Step S530: In event-triggered update mode, dynamically adjust the data sampling interval according to the elevation deviation change rate. 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 adaptively optimize the data acquisition frequency in event-triggered update mode. The elevation deviation change rate is defined as the time derivative of the current absolute value of the elevation deviation with respect to the previous measurement, and is calculated using the formula v=Δh. abs (t2)-Δh abs (t1) / (t2-t1), where the 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, β = 0.5 cm / hour in a general monitoring scenario and β = 2 cm / hour in a landslide early warning scenario. The dynamic adjustment strategy includes: maintaining the basic sampling interval T when v ≤ β. base(e.g., 2 hours); when v>β, shorten the sampling interval according to the rate exceeding the limit, the formula is T new =T base The data output frequency of the GNSS-RTK receiver is increased from the conventional 1Hz to ceil(v / β)×1Hz (ceil is the floor function). During implementation, the data transmission bandwidth and storage buffer size need to be adjusted synchronously to ensure stable processing of the high-frequency data stream. Furthermore, the system continuously monitors the rate change trend; if v falls below β and remains below it for more than 3 sampling periods, it is gradually restored to the basic sampling parameters.

[0068] Step S540: In real-time streaming update mode, the 3D terrain surface data is divided into multiple data blocks. Each data block is independently verified and interpolated by parallel computing nodes to generate processed data blocks.

[0069] Step S540 aims to achieve real-time distributed processing of large-scale terrain data. Data segmentation employs a spatial partitioning strategy, specifically: First, based on the spatial extent of the target area (e.g., longitude span Δλ, latitude span Δφ) and the number of computation nodes N, the 3D terrain surface data is divided into N sub-blocks, each containing a longitude range [λ]. start +(i-1)×Δλ / N, λ start +i×Δλ / N] and latitude range [φ 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 for longitude. startThe starting value for latitude is determined; secondly, a unique spatial identifier is assigned to each data block, with the encoding rule being "lower longitude limit_upper longitude limit_lower latitude limit_upper latitude limit". For example, the identifier "116.300E_116.305E_39.900N_39.905N" represents a data block with longitudes from 116.300° to 116.305° and latitudes from 39.900° to 39.905°. Parallel computing nodes load a pre-defined terrain data verification rule base, which includes rules for detecting elevation abrupt changes (such as a single point elevation change exceeding three times the standard deviation of the average of adjacent points), rules for contour line continuity (such as the elevation difference between adjacent contour lines must not exceed a preset interval), and rules for slope consistency (such as the deviation between local slope and regional average slope must not exceed 20%). The verification and processing flow is as follows: Each node uses a streaming data processing engine (such as Apache Flink) to scan the elevation points in the data block point by point. When a point with a sudden change in elevation is detected, it is marked as an anomaly and an interpolation compensation mechanism is triggered. The interpolation method uses the moving average of the elevation points in the surrounding 8 neighborhoods to calculate the compensation value. After processing, the node sends the verification log and the corrected data block to the central node.

[0070] Step S550: Spatial stitching and timestamp alignment are performed on the processed data blocks to output real-time terrain mapping results with spatiotemporal consistency markers.

[0071] Step S550 is used to integrate the distributed processing results and ensure the consistency of the spatiotemporal attributes of the data. The spatial stitching operation includes: First, sorting the processed data blocks according to their spatial identifiers, arranging them in ascending order of longitude and descending order of latitude; Second, using a seamless stitching algorithm to eliminate differences in the boundaries between blocks, specifically by extracting a 5-pixel-wide overlapping area at the boundaries of adjacent data blocks and calculating the average elevation difference Δh in the overlapping area. overlap If Δh overlap If the splicing threshold γ is less than or equal to the splicing threshold γ (e.g., γ = 0.1 meters), splice directly; if Δh overlap If the value is greater than γ, a Gaussian weighted average fusion is performed on the overlapping area, with the weights inversely proportional to the distance from the boundary. Timestamp alignment must ensure that the processing time deviation of all data blocks does not exceed the system clock precision (e.g., 1 millisecond). This is implemented by recording the processing start time t in the data block metadata. start With end time t end and with t endThis serves as the unified timestamp for the data block. The generation process of the spatiotemporal consistency marker is as follows: the spatial identifier, timestamp, and processing pipeline number (such as "Node03_Phase2") of each data block are combined into a globally unique identifier, for example, "116.300E-116.305E_39.900N-39.905N_20231012T153045Z_Node03_Phase2", and a fixed-length digest value (such as a 64-bit hexadecimal string) is generated using the SHA-256 hash algorithm and embedded into the metadata segment of the data file. The final output real-time terrain mapping result is in the format of GeoPackage or LASer (LAS) point cloud file, containing complete 3D coordinates, elevation values, timestamps, and hash markers, supporting direct loading and integrity verification by GIS platforms.

[0072] As one implementation method, step S510: Calculating the absolute value of the elevation deviation between the three-dimensional terrain surface data and the corresponding location in the historical terrain mapping results, specifically may include:

[0073] Step S511: Extract the historical elevation dataset corresponding to the target area from the historical topographic mapping results. The historical elevation dataset contains the coordinates of elevation points collected at multiple historical time points and their corresponding elevation values.

[0074] Step S511 involves spatiotemporal retrieval and standardization of historical topographic data. Historical topographic mapping results are verified topographic datasets obtained in the past through lidar, photogrammetry, or GNSS-RTK technology. They are stored in a spatiotemporal four-dimensional database (longitude, latitude, elevation, and timestamp), with each data entry containing the elevation value of a unique coordinate point and its acquisition time. The extraction process involves the following steps: First, a spatial range query is performed in the historical database based on the boundary range of the target area (e.g., the coordinates of the vertices of the minimum bounding rectangle or the geofence polygon), filtering out all historical elevation points whose spatial locations are within the target area. Second, the filtered data is sorted by timestamp, retaining datasets from at least three different historical time points to support trend analysis, such as selecting three survey results spanning January 2020, June 2021, and March 2023. Subsequently, coordinate system unification is performed on the multi-period historical elevation data, transforming it to the same plane coordinate system (e.g., WGS84 UTM Zone 50N) and elevation datum (e.g., EGM2008 geoid) as the current 3D terrain surface data, eliminating systematic biases caused by datum differences. Finally, a spatiotemporal index structure (e.g., R-tree or quadtree) is constructed to store the historical elevation dataset hierarchically by spatial location and timestamp, supporting efficient spatial matching and time series analysis.

[0075] Step S512: Spatial match the coordinates of elevation points in the 3D terrain surface data with the coordinates in the historical elevation dataset to determine the successfully matched coordinate pairs.

[0076] Step S512 aims to achieve precise spatial alignment between current terrain data and historical data. Spatial matching refers to determining whether elevation points in the 3D terrain surface data and points in the historical elevation dataset are in the same spatial location, using a coordinate tolerance threshold, within the same coordinate system. Specifically, the operation involves: first, aligning the coordinates (X, Y, X) of the elevation points in the 3D terrain surface data... current ,Y current ) and coordinates (X) in historical elevation dataset historical Y historical The point-to-point distance is calculated using the formula D=√[(X current - X historical ) 2 +(Y current - Y historical ) 2 Secondly, set a spatial tolerance threshold δ (usually half the data resolution, such as 0.25 meters). If D ≤ δ, then the match is considered successful, and a coordinate pair ((X...) is generated. current Y current ),(X historical Y historical If D > δ, then mark it as an unmatched point and exclude it from subsequent calculations. For datasets with inconsistent resolutions, the low-resolution data needs to be resampled using bilinear interpolation to match the high-resolution data. For example, when the current terrain data resolution is 0.2 meters and a certain historical data resolution is 1 meter, the historical data is interpolated to a 0.2-meter grid to ensure that all coordinate points are comparable. Successfully matched coordinate pairs are stored in a list structure, with each entry containing the current elevation value h. current Historical elevation value h historical And 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 range of a single point. The formula for calculating the absolute value of the difference is Δh_ i =|h current _ i - h historical _ i |, where h current _ i h is the elevation value of the i-th matching point in the 3D terrain surface data. historical _ iThis refers to the elevation values ​​corresponding to historical elevation points. During the calculation process, if the historical elevation dataset contains data from multiple periods, the historical elevation value with the closest time interval to the current data should be selected for calculation. For example, if the current data collection time is October 2023, then historical data from March 2023 should be used instead of data from 2020. Outliers caused by sensor malfunctions or environmental interference (such as points where the elevation value exceeds the geographically reasonable range or abrupt changes exceed 10 meters) should be excluded.

[0079] Step S514: Calculate the weighted average of the absolute values ​​of the differences of all successfully matched coordinate pairs. The weights are determined based on the terrain stability coefficient of the area where the coordinate points are located. In geologically active areas, a first weight is set, and in flat areas, a second weight is set. The first weight is less than the second weight.

[0080] Step S514 eliminates the interference of inherent topographic characteristics on the difference assessment by using a weighted average. The topographic stability coefficient is a pre-set weight value based on the geological activity classification of the target area, defined as follows: geologically active areas (such as fault zones and landslide-prone areas) are assigned a weight w1 (e.g., 0.3), flat areas (such as plains and stable sedimentary areas) are assigned a weight w2 (e.g., 0.7), and transitional areas (such as hills and gentle slopes) are assigned a moderate weight w3 (e.g., 0.5). The formula for calculating the weighted average is Δh. avg =(Σ(w i ×Δh i )) / Σw i , where w i Let Δh be the terrain stability weight for the i-th coordinate point. i This corresponds to the absolute value of the difference. The implementation process includes: First, querying the classification label of the region where each matching point is located from the preset terrain stability zoning map and mapping it to the corresponding weight value; Second, traversing all matching points, accumulating the absolute value of the weighted difference and the sum of the weights; Finally, calculating the weighted average and storing it as a scalar value.

[0081] Step S515: Multiply the weighted average value 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 using a sensitivity coefficient. The terrain change sensitivity coefficient k... sense This is a dimensionless parameter, typically ranging from 1.0 to 3.0, and is calculated using steps S5151 to S5156. The absolute value of the elevation deviation Δh abs The calculation formula 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 assessment results based on regional geological risk, ensuring that small changes in sensitive areas (such as around volcanoes) are effectively amplified, while larger changes in stable areas are appropriately suppressed.

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

[0084] The implementation of step S5151 involves the acquisition and preprocessing of multi-source geological monitoring data. Geological tectonic activity monitoring data includes: 1) Surface displacement rate, acquired through GNSS Continuously Operating Reference Stations (CORS) or Synthetic Aperture Radar Interferometry (InSAR) technology, measured in millimeters per year, reflecting the intensity of crustal deformation; 2) Seismic activity frequency, extracted from the historical seismic event catalog of the target area from the seismic network database, and the number of earthquakes with a magnitude ≥ 2.0 within a unit of time (e.g., one year); 3) Groundwater level change parameters, calculated from the pressure sensor data of groundwater monitoring wells, measured in meters per month. Data preprocessing includes: removing seasonal trends from 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 (e.g., the number of earthquakes per month); and 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 degree of anomaly in surface deformation. Reference displacement rate v base This is a reference value set based on the regional geological background, typically taken from the long-term (e.g., 10-year) average deformation rate. Displacement deviation index I disp The calculation formula is I disp =(v current -v base ) / v base , where v current This represents the current monitored displacement rate. For example, if a certain region v base =5 mm / year, v current =12.5 mm / year, then I disp =(12.5-5) / 5=1.5. An index greater than 0 indicates accelerated deformation, while an index less than 0 indicates decelerated deformation.

[0087] Step S5153: Select the sensitivity adjustment factor according to the level range of the seismic activity frequency. The level range includes low frequency range, medium frequency range and high frequency range.

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

[0089] Step S5154: Calculate the hydrological impact factor based on the groundwater level change parameters. The hydrological impact factor is positively correlated with the magnitude of water level change.

[0090] Step S5154 is used to quantify the impact of groundwater level changes on surface stability. Hydrological Influence Factor I hydro The calculation formula is I hydro =1+0.2×|Δh water |, where Δh water This represents the average monthly variation in water level (unit: meters). For example, if Δh water =0.2 meters, 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] Step S5155 involves eliminating dimensional differences and integrating the effects 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 of α; I hydro _ norm =(I hydro -1.0) / (2.0-1.0), is the normalized I. hydro Value. Overall 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 nonlinear transformation on the comprehensive influence coefficient using an exponential function to generate the terrain change sensitivity coefficient.

[0094] Step S5156 maps the comprehensive influence coefficient to the sensitivity coefficient range. 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 changes in areas of high influence. For example, C comb When k = 0.4705, 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 step S515 to calculate the absolute value of the elevation deviation, thereby dynamically amplifying the sensitivity to terrain changes.

[0095] As one implementation method, in step S540, each data block is independently verified and interpolated using parallel computing nodes, which may specifically include:

[0096] Step S540: In real-time streaming update mode, the 3D terrain surface data is divided into multiple data blocks. Each data block is independently verified and interpolated by 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. Data blocks are subsets of 3D terrain surface data divided by spatial extent, with the division rules based on the longitude, latitude, and elevation distribution characteristics of the target area. Specifically, firstly, the maximum allowable size of a single data block is determined based on 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 area of ​​0.1° × 0.1° (approximately 11 km × 11 km) and an elevation range of -100 meters to 9000 meters. Secondly, a spatial grid partitioning algorithm is used to cut the original 3D terrain surface data into multiple non-overlapping regular cubic data blocks. During the cutting process, it is ensured that the boundaries of the data blocks are aligned with the geographic coordinate grid lines to avoid cross-blocks. Data redundancy is addressed. Finally, a unique spatial identifier is generated for each data block. This identifier consists of a six-tuple: lower longitude limit, upper longitude limit, lower latitude limit, upper latitude limit, lower elevation limit, and upper elevation limit. For example, the identifier "E116.300-E116.400_N39.900-N40.000_H-100-H9000" represents a data block with longitude from 116.300° to 116.400°, latitude from 39.900° to 40.000°, and elevation from -100 meters to 9000 meters. After receiving a data block, the parallel computing nodes load a pre-defined terrain data verification rule base. This rule base includes rules for detecting elevation abrupt changes (e.g., the difference between a single point's elevation and the average of its eight neighboring points exceeds three times the standard deviation), rules for contour line continuity (e.g., the distance between adjacent contour lines must not exceed 1.5 times the preset threshold), and rules for slope consistency (e.g., the absolute value of the deviation between the local slope and the regional average slope is less than 20%). In the verification and interpolation process, the streaming data processing engine (e.g., Apache Kafka Streams) scans the data blocks point by point according to the timestamps of the elevation points. When an elevation point is detected to violate the elevation abrupt change detection rule, it is immediately marked as an anomaly and the interpolation compensation mechanism is triggered. For example, if a point's elevation is 152.3 meters and the average of its eight neighboring points is 148.1 meters with a standard deviation of 1.2 meters, then the three-times-standard-deviation threshold is 148.1 ± 3.6 meters. This point is marked as an anomaly because it exceeds the upper limit. Interpolation compensation uses a neighborhood moving average algorithm. For example, taking a 5×5 window of neighborhood elevation points centered on the outlier, removing the highest and lowest 10% of outliers, and calculating the arithmetic mean of the remaining points as the compensation value. After processing, all outliers within the data block are corrected. The verification log records the outlier type, location, and correction value. The corrected data block is then sent to the central node with an appended spatial identifier and processing timestamp.

[0098] Step S541: Assign a unique spatial identifier to each data block. The spatial identifier contains 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 identifier is a string encoding conforming to the ISO 19115 geographic information metadata standard. Its generation rule is as follows: the lower and upper limits of longitude, latitude, and elevation of the data block are concatenated in a fixed format, with each field separated by an underscore. Numerical precision is retained to three decimal places. Longitude and latitude are represented in decimal units, and elevation is in meters. For example, a data block covering longitude 116.300° to 116.400°, latitude 39.900° to 40.000°, and elevation -100 meters to 500 meters has a spatial identifier code of "E116.300-E116.400_N39.900-N40.000_H-100-H500". In practice, the allocation operation needs to be performed synchronously with the data segmentation process: when the 3D terrain surface data is cut into cubic data blocks, the boundary coordinates of each block are calculated in real time and a corresponding spatial identifier code is generated. This identifier code is embedded in the metadata header of the data block and registered with the spatial index service of the central node. The functions of the spatial identifier code include: 1) quickly locating the geographical extent of the data block in parallel computing nodes; 2) guiding the sorting and spatial relationship reconstruction of data blocks during the data stitching stage; and 3) supporting the rapid association of the original data location when tracing the source of abnormal data.

[0100] Step S542: Load the preset terrain data verification rule base into each parallel computing node. The rule base includes elevation change detection rules, contour line continuity rules, and slope consistency rules.

[0101] Step S542 is used to deploy a unified data quality control standard in the distributed computing nodes. The terrain data verification rule base is a configuration file in XML or JSON format, and its content is defined as follows: 1) The elevation change 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 μ8 of its 8 neighboring points exceeds 3 times the neighborhood standard deviation σ8, that is, |h(x,y)-μ8|>3σ8, then P is determined to be an outlier point; 2) The contour line continuity rule requires that the elevation difference Δh between adjacent contour lines be greater than 3 times the neighborhood standard deviation σ8. contour It must not exceed 1.5 times the contour interval D, i.e., Δh contour ≤1.5D, if a section of contour line is detected to violate this rule, then that section is marked as a fracture area; 3) The slope consistency rule is defined as follows: the local slope value S(x,y) and the regional average slope S avg The absolute value of the deviation must not exceed 20%, i.e., |S(x,y)-S avg | / S avg≤0.2. The rule base 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 an independent data detection thread for each rule. For example, the elevation change detection thread scans the input data stream in real time, calls the neighborhood statistics 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 the elevation points in the data block point by point. When any elevation point is detected to violate the elevation change detection rule, the elevation point is marked as an anomaly and the interpolation compensation mechanism is triggered.

[0103] The specific implementation of step S543 can leverage the high-throughput real-time computing capabilities of streaming data processing engines (such as Apache Flink or Spark Streaming). The streaming verification process is as follows: Data blocks are split into a stream of elevation points arranged in row-major order. Each point contains longitude, latitude, elevation value, and timestamp. The engine reads each point sequentially, retrieves its neighboring points (e.g., a 5×5 window) based on spatial coordinates, and calculates neighborhood statistics (mean, standard deviation). Subsequently, the elevation value and statistics of the current point are substituted into the elevation change detection rules for judgment. If the rules are violated, an anomaly flag is inserted into the metadata of that point (e.g., setting the anomaly flag to 1), and an interpolation compensation event is triggered. After the event is triggered, the engine pauses the processing of the current data stream, calls the interpolation compensation module to generate compensation values, and writes the corrected elevation values ​​back into the data stream.

[0104] Step S544: For the locations marked as outliers, interpolation compensation is performed using the moving average of surrounding elevation points to generate compensated elevation points.

[0105] The interpolation compensation algorithm in step S544 can employ robust statistical methods to reduce the impact of outliers. Specifically, the operation could be as follows: Select a circular neighborhood with a radius R (e.g., R = 5 pixels) centered on the outlier, and extract all elevation points within this area; sort the neighborhood points by elevation value, and remove the top 10% and bottom 10% extreme values; calculate the arithmetic mean of the remaining points as the compensation value, using the formula h. comp =Σh i / (N-2k), where N is the total number of neighborhood points, k=floor(0.1×N), h comp h is the compensated elevation value. iThis represents the elevation values ​​of each point within the neighborhood. For example, a 5×5 neighborhood of an outlier contains 25 points. After removing the two highest values ​​(e.g., 155.6 meters, 154.9 meters) and the two lowest values ​​(e.g., 145.2 meters, 146.0 meters), the average of the remaining 21 points is 148.3 meters. This average is used as the compensation value to replace the original outlier value of 152.3 meters. The compensated elevation points retain their original coordinates and timestamps, and the compensation operation type (e.g., "moving average compensation_R5") and the number of neighborhood points involved in the calculation are recorded 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 identifier code.

[0107] The data transmission protocol in step S545 must ensure the integrity and timing consistency of the processed data blocks. Specifically, each parallel computing node, after processing a data block, encapsulates it into a data packet according to the lexicographical order of its spatial identifier (longitude ascending → latitude ascending → elevation ascending) and sends it to the central node via TCP. Upon receiving the data packet, the central node parses the spatial identifier and registers it in the global spatial index table, which 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 based on the longitude and latitude range of the spatial identifier. For example, the identifier "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, awaiting splicing instructions.

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

[0109] Step S546 is used to eliminate inter-block seam problems caused by distributed processing. The preset splicing threshold γ is set according to data accuracy requirements; for example, γ = 0.1 meters means that the elevation difference between corresponding points at the boundaries of adjacent blocks must not exceed 0.1 meters. The detection process is as follows: extract a 2-pixel-wide overlapping region at the boundaries 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| at all corresponding points within the overlapping region. A(x,y) -h B(x,y) The percentage of points with a Δh exceeding γ is counted. If this percentage exceeds 5%, it is considered a severe seam problem, triggering the boundary smoothing algorithm. For example, if the overlapping area between adjacent blocks A and B contains 100 points, and 8 of these points have a Δh > 0.1 meters (8% > 5%), smoothing processing is required.

[0110] As one implementation, step S546, which involves initiating a boundary smoothing algorithm to re-interpolate the boundary region, may include:

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

[0112] The preset width in step S5461 can be set to, for example, 2 to 5 pixels wide of the data block boundary. For example, for 0.5-meter resolution data, the preset width is 2 pixels, or 1.0 meter. The extraction operation includes: extracting all elevation points from the column containing the maximum longitude and the column preceding it (a total of 2 columns) from the eastern boundary of data block A; simultaneously, extracting all points from the column containing the minimum longitude and the column following it (a total of 2 columns) from the western boundary of data block B, and merging them to generate a set of boundary elevation points. For example, if the longitude of the eastern boundary of data block A is 116.400°, then all points from longitude 116.395° to 116.400° (2-pixel width) are extracted; if the longitude of the western boundary of data block B is 116.400°, then points from 116.400° to 116.405° are extracted. 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 standard deviation of elevation of all points in the boundary elevation point set.

[0114] The statistical calculations in step S5462 aim 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; elevation standard deviation σ boundary =√[Σ(h i -μ boundary ) 2 / (N-1)],h i This refers to the elevation value of each point in the boundary elevation point set. For example, the boundary set contains 200 points with an average elevation μ = 150.2 meters and a standard deviation σ = 0.3 meters. This statistical value is used for subsequent Gaussian model construction and confidence weight calculation.

[0115] Step S5463: Based on the average elevation value, construct a Gaussian distribution model and generate the elevation confidence weight for 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π)), confidence weight w i =P(h i ) / P max , where P maxLet P(μ) be the peak value of the distribution (P(μ) = 1 / (σ√(2π))), h be the elevation value, and σ be the standard deviation. The weights are normalized to the range of 0 to 1, with points closer to μ having higher weights.

[0117] Step S5464: Weighted fusion of elevation points in the boundary area based on elevation confidence weights to generate fused boundary elevation values.

[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 with w Bi h represents the weights of the corresponding points from data blocks A and B, respectively. Ai h is the elevation value of data block A. Bi This is the elevation value of data block B.

[0119] Step S5465: Use bilinear interpolation to smooth the transition of the merged boundary elevation values ​​and eliminate splicing gaps.

[0120] Step S5465, bilinear interpolation, constructs a continuous elevation transition surface within the boundary region. Specifically, the merged boundary elevation points are used as control points, and an interpolation grid is generated along the east-west (longitude) and north-south (latitude) directions. For any point (x, y) to be interpolated, its elevation h(x, y) = a×x + b×y + c×x×y + d. The coefficients a, b, c, and d are solved by fitting the control points using the least squares method. For example, an interpolation grid with a resolution of 0.2 meters is generated within a 1.0-meter-wide boundary region, allowing the elevation value of data block A to smoothly transition from 150.4 meters to 150.1 meters in data block B.

[0121] Step S5466: Update the processed boundary elevation values ​​to the adjacent data blocks and recalculate the elevation differences of the spliced ​​area.

[0122] For example, the interpolated boundary elevation values ​​are written back to the corresponding positions in data blocks A and B, respectively, overwriting the original boundary data; then, the elevation points of the boundary region are extracted again, and Δh is calculated. new =|h_A new (x,y)-h_B new (x,y)|,h_A new (x,y) and h_B new (x, y) represent the updated elevation values ​​of data block A and data block B at coordinates (x, y), respectively, to verify whether the γ threshold is met. For example, the original Δh max=0.3 meters, after smoothing Δh max =0.05m < γ = 0.1m, the seam is eliminated. The updated data block is resent to the central node to participate in the construction of the global terrain model.

[0123] As one implementation method, the process of generating spatiotemporal consistency markers may specifically include:

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

[0125] Step S551 is used to assign precise spatiotemporal attribute identifiers to the terrain data blocks generated by distributed processing. The data acquisition timestamp is in Coordinated Universal Time (UTC) encoding conforming to the ISO 8601 standard, with no limit on precision, such as milliseconds. The format can be set to "YYYYMMDDThhmmss.sssZ", for example, "20231015T083045.123Z" represents October 15, 2023, 08:30:45.123. The spatial coordinate range identifier is generated based on the geographical coverage of the data block, including a six-tuple of lower longitude, upper longitude, lower latitude, upper latitude, lower elevation, and upper elevation. The numerical precision is retained to three decimal places. Longitude and latitude are represented in decimal units, 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 longitude from 116.300° to 116.400°, latitude from 39.900° to 40.000°, and elevation from -100 meters to 500 meters. In practice, the timestamp is synchronously generated by the GPS timing module of the data acquisition equipment, while the spatial coordinate range identifier is automatically calculated during the data block segmentation stage using a spatial grid partitioning algorithm. Additional operations are implemented by modifying the metadata segment of the data block: the timestamp is written to the "AcquisitionTime" field, and the spatial coordinate range identifier is written to the "Spatial Extent" field, both of which are encoded in binary to improve storage efficiency.

[0126] Step S552: 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.

[0127] Step S552 aims to ensure the spatial integrity and seamlessness of data blocks acquired at the same time. For example, topology verification includes two checks: 1) no spatial coverage omissions, meaning the union of the spatial coordinate range identifiers of all data blocks must completely cover the target area; 2) no spatial overlap, meaning the intersection of the spatial coordinate range identifiers of any two data blocks is empty. In implementation, firstly, the spatial coordinate range identifiers of all data blocks under the same timestamp are extracted from the spatiotemporal index library of the central node and arranged according to the rules of ascending longitude and then ascending latitude; then, a computational geometry library (such as GEOS) is used to perform topology 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 range is continuous, it is determined whether there are gaps or overlaps. If a gap is detected (e.g., the upper longitude limit of data block A is 116.400°, while the lower longitude limit of data block B is 116.405°), a region omission alarm is generated and the data re-sampling process is triggered; if an overlap is detected (e.g., the upper longitude limit of data block A is 116.400°, while the lower longitude limit of data block B is 116.395°), it is marked as a conflict region and the data block is re-segmented.

[0128] Step S553: ​​Generate a processing pipeline number based on the processing order of the data blocks. 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 rules for processing pipeline numbers are as follows: the node identifier uses a globally unique node ID assigned by the distributed computing framework (e.g., "Node07"). The data processing stage code consists of the English abbreviation of the processing stage and its sequence number (e.g., "VAL1" represents the first verification, "INT2" represents the second interpolation). The numbering logic is as follows: when a data block enters a 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 sequentially passes through the verification stage (VAL1) of node Node03 and the interpolation stage (INT2) of node Node12, then the pipeline numbering sequence is "Node03_VAL1_20231015T083045.123Z→Node12_INT2_20231015T083102.456Z". The number is stored in the "PipelineID" field of the data block metadata and records the multi-stage processing history in JSON array format. This numbering mechanism allows the problematic stage to be located by tracing back the pipeline number when data anomalies occur. For example, if a data block is marked as an anomaly in the VAL1 stage of node Node05, the logs of that 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 spatiotemporal consistency marker.

[0131] The implementation process of step S554 ensures the global uniqueness and information integrity of the tag through structured coding. The combination coding rule is as follows: the timestamp, spatial coordinate range identifier and processing pipeline number are concatenated in a fixed order, and the fields are separated by "#", 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 using the Base64 algorithm. For example, the original string is converted to "MjAyMzEwMTVUMDgzMDQ1LjEyM1ojRTExNi4zMDAtRTExNi40MDBfTjM5LjkwMC1OMzkuOTAwX0gtMTAwLUg1MDAjTm9kZTAzX1ZBTDEtPk5vZGUxMl9JTlQy". Global uniqueness is guaranteed by the following mechanisms: 1) Timestamps are accurate to the millisecond level to avoid conflicts between different batches of data; 2) Spatial coordinate range identifiers cover unique geographical areas; 3) Processing pipeline numbers include node IDs and timestamps to ensure differences in processing paths. For example, data blocks collected from the same geographical area at different times will inevitably 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 output of real-time terrain mapping results, the spatiotemporal consistency marker is embedded into the metadata segment of the data file and hash-bound with the three-dimensional terrain surface data.

[0133] Step S555 is used to ensure data integrity and tamper-proofing through digital signature technology. The embedding operation includes: writing the spatiotemporal consistency tag into the metadata segment of a standard geographic data format (such as GeoTIFF or LAS), specifically in the "CustomTags" field of the file header. The hash binding process is as follows: First, the SHA-256 algorithm is used to hash all elevation points of the 3D terrain surface data to generate a 64-bit hexadecimal digest value (such as "a1b2c3d4e5f6…"); then, the digest value is encrypted and merged with the spatiotemporal consistency tag to generate a digital signature string; finally, the signature is written into the "DigitalSignature" field of the data file. For example, the metadata of a 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 data block is maliciously modified during transmission (e.g., the elevation value of a point is changed from 150.3 meters to 160.3 meters), its hash value will change significantly, causing signature verification to fail and triggering a security alarm.

[0134] In an optional derivative implementation, after outputting the real-time terrain mapping results in step S500, the method provided in this embodiment of the invention may further include:

[0135] Step S600: Detect terrain anomalies in the real-time terrain mapping results and generate boundary coordinates and anomaly type identifiers for the anomalies.

[0136] Step S600 is used to identify anomalous areas in real-time terrain data that do not conform to the natural evolution patterns of landforms. The real-time terrain mapping result is three-dimensional terrain surface data with spatiotemporal consistency markers output by step S500, in the format of a gridded elevation point set or a triangular mesh model. Anomaly detection is achieved through multi-dimensional terrain feature analysis, specifically including elevation gradient change rate calculation, historical terrain stability assessment, and anomaly probability modeling. In practice, the elevation gradient change rate is first extracted from the real-time terrain data. This parameter is defined as the change in elevation value per unit horizontal distance, calculated using the formula G=√[(∂h / ∂x)]. 2 +(∂h / ∂y) 2 ], where ∂h / ∂x and ∂h / ∂y are the east-west and north-south elevation gradient components, respectively. Preset threshold Gthreshold The elevation gradient is dynamically set based on regional topographic features; for example, it might be set to 0.5 meters per meter in plains and 2.0 meters per meter in mountainous areas. When the elevation gradient change rate of a continuous region exceeds G... threshold When an anomaly is detected, it is marked as a candidate anomaly region. For example, if a region has a detected G=2.5 m / m within a 5×5 pixel area, exceeding the mountainous threshold, it is listed as a candidate anomaly. Subsequently, the historical terrain stability index S of the candidate region is obtained. index Its calculation formula is S index =1 / (σ 2 ×e (λΔt) ), where σ 2 Here, Δt represents the variance 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 compares the elevation gradient change rate G with S. index Combined, generate the regional anomaly probability value P. anomaly =α×G / G max +(1-α)×(1-S index ), where α is the weighting factor (e.g., α=0.6), G max This represents the maximum gradient value in the region. When P... anomaly Exceeding the dynamic threshold P threshold When the gradient is 0.7, it is determined to be a valid anomaly area, and anomaly type identifiers are generated by comparing the gradient change pattern with historical data: natural erosion identifier (gradient change and matching historical erosion pattern), artificial construction identifier (gradient change and matching building outline), and geological collapse identifier (gradient is distributed in concentric circles and accompanied by negative elevation change).

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

[0138] Step S700 is used to enhance the analysis of surface attributes in the anomalous area using multispectral remote sensing data. The multispectral remote sensing image data is a high-resolution image (such as Sentinel-2 or Landsat-8 data) containing visible, near-infrared, and shortwave infrared bands, and its spatial resolution must match the boundary coordinates of the anomalous area (e.g., 10 meters). Surface reflectance feature extraction includes: 1) Calculating the Normalized Difference Vegetation Index (NDVI) = (NIR - Red) / (NIR + Red), where NIR is the near-infrared reflectance and Red is the red light reflectance; 2) Calculating the Modified Soil Adjusted Vegetation Index (MSAVI) = (2 × NIR + 1 - √(2 × NIR + 1)) 2-8×(NIR-Red) )) / 2, used to reduce soil background interference; 3) Calculate the Normalized Difference Water Index (NDWI) = (Green-NIR) / (Green+NIR) to identify water body distribution. Vegetation cover density parameters are calculated through a pixel decomposition model, decomposing mixed pixels into the proportions of endmembers such as vegetation, bare soil, and water bodies. Vegetation cover density V cover =Vegetation end-member ratio × 100%.

[0139] Step S800: Match the surface reflectance features with the preset geological hazard feature database. When the match is successful, activate the terrain verification mode, send a high-density sampling command to the GNSS-RTK receiver, and trigger UAV-borne lidar collaborative mapping.

[0140] Step S800 aims to identify geological hazard risks and initiate a review mechanism through multispectral feature matching. The geological hazard feature database is a database containing spectral features of typical hazards (landslides, debris flows, collapses), stored in the form of a multidimensional vector set, where each vector contains NDVI, MSAVI, NDWI, and V. cover Historical statistical values ​​of parameters such as... The matching algorithm uses cosine similarity calculation: Sim=Σ(F... i ×F' i ) / √(ΣF i 2 ×ΣF' i 2 ), where F i Given the current surface reflectance characteristics, F' i This is a reference vector in the feature library. A successful match is determined when Sim ≥ 0.85, activating the terrain verification mode. For example, if the current feature vector has a similarity Sim = 0.89 with the "landslide" feature in the library, the verification process is triggered. The high-density sampling command increases the data acquisition frequency of the GNSS-RTK receiver from 1Hz to 10Hz, and simultaneously, the UAV-borne lidar starts scanning, setting the scanning parameters as follows: pulse frequency 500kHz, scan angle ±30°, and point density ≥ 50 points / square meter.

[0141] Step S900: Based on the spatial registration results of high-density sampling data and lidar point cloud data, the position coordinates of the boundary of the abnormal area are corrected, the corrected terrain mapping results are generated and updated to the three-dimensional terrain surface data.

[0142] Step S900 is used to improve the positioning accuracy of anomaly boundaries through multi-source data fusion. Spatial registration employs the Iterative Closest Point (ICP) algorithm, rigidly transforming (translation and rotation) 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 surface is generated, with its coordinate system unified to WGS84 UTM Zone 50N and its elevation reference to EGM2008. The dual-constraint adjustment calculation incorporates the GNSS-RTK positioning accuracy model (horizontal ±1cm, elevation ±2cm) and the LiDAR scanning error model (angular error ±0.01°), constructing 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 This refers to the coordinate values ​​measured by GNSS (Global Navigation Satellite System), x adj The adjusted coordinate values, w gnss w lidar These are weights calculated based on their respective error models. For example, after adjustment, the horizontal error of a boundary point decreased from ±5cm to ±1.5cm. Morphological closing operations use 3×3 circular structuring elements to perform dilation-erosion operations on the boundary point set, filling pores with a diameter of less than 2 meters. Anisotropic diffusion filtering is applied to the slope continuity transition processing, smoothing the gradient along the boundary normal direction to ensure a natural connection between the corrected boundary line and the surrounding terrain. Finally, the updated 3D terrain surface data incrementally replaces the elevation values ​​of the original anomaly areas and adds terrain classification labels (such as "Geological Subsidence Area_L2").

[0143] As one implementation method, step S600, which involves detecting terrain anomalies in the real-time terrain mapping results and generating boundary coordinates and anomaly type identifiers for the anomalies, may include:

[0144] Step S610: Extract continuous areas from the real-time terrain mapping results where the elevation gradient change rate exceeds a preset threshold, and mark them as candidate anomaly areas.

[0145] For example, performing Sobel operator convolution on real-time terrain data to calculate elevation gradients, and setting a dynamic threshold G. threshold =μ G +3σ G, where μ G For the region-averaged gradient, σ G This represents the standard deviation. For example, μ for a certain region. G =0.8m / m, σ 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 region of at least 5 × 5 pixels with an out-of-range gradient. Candidate abnormal regions are marked using a region growing algorithm. The coordinates of the candidate regions are stored as polygon vectors, and the attribute table records the mean gradient, area, and centroid position.

[0146] Step S620: Obtain the historical topographic stability index of the candidate anomaly area. The stability index is calculated by the variance of elevation changes and the time decay coefficient in the historical mapping data.

[0147] In step S620, the variance σ of historical elevation changes 2 This is obtained through multi-period DEM difference calculations. For example, by selecting three DEM data periods from the most recent five years, the elevation change Δh for each period can be calculated. t variance σ 2 =Σ(Δh t -μ Δh ) 2 / (n-1), μ Δh This refers to the average historical elevation change. The time decay coefficient λ = 0.01 / day. For example, σ in a certain area... 2 =0.25 meters 2 If Δt = 365 days, then S index =1 / (0.25×e (0.01×365) The result is 1 / (0.25×38.5)=0.104, indicating low stability.

[0148] Step S630: Compare the elevation gradient change rate with the historical terrain stability index in a weighted manner to generate the regional anomaly 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, S index =0.104, then P anomaly =0.6×(2.5 / 5.0)+0.4×(1-0.104)=0.3+0.358=0.658, no alarm is triggered when it is below the threshold of 0.7.

[0150] Step S640: When the probability value of an anomaly in a region exceeds the dynamic threshold, it is determined to be a valid anomaly region and an anomaly type identifier is generated.

[0151] Dynamic threshold P threshold =0.7, when P anomaly When the value is ≥0.7, an identifier is generated by combining the gradient direction distribution with the spectral feature matching results. For example, for a certain region P... anomaly =0.82, the gradient is radially distributed, and NDVI decreases by 0.2, which is identified as a geological collapse indicator.

[0152] Step S650: Determine the corresponding verification strategy based on the anomaly type identifier.

[0153] For example, verification strategies may include: using historical satellite imagery (e.g., at 5-year intervals) to detect changes in areas marked by natural erosion; querying local government construction permit databases to match coordinates with permit limits in areas marked by artificial construction; and activating ground-penetrating radar (GPR) or resistivity imaging equipment to scan underground cavities in areas marked by geological subsidence.

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

[0155] Step S910: Perform spatial coordinate system unification processing on the GNSS-RTK elevation points and lidar point cloud data in the high-density sampling data to generate a fused positioning reference surface.

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

[0157] Step S920: Perform dual-constraint adjustment calculations on the boundary coordinates of the abnormal area on the fusion positioning reference plane to eliminate GNSS-RTK signal jitter error and 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.02m, where T is the coordinate transformation matrix. After adjustment, the horizontal accuracy of the boundary points is improved to ±0.015m. gnss These are the coordinate values ​​observed by the GNSS-RTK receiver (planar coordinates x, y or three-dimensional coordinates x, y, z); x is the "true" coordinate value to be determined (the optimal estimate after adjustment); x lidar These are the original point cloud coordinates obtained by the lidar scan; 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 σ lidar It is the standard deviation of GNSS-RTK and lidar observations, representing the accuracy of the data.

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

[0160] For example, the morphological closing operation uses a 3×3 circular structuring element, first expanding (filling the 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 surrounding normal area terrain data to generate a smooth transition 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 value is 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, where d is the distance from the boundary, and D = 10 meters. S(x,y) is the slope value after mixing; S b S represents the slope value of the boundary region. n The slope value for the normal area; w b and w n represents the mixing ratio of the boundary slope (Sb) and the normal area slope (Sn), respectively; 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: Based on the corrected boundary line, re-divide the terrain into zones and update the elevation values ​​and terrain classification labels of the corresponding areas in the 3D 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 anomaly 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 is corrected by Δh = +0.3 meters.

[0165] Figure 2 This is a schematic diagram of the hardware entity of a topographic mapping system provided in an embodiment of the present invention, such as... Figure 2 As shown, the hardware entity of the terrain mapping system 1000 includes a processor 1001 and a memory 1002, wherein the memory 1002 stores a computer program that can run on the processor 1001, and the processor 1001 executes the program to implement the steps in the method of any of the above embodiments.

Claims

1. A topographic mapping method based on GNSS-RTK, characterized in that, The method includes: Acquire multi-source positioning signal data within the target area, including raw satellite observations, receiver antenna phase center deviation parameters, and real-time dynamic differential correction information; The multi-source positioning signal data is subjected to multipath interference suppression processing to generate a phase observation sequence after interference suppression; The phase observation sequence after interference suppression is input into the terrain feature calculation model to extract the terrain elevation change feature parameters of the target area; the terrain feature calculation model is a multimodal data fusion architecture based on deep learning; The baseline elevation parameters corresponding to the target area are extracted from the preset terrain database. The baseline elevation parameters include contour line distribution data, slope parameters, and land cover type parameters. The elevation increment data in the terrain elevation change characteristic parameters are spatially superimposed with the contour line distribution data to generate a superimposed elevation distribution map. Based on the slope parameters, the terrain continuity of the superimposed elevation distribution map is corrected to eliminate areas of abrupt elevation changes, resulting in a corrected elevation distribution map. Based on the land cover type parameter, the elevation of the vegetation cover area is compensated on the corrected elevation distribution map to generate compensated elevation distribution data. The compensated elevation distribution data is converted into a gridded elevation point set, and three-dimensional terrain surface data is generated by cubic spline interpolation algorithm. The topographic map update frequency is dynamically adjusted based on the difference between the three-dimensional topographic surface data and historical topographic mapping results, and real-time topographic mapping results are output.

2. The method according to claim 1, characterized in that, The step of performing multipath interference suppression processing on the multi-source positioning signal data to generate a phase observation sequence after interference suppression includes: Obtain the carrier phase measurement value and pseudorange measurement value from the original satellite observation values, and construct the carrier phase fluctuation sequence and pseudorange noise distribution sequence respectively; Based on the receiver antenna phase center deviation parameter, antenna phase deviation compensation is performed on the carrier phase fluctuation sequence to generate a compensated carrier phase sequence. Based on the horizontal positioning error parameters in the real-time dynamic differential correction information, a correlation model is established between the multipath effect influence factor and the pseudorange noise distribution sequence. Dynamic threshold segmentation of the correlation model is performed using a sliding time window to identify abrupt noise points in the pseudorange noise distribution sequence. The compensated carrier phase sequence and the pseudorange measurement value after removing abrupt noise points are jointly adjusted to generate the phase observation sequence after interference suppression.

3. The method according to claim 2, characterized in that, The method also includes a training process for the terrain feature calculation model, including: Multiple sets of training data were collected from historical topographic mapping tasks. Each set of training data included phase observation sequence samples processed by multipath interference suppression and lidar elevation verification data corresponding to the phase observation sequence samples. A deep neural network model is constructed, which includes a signal feature extraction layer, an elevation mapping layer, and a terrain surface generation layer. The phase observation sequence samples are input into the signal feature extraction layer to extract phase fluctuation features, satellite elevation angle distribution features, and receiver motion state features. The phase fluctuation characteristics and the satellite elevation angle distribution characteristics are fused to generate a spatial correlation feature vector, and the spatial correlation feature vector is input into the elevation mapping layer for nonlinear transformation to obtain the nonlinear transformation result. The receiver motion state characteristics and the nonlinear transformation results are weighted and superimposed to generate preliminary terrain elevation parameters. The preliminary terrain elevation parameters are smoothed using the terrain surface generation layer to output predicted terrain surface data. The elevation difference loss value between the predicted terrain surface data and the lidar elevation verification data is calculated, and the parameters of the deep neural network model are optimized using the backpropagation algorithm until the elevation difference loss value is less than a preset threshold.

4. The method according to claim 1, characterized in that, The step of dynamically adjusting the topographic map update frequency based on the difference between the three-dimensional topographic surface data and historical topographic mapping results, and outputting real-time topographic mapping results, includes: Calculate the absolute value of the elevation deviation between the three-dimensional terrain surface data and the corresponding location in the historical terrain mapping results; The topographic map update mode is selected based on the preset range of the absolute value of the elevation deviation. Specifically, when the absolute value of the elevation deviation is less than a first threshold, a timed polling update mode is used; when the absolute value of the elevation deviation is between the first threshold and a second threshold, an event-triggered update mode is used; and when the absolute value of the elevation deviation is greater than the second threshold, a real-time streaming update mode is activated. In the event-triggered update mode, the data sampling interval is dynamically adjusted according to the elevation deviation change rate. When the elevation deviation change rate exceeds a preset rate threshold, the data sampling interval is shortened and the data output frequency of the GNSS-RTK receiver is increased. 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 parallel computing nodes to generate processed data blocks; The processed data blocks are spatially stitched and timestamp aligned to output real-time terrain mapping results with spatiotemporal consistency markers.

5. The method according to claim 4, characterized in that, The calculation of the absolute value of the elevation deviation between the three-dimensional terrain surface data and the corresponding location in the historical terrain mapping results includes: The historical elevation dataset corresponding to the target area is extracted from the historical topographic mapping results. The historical elevation dataset contains the coordinates of elevation points collected at multiple historical time points and their corresponding elevation values. Spatial matching is performed between 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. For each successfully matched coordinate pair, calculate the absolute value of the difference between the current elevation value and the historical elevation value; The weighted average of the absolute values ​​of the differences of all successfully matched coordinate pairs is calculated. The weights are determined based on the terrain stability coefficient of the region where the coordinate points are located. A first weight is set in geologically active areas, and a second weight is set in flat areas. The first weight is less than the second weight. The absolute value of the elevation deviation is obtained by multiplying the weighted average value by the terrain change sensitivity coefficient. The process of determining the terrain change sensitivity coefficient includes: Acquire geological tectonic activity monitoring data of the target area, including surface displacement rate, seismic activity frequency, and groundwater level change parameters; The surface displacement rate is compared with the reference displacement rate to generate a displacement deviation index; A sensitivity adjustment factor is selected based on the level range of the earthquake activity frequency, which includes low frequency range, medium frequency range and high frequency range; Hydrological influencing factors are calculated based on the groundwater level change parameters, and the hydrological influencing factors are positively correlated with the magnitude of water level change. The displacement deviation index, sensitivity adjustment factor, and hydrological influence factor are normalized to generate a comprehensive influence coefficient. The terrain change sensitivity coefficient is generated by performing a nonlinear transformation on the comprehensive influence coefficient using an exponential function.

6. The method according to claim 5, characterized in that, The step of independently verifying and interpolating each data block using parallel computing nodes includes: Each data block is assigned a unique spatial identifier code, which includes information on longitude range, latitude range, and elevation range. A preset terrain data verification rule base is loaded into each parallel computing node. The rule base includes rules for detecting elevation abrupt changes, rules for contour line continuity, and rules for slope consistency. A streaming data processing engine is used to verify the elevation points in the data block point by point. When any elevation point is detected to violate the elevation change detection rule, the elevation point is marked as an anomaly and an interpolation compensation mechanism is triggered. For locations marked as outliers, interpolation compensation is performed using the moving average of surrounding elevation points to generate compensated elevation points. The verified and compensated data blocks are sent to the central node, so that the central node can sort and splice the data blocks according to the spatial identifier code; During the stitching process, the elevation difference at the boundary of adjacent data blocks is detected. If the elevation difference exceeds the preset stitching threshold, the boundary smoothing algorithm is activated to re-interpolate the boundary area.

7. The method according to claim 6, characterized in that, The initiated boundary smoothing algorithm re-interpolates the boundary region, including: Extract elevation points within a preset width range at the boundary of adjacent data blocks to generate a set of boundary elevation points; Calculate the average elevation value and standard deviation of elevation for all points in the set of boundary elevation points; Based on the average elevation value, a Gaussian distribution model is constructed to generate the elevation confidence weight for each elevation point; The elevation points in the boundary area are weighted and fused according to the elevation confidence weight to generate the fused boundary elevation value. Bilinear interpolation was used to smooth the transition of the merged boundary elevation values ​​and eliminate splicing gaps. The processed boundary elevation values ​​are updated in adjacent data blocks, and the elevation differences in the spliced ​​area are recalculated.

8. The method according to claim 7, characterized in that, The process of generating the spatiotemporal consistency marker includes: Add a data acquisition timestamp and 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. A processing pipeline number is generated based on the processing order of the data blocks, and the number includes a node identifier and a data processing stage code; The timestamp, spatial coordinate range identifier, and processing pipeline number are combined and encoded to generate a globally unique spatiotemporal consistency marker. During the output of real-time terrain mapping results, the spatiotemporal consistency marker is embedded in the metadata segment of the data file and hash-bound with the three-dimensional terrain surface data.

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

Images