High-precision DEM modeling method for faulted terrain area considering spatial heterogeneity

Abstract

The application discloses a high-precision DEM modeling method for a faulted terrain area considering spatial heterogeneity, and in view of the problem that the existing interpolation algorithm has low interpolation accuracy at discontinuous terrain features such as fault lines, a multi-element radial basis function interpolation method considering spatial heterogeneity is provided. Compared with the traditional interpolation algorithm, the application provides a multi-element radial basis kernel function coupling three kinds of terrain feature information of distance, height difference and normal vector, which weakens the influence of sampling points on the opposite side of the interpolation point on the interpolation result, and thus has strong maintaining ability for fold faulted terrains such as ridges and valleys, and jump faulted terrains such as cliffs and terraces. In addition, the application also provides specific experiments, and the experiments show that compared with the interpolation method constrained by the structure tensor and several classical interpolation methods, the average total error of the application method is the smallest, the interpolation performance is the best, and the faulted terrain features can be well maintained, thereby proving the effectiveness of the application method.

Description

Technical Field

[0001] This invention relates to a high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity. Background Technology

[0002] Currently, high-resolution and high-precision DEMs have become indispensable key supports for major national needs such as refined urban management, natural resource management and disaster monitoring, and the construction of a realistic 3D China. DEMs are mainly obtained by interpolating discrete sampling points using spatial interpolation methods. Due to the complexity and polymorphism of the Earth's surface, classical interpolation methods, based on the assumptions of surface continuity or smoothness, are affected by fault terrain features such as steep slopes, gullies, and terraces. This results in the elevation around the fault lines being smoothed, leading to distortion in the constructed DEM. Therefore, preserving the characteristics of fault terrain is crucial for achieving high-quality DEM modeling.

[0003] To address the above issues, scholars both at home and abroad have proposed various modeling algorithms, which can be broadly categorized into two types.

[0004] The first approach utilizes known fault lines as constraints to assist in DEM interpolation. For example, some literature has proposed a two-step constrained Delaunay triangulation method. The idea is to first generate an unconstrained Delaunay triangulation, and then adjust the triangulation by embedding constrained fault lines to construct a high-precision DEM. Other literature proposes incorporating fault lines into multi-resolution thin-plate spline surface fitting by using the square derivative, thereby constructing a high-quality DEM with fault terrain features.

[0005] However, since there is still a lack of effective technical means to automatically and accurately extract various types of terrain fault lines from raw data, the whole process relies on a high degree of human intervention, resulting in the limited practicality of such methods.

[0006] The second approach fully considers the spatial distribution characteristics of the fracture terrain when the fault line is unknown, directly interpolating the original ground points into the DEM. The key to this method is better perceiving the spatial distribution of terrain features and reducing the influence of sampling points on opposite sides of the interpolation point on the calculation results. One paper proposed a point cloud feature-preserving denoising algorithm suitable for DEM construction, which rationally allocates the weights of sampling points by judging the difference in normal vectors on both sides of the fault line. This method shows good ability in preserving folded fracture terrain such as ridges and valleys; however, it is not effective in preserving abrupt fracture terrain such as cliffs and terraces. Another paper proposed a weighted radial basis function method that considers the structure tensor. By integrating the structure tensor into the radial basis function, it fully utilizes the gradient and direction information of sampling points near the fault line, achieving high-quality DEM modeling for various terrain features. However, this method is significantly affected by the accuracy of the normal vector and easily loses fracture terrain features with large elevation changes.

[0007] In summary, existing interpolation algorithms suffer from low interpolation accuracy in discontinuous terrain features such as fault lines. Summary of the Invention

[0008] The purpose of this invention is to propose a high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity. This method couples three types of terrain information, such as spatial distance, elevation difference, and normal vector, and fully considers the spatial correlation and terrain feature heterogeneity between sampling points and interpolation points to ensure high-precision DEM modeling of fractured terrain areas.

[0009] To achieve the above objectives, the present invention adopts the following technical solution:

[0010] A high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity includes the following steps:

[0011] Step 1. Construct a DEM grid with resolution r using the maximum and minimum coordinates of the ground point cloud;

[0012] Step 2. Calculate the elevation of each point to be interpolated using the nearest neighbor interpolation method;

[0013] Step 3. Calculate the spatial distance, elevation difference, and normal vector between each interpolation point and its neighboring sampling points;

[0014] Step 4. Construct a multivariate kernel function based on spatial distance, elevation difference, and the angle between normal vectors, and substitute the multivariate kernel function into the radial basis function to calculate the elevation of each interpolation point, thereby realizing surface simulation to obtain the DEM;

[0015] Step 5. Determine whether the maximum value of the DEM difference obtained in this iteration and the previous iteration is less than the preset threshold; if not, return to step 3, and recalculate the elevation difference between the interpolation point and the neighboring sampling points and the normal vector of each interpolation point and the neighboring sampling points based on the elevation value of the interpolation point obtained in this iteration; otherwise, output the final interpolation result.

[0016] The present invention has the following advantages:

[0017] As described above, this invention relates to a high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity. The method involves a multivariate radial basis function interpolation method that takes into account spatial heterogeneity. Compared with traditional interpolation algorithms, by proposing a multivariate kernel function that couples three terrain feature information—distance, elevation difference, and normal vector—it effectively reduces the influence of sampling points on opposite sides of the interpolation point on the interpolation result. Therefore, the high-precision DEM modeling method for fractured terrain areas described in this invention has a strong ability to preserve fractured terrain such as ridges and valleys, as well as jumping fractured terrain such as cliffs and terraces. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of a high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity in an embodiment of the present invention.

[0019] Figure 2 This is a schematic diagram of the point cloud distribution at the fault terrain.

[0020] Figure 3 This is a schematic diagram showing the comparison of normal vector estimation.

[0021] Figure 4 3D mountain shadow maps generated by various interpolation methods when processing data s23.

[0022] Figure 5 This diagram illustrates the comparison of RMSE (Recovery Mean Squared Error) for various methods at different point cloud densities when processing surface fault data.

[0023] Figure 6 This is a schematic diagram comparing the MAE of various methods under different point cloud densities when processing surface fault data.

[0024] Figure 7 Mountain shadow maps generated by various interpolation methods when processing surface fault data. Detailed Implementation

[0025] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0026] When modeling a Digital Elevation Model (DEM), traditional interpolation methods assume a smooth and continuous surface and only consider the spatial correlation between sampling points and interpolation points. They ignore the spatial heterogeneity caused by discontinuous terrain features such as fault lines, which results in the elevation around the fault lines being smoothed, leading to distortion of the constructed DEM.

[0027] To address the aforementioned technical issues, this embodiment proposes a high-precision DEM modeling method for fractured terrain areas that takes into account spatial heterogeneity. The method of this invention proposes a multivariate radial basis function interpolation method that takes into account spatial heterogeneity. By coupling three types of terrain information, such as spatial distance, elevation difference, and normal vector, it fully considers the spatial correlation and terrain feature heterogeneity between sampling points and interpolation points, thereby ensuring high-precision DEM modeling of fractured terrain areas.

[0028] like Figure 1 As shown, the high-precision DEM modeling method for fractured terrain areas that take into account spatial heterogeneity includes the following steps:

[0029] Step 1. Construct a DEM grid with resolution r using the maximum and minimum coordinates of the ground point cloud.

[0030] Step 2. Calculate the elevation of each point to be interpolated using the nearest neighbor interpolation method.

[0031] Step 3. Calculate the spatial distance, elevation difference, and normal vector between each interpolation point and its neighboring sampling points.

[0032] Step 4. Construct a multivariate kernel function based on spatial distance, elevation difference, and the angle between normal vectors. Substitute the multivariate kernel function into the radial basis function to form a multivariate radial basis function. Calculate the elevation of each interpolation point to achieve surface simulation and obtain the DEM.

[0033] Based on their properties, fracture lines can be divided into two types: skip fracture lines and crease fracture lines.

[0034] When the point to be interpolated is located around a jump break line (e.g., a cliff, a curb) (e.g.) Figure 2 a) The elevation of its neighboring points changes drastically; similarly, when the point to be interpolated is located around the fold line (e.g., valley, ridge) (e.g.) Figure 2 (b) There are large abrupt changes in the normal vector within its neighborhood. Therefore, if all these neighborhood points are used for the calculation of the interpolation point, it can easily lead to DEM distortion.

[0035] To address the aforementioned problems, this invention constructs a multivariate kernel function Φ that takes into account terrain features by coupling three types of terrain information: distance, elevation difference, and normal vector. The construction process of the multivariate kernel function Φ is as follows:

[0036] Φ=Φ(d ij )×Φ(h ij )×Φ(n ij (1)

[0037] In the formula, Φ(d) ij ), Φ(h ij ), Φ(n ij ) represent spatial distance, elevation difference, and normal vector kernel function, respectively.

[0038]

[0039] In the formula, d ij σ is the Euclidean distance between the point to be interpolated and its neighboring sampling points. d is the kernel radius of the kernel function, used to control the smoothness of the basis functions. Equation (2) shows that the closer the sampling point is to the point to be interpolated, the greater the weight assigned, and thus the greater the influence on the interpolation; conversely, the farther the sampling point is from the point to be interpolated, the smaller the weight assigned, and the smaller the influence on the interpolation.

[0040]

[0041] In the formula, h i h is the elevation value of the i-th point to be interpolated. jLet σ be the elevation value of the j-th neighboring sampling point. h The kernel radius is used to adjust the smoothness of the basis functions.

[0042] Equation (3) shows that the weight assigned to a neighboring sampling point is inversely proportional to the elevation difference between it and the point to be interpolated. That is, the greater the elevation difference between a neighboring sampling point and the point to be interpolated, the smaller the weight assigned to that point, and vice versa.

[0043] Generally, sampling points located on opposite sides of the feature line within the point cloud neighborhood have a significant difference in elevation from their corresponding points. Figure 2 As shown in (a).

[0044]

[0045] n i ·n j =|n i |×|n j |×cosθ (5)

[0046] In the formula, n i Let n be the normal vector of the i-th interpolation point. j Let θ be the normal vector of the j-th neighboring sampling point, and θ be the angle between the normal vectors of the point to be interpolated and the normal vectors of the neighboring sampling points, where 0 ≤ θ ≤ π, and σ n The kernel radius is used to adjust the smoothness of the basis functions.

[0047]

[0048] In the formula, k is the number of neighboring sampling points of the i-th interpolation point, and N(i) is the set of neighboring sampling points of the i-th interpolation point.

[0049] Equation (4) shows that the larger the angle between the normal vector of the sampling point and the normal vector of the interpolation point, the smaller the weight assigned; conversely, the smaller the angle between the normal vector of the sampling point and the normal vector of the interpolation point, the larger the weight assigned.

[0050] Generally, sampling points on the same side of the feature line within the neighborhood of a point cloud have a smaller angle with it, while those on opposite sides have a larger angle. Figure 2 As shown in (b).

[0051] As can be seen from equation (4), the accuracy of point cloud normal vector estimation directly affects the final interpolation accuracy.

[0052] Conventional normal vector estimation methods treat the point cloud neighborhood as a smooth surface, assigning equal weights to points within the neighborhood for normal vector estimation. However, they neglect the influence of anisotropy in point cloud distribution under special conditions such as fractured terrain, leading to inaccurate normal vector calculations at terrain features, such as... Figure 3 As shown in (a).

[0053] This invention utilizes an iterative weighted PCA method based on the M-estimator, employing dual automatic initialization to address the anisotropy problem near fractured terrain features, thereby improving the accuracy of normal vector estimation at different terrain features. Figure 3 As shown in (b).

[0054] Let the estimated normal vector of any point p0 in the grid be n, and its residual relative to its neighboring points be r. i (n)=n·(p i –p0). Where p i This represents the normal vector of the neighboring sampling points.

[0055] The choice of the estimator kernel function is crucial to ensuring good convergence of the optimization. The chosen kernel function is:

[0056]

[0057] In the formula, μ is a quantitative parameter that decreases with each iteration to control the number of iterations. Newton's method is used to solve the optimization problem, followed by IRLS minimization to obtain the weight w for each point. i (n):

[0058]

[0059] The normal vector n for each point is obtained by combining the classic weighted PCA method. * for:

[0060] n * =argmin∑ i w i (n).r i (n) 2 (III)

[0061] First, we perform the initialization step by calculating the initial value n at p0 using classic PCA. init1 Then, the two steps of equation (II) and equation (III) are executed alternately: the weights are updated relative to the currently estimated normal vector, and then a new normal vector estimate is calculated by weighted PCA. With continuous iteration, the direction of the normal vector is gradually optimized, and finally the normal vector n1 is obtained.

[0062] Then, the second initialization step is performed. When p0 is near the break line, the normal vector n1 resulting from the first step may be attracted to the side with more points / or the side with less curvature. In this case, it is necessary to use n init2 =n1^v1(v1:n init1 The cross product between n1 and n2 is used to set the initial normal n in the second step. init2 Repeat the iterative process in the first step to finally obtain the normal vector n2.

[0063] After obtaining two normal vectors n1 and n2, select the normal vector of the face closest to the point as the final normal vector.

[0064] The stopping of the iteration is controlled by a quantitative parameter μ.

[0065] The multivariate kernel function in this invention is a nonlinear kernel function that defines the contribution of neighboring points to the interpolation point by coupling the spatial distance, elevation difference, and normal vector angle between the sampling point and the interpolation point. For example, in a smooth and continuous flat region, the elevation and normal vectors of points within the point cloud neighborhood are not significantly different; in this case, the multivariate kernel function approximates the result of a regular Gaussian kernel function. When the interpolation point is located around abruptly fractured terrain such as cliffs or terraces, the elevation changes drastically while the normal vector remains relatively constant; in this case, the multivariate kernel function can be considered as a kernel function coupling distance and elevation. Similarly, in fractured terrain regions such as ridges and valleys, the normal vector undergoes abrupt changes while the elevation remains continuous; in this case, the multivariate kernel function can be considered as a kernel function coupling distance and normal vector. Therefore, the improved multivariate kernel function of this invention has good universality for various terrain features.

[0066] For the designed multivariate kernel function Φ, the process of multivariate radial basis function interpolation calculation is given below:

[0067] For any point x to be interpolated, there are n neighboring sampling points. If interpolation is involved, then the RBF interpolation model is:

[0068]

[0069] In the formula, and λ i These are the RBF kernel function and the weights, p k (x) and b i These are the polynomial function and its corresponding coefficients.

[0070] The matrix form of equation (7) is as follows:

[0071] f=Φλ+Pb (8)

[0072] The multivariate kernel function Φ obtained in formula (1) is used as the Φ value in formula (8).

[0073]

[0074] Therefore, when using RBF for interpolation calculations, the coefficients λ and b must first be calculated.

[0075] Using the biased estimation principle of ridge regression, the objective function of RBF is obtained as follows:

[0076]

[0077] In the formula, ρ is a smoothing parameter, the value of which can be obtained through cross-validation.

[0078] Let F = (Φλ + Pb - f) T (Φλ+Pb-f)+ρλ T Taking the derivative of F with respect to λ, we get:

[0079]

[0080] That is: (Φ+ρI)λ+Pb=f (11)

[0081] Taking the derivative of F with respect to b, we get:

[0082]

[0083] Substituting equation (11) into equation (12) yields:

[0084] P T λ=0 (13)

[0085] From equations (11) and (13), we obtain:

[0086]

[0087] The coefficients of the RBF model are obtained by solving equation (14).

[0088] Substituting the RBF model coefficients into equation (8) estimates the elevation of the point to be inserted, thereby enabling surface simulation to obtain the DEM.

[0089] Step 5. Determine whether the maximum value of the DEM difference obtained in this iteration and the previous iteration is less than the preset threshold; if not, return to step 3, and recalculate the elevation difference between the interpolation point and the neighboring sampling points and the normal vector of each interpolation point and the neighboring sampling points based on the elevation value of the interpolation point obtained in this iteration; otherwise, output the final interpolation result.

[0090] Inaccurate initial values ​​of the points to be interpolated lead to deviations in the estimation of normal vectors, resulting in inaccurate interpolation results. To achieve high-precision DEM modeling, this invention employs a multivariate radial basis function stepwise update strategy to calculate the elevation of the points to be interpolated.

[0091] For example, let the current iteration number be the i-th iteration, where i ≥ 2, such as... Figure 1 As shown.

[0092] Define the result of the i-th interpolation as DME i Define the interpolation result of the (i-1)th iteration as DME i-1The maximum value of the DEM difference obtained in this (i-th) iteration and the previous iteration is defined as max(DEM). i -DME i-1 ).

[0093] If max(DME) i -DME i-1 If ) < tol, then the final interpolation result will be output.

[0094] Otherwise, if max(DME) i -DME i-1 If )≥tol, based on the elevation values ​​obtained in this iteration (the i-th iteration result), recalculate the elevation difference between the insertion point and the sampling point, as well as the normal vector between each insertion point and the sampling point.

[0095] This invention iterates through the calculation of elevation difference, estimation of normal vector, and interpolation until the interpolation result tends to stabilize.

[0096] In addition, the following experiments are presented to demonstrate the effectiveness of the method of the present invention.

[0097] First, 10 sets of airborne LiDAR point cloud benchmark datasets provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) were selected as the research objects. These datasets contain rich terrain features, enabling comprehensive accuracy verification of the method presented in this invention.

[0098] For each public dataset, all ground points are randomly divided into training points and check points in a ratio of 9:1. Training points are used to construct the DEM, and check points are used to evaluate the accuracy of the DEM. Table 1 lists in detail the number of training points and check points, survey area, point cloud density, elevation standard deviation, and other information for the 10 datasets.

[0099] Table 1. Statistical data of the study area

[0100]

[0101] In addition, this invention also selected a set of complex surface fault point cloud data from Wellington, New Zealand for accuracy verification. This dataset has obvious fault characteristics and can be used to test the ability of various interpolation methods to preserve terrain features.

[0102] The study area covers 25 hectares and has an average slope of 17.4°.

[0103] To ensure the quality of the experimental data, Terrascan software was first used to automatically obtain the initial filtering results of the point cloud. Then, the filtered ground point cloud was manually checked and corrected, and finally a total of 384,906 ground points in the survey area were obtained.

[0104] The elevation range of this ground point is [101.6, 211.2m], the standard deviation of elevation is 22.8m, and the average distance between points is 0.48m.

[0105] To verify the DEM accuracy of the interpolation method under different sampling densities, the training dataset was randomly reduced to 70%, 50%, and 30% of the full density (90% of the total sample data). To verify the efficiency of the method of this invention, its calculation results were compared with the accuracy of the structure tensor constrained interpolation method (ST-RBF) and four classical interpolation methods (including inverse distance weighted interpolation (IDW), ordinary kriging (OK), ANUDEM, and standard radial basis interpolation (RBF)).

[0106] The optimal parameters for each method were obtained through cross-validation.

[0107] The selected accuracy metrics include Mean Square Error (RMSE) and Mean Absolute Error (MAE), calculated using the following formulas:

[0108]

[0109]

[0110] In the formula: Z i * Z is the predicted value of the known sampling point i. i is the reference value for sampling point i, and n is the number of checkpoints.

[0111] The calculation results for ISPRS data (Table 2) show that all methods have the lowest and highest accuracy in s53 and s31, respectively. This is mainly because the surface of s53 is rough and the point cloud density is low, while the surface of s31 is flat (Table 1). In addition, except for s21 and s31 (these two sets of data are mostly flat areas), the DEM generated by the new method has better accuracy than the other five methods.

[0112] Table 2 shows the RMSE and MAE (m) of the new method compared to five other methods for 10 sets of ISPRS data.

[0113]

[0114] Table 2 also shows the average RMSE and MAE of the method of this invention compared to five other methods when processing ten sets of ISPRS data. It can be seen that ANUDEM performs the worst under both accuracy metrics because it employs a global smoothness assumption during interpolation. In contrast, the method of this invention achieves the lowest average RMSE (0.341m) and MAE (0.129m). Specifically, the average RMSE of the method of this invention is reduced by 16.0%, 17.2%, 19.8%, and 14.8% compared to the classic interpolation methods IDW, OK, ANUDEM, and RBF, respectively, and the average MAE is reduced by 24.6%, 24.1%, 25.0%, and 24.1%, respectively. Compared to ST-RBF, the new method reduces the average RMSE by 6.1% and the average MAE by 7.2%. The higher interpolation accuracy of the method of this invention is due to the consideration of the spatial heterogeneity of discontinuous terrain features such as break lines during the interpolation process.

[0115] Furthermore, this invention also compares the mountain shadow map after s23 interpolation using six different methods, such as... Figure 4 As shown.

[0116] The results show that IDW has obvious pseudo-pits (such as...). Figure 4 (circular symbol in a) OK (e.g., circle) Figure 4 b) circular mark) and standard RBF (such as Figure 4 The circular and rectangular markers in section d exhibit some missing terrain features and have a certain smoothing effect on fractured terrain. This is particularly noticeable in the OK section. Figure 4 The circular marker in b) is mainly due to the fact that the second-order stationarity assumption used by OK is difficult to hold in fault terrain areas. Among the six methods, ANUDEM has the worst visualization effect, with an abnormally smooth surface and a large number of irregular protrusions. Figure 4 (circular marker in c). Compared to classic interpolation methods, ST-RBF (such as...) Figure 4 (As shown in e) shows a significant improvement in overall interpolation performance; however, it is not as effective as the interpolation method proposed in this invention (e.g., in areas with abrupt changes in fractured terrain) Figure 4 (e.g., oval and rectangular icons in f).

[0117] Figure 5 and Figure 6The diagrams show a comparison of RMSE and MAE for various methods under different point cloud densities. As can be seen, the accuracy of each interpolation method decreases with decreasing point cloud density. Therefore, the interpolation error of each algorithm is negatively correlated with point cloud density. Furthermore, IDW consistently exhibits the worst interpolation accuracy because its estimates cannot exceed the maximum and minimum values ​​of the sampled points, leading to significant errors at lower point cloud densities. OK and AUNDEM show similar interpolation accuracy and are less affected by point cloud density. Moreover, regardless of the point cloud density, the interpolation accuracy of the method described in this invention remains optimal.

[0118] Table 3 lists the average interpolation errors (RMSE and MAE) for various interpolation methods. The results show that IDW has the worst accuracy among the classical interpolation methods, while OK has the highest accuracy. The interpolation accuracy of the method in this invention is significantly improved compared to the five interpolation methods. Its average RMSE is reduced by 46.9%, 41.8%, 41.1%, 43.7%, and 6.7% compared to IDW, OK, ANUDEM, standard RBF, and ST-RBF, respectively, and its average MAE is reduced by 61.0%, 58.3%, 58.6%, 59.9%, and 8.3%, respectively.

[0119] Table 3 shows the average RMSE and MAE (m) of the new method compared to five other methods when processing surface fault data.

[0120]

[0121] Furthermore, this invention also compares various interpolation methods for processing DEM hillside shadow maps (point cloud density of 90%) from surface fault data, such as... Figure 7 As shown. The results show that the limitation of IDW is the presence of abnormal protrusions on the hillside (such as...). Figure 7 (e.g., oval icon in a) OK (e.g., ...) Figure 7 (as shown in b) and ANUDEM (as shown in b) Figure 7 As shown in c), it is very smooth, losing a lot of terrain detail. RBF retains terrain detail better than the previous two, but has many unreasonable pseudo-dimples (such as...). Figure 7 (d) Circular markers. In comparison, the method of this invention can better preserve the characteristics of the fractured terrain (such as...). Figure 7 (Medium-sized rectangular icon).

[0122] Using 10 sets of public data provided by ISPRS and 1 set of surface fault airborne lidar point cloud data, a comparison of the method of this invention with five other interpolation methods (IDW, OK, ANUDEM, standard RBF, and ST-RBF) shows that the method of this invention has the highest accuracy, with an average RMSE reduction of at least 14.8% and an average MAE reduction of at least 24.1% compared to the four classic interpolation algorithms. Furthermore, the method of this invention significantly outperforms the other five methods in terms of preserving terrain details and fault locations.

[0123] Of course, the above description is only a preferred embodiment of the present invention. The present invention is not limited to the above-described embodiments. It should be noted that any equivalent substitutions or obvious modifications made by those skilled in the art under the guidance of this specification fall within the scope of this specification and should be protected by the present invention.

Claims

1. A high-precision DEM modeling method for fractured terrain areas considering spatial heterogeneity, characterized in that, Includes the following steps: Step 1. Construct a DEM grid with resolution r using the maximum and minimum coordinates of the ground point cloud; Step 2. Calculate the elevation of each point to be interpolated using the nearest neighbor interpolation method; Step 3. Calculate the spatial distance, elevation difference, and normal vector between each interpolation point and its neighboring sampling points; Step 4. Construct a multivariate kernel function based on spatial distance, elevation difference, and the angle between normal vectors, and substitute the multivariate kernel function into the radial basis function to calculate the elevation of each interpolation point, thereby realizing surface simulation to obtain the DEM; Step 5. Determine whether the maximum value of the DEM difference obtained in this iteration and the previous iteration is less than the preset threshold; if not, return to step 3, and recalculate the elevation difference between the interpolation point and the neighboring sampling points, as well as the normal vector between each interpolation point and the neighboring sampling points, based on the elevation value of the interpolation point obtained in this iteration; otherwise, output the final interpolation result. Multivariate kernel function The construction process is as follows: (1) In the formula, , , These represent spatial distance, elevation difference, and normal vector kernel function, respectively. (2) In the formula, d ij The Euclidean distance between the point to be interpolated and its neighboring sampling points is given. The kernel radius is used to adjust the smoothness of the basis functions; (3) In the formula, h i h is the elevation value of the i-th point to be interpolated. j Let σ be the elevation value of the j-th neighboring sampling point. h The kernel radius is used to adjust the smoothness of the basis functions; (4) (5) In the formula, n i Let n be the normal vector of the i-th interpolation point. j Let j be the normal vector of the j-th neighboring sampling point. The angle between the normal vectors of the point to be interpolated and the neighboring sampling points. , The kernel radius is used to adjust the smoothness of the basis functions; (6) In the formula, k is the first... The number of neighboring sampling points of each interpolation point. Let be the set of neighboring sampling points of the i-th interpolation point.

2. The high-precision DEM modeling method for fractured terrain areas according to claim 1, characterized in that, In step 4, the radial basis function (RBF) interpolation calculation process is as follows: For any point x to be interpolated, there are n neighboring sampling points. If interpolation is involved, then the RBF interpolation model is: (7) In the formula, and These are the RBF kernel function and the weights, respectively. and These are the polynomial function and its corresponding coefficients; The matrix form of equation (7) is as follows: (8) The multivariate kernel function obtained in formula (1) As in formula (8) value; ； ； ； ； ； Therefore, when using RBF for interpolation calculations, the coefficients must first be calculated. and ; Using the biased estimation principle of ridge regression, the objective function of RBF is: (9) In the formula, For smoothing parameters; set up ,right Seeking information about The derivative is obtained as follows: (10) Right now: (11) Taking the derivative of F with respect to b, we get: (12) Substituting equation (11) into equation (12), we get: (13) From equations (11) and (13), we obtain: (14) The RBF model coefficients are obtained by solving equation (14), and the elevation of the insertion point is estimated by substituting the coefficients into equation (8).

3. The high-precision DEM modeling method for fractured terrain areas according to claim 1, characterized in that, In step 3, the normal vectors of the point to be inserted and the neighboring sampling points are calculated using the iterative weighted PCA method based on the M-estimator.

Images