High-precision DEM modeling and feature line extraction method based on second-order variational model

Abstract

The invention belongs to the technical field of DEM modeling, and discloses a high-precision DEM modeling and feature line extraction method based on a second-order variational model. The method comprises the following steps: firstly, establishing a second-order variation model containing a fidelity item, a smooth item and a discontinuity detection item according to characteristics of a fracture terrain, then solving the second-order variation model by adopting a vector block coordinate descent algorithm based on finite difference, and realizing high-precision DEM modeling while maintaining the characteristics of the fracture terrain. Finally, the fracture feature points recognized by the second-order variation model are converted into feature lines through operations such as clustering denoising and contraction refining, and therefore accurate extraction of the topographic feature lines is achieved. The method well solves the problem that a traditional curved surface modeling method cannot effectively maintain fracture characteristics for complex terrains.

Description

Technical Field

[0001] The present invention belongs to the technical field of DEM modeling, and specifically relates to a high-precision DEM modeling and feature line extraction method based on a second-order variational model. Background Art

[0002] The Digital Elevation Model (DEM), as the spatio-temporal basis of the real-scene three-dimensional map, provides reliable basic data for geological disaster prevention and control, national territorial space planning, etc. During the construction of DEM, traditional spatial interpolation methods are based on the basic assumption of surface continuity, and are prone to losing discontinuous fracture terrain features such as cliffs, ridges, and slopes, resulting in DEM distortion. Therefore, how to effectively maintain terrain features in complex areas is the key to achieving high-precision DEM modeling.

[0003] In response to this, researchers have proposed a variety of methods, which can be roughly classified into two categories: one is spatial interpolation based on known fracture line constraints. For example, a feature-preserving DEM is constructed, which accurately depicts complex fracture terrain features by embedding irregular triangular meshes near the fracture lines and dynamically adjusting the layout and shape of the triangles. Or, a conditional generative adversarial network is used, and key feature lines are extracted from a coarse-resolution open-source terrain dataset as generation conditions to construct a deep learning model to generate high-precision DEM. Another example is a high-precision DEM interpolation method constrained by terrain features. First, an unconstrained TIN model is established, and then an edge detection algorithm is used to extract fracture feature lines, which are added to the TIN model as constraint conditions. Although the above methods show great advantages in fracture terrain modeling, they still face the challenge of efficiently and accurately extracting terrain fracture lines from ground point data.

[0004] The other category is to construct DEM using a spatial interpolation method that takes into account terrain features when the fracture lines are unknown. For example, a weighted radial basis function interpolation method integrating structure tensors effectively enhances the DEM's ability to reproduce complex terrain features (such as cliffs, ridges, etc.) by mining the gradient vector information of sampling points around the fracture lines. A multi-radial basis function interpolation method that fuses spatial distance, elevation difference, and normal vector information realizes high-precision characterization of fracture terrain areas by considering the spatial dependence relationship between sampling points and points to be interpolated and the heterogeneity of terrain features. Although this type of method shows good performance in capturing and retaining specific fracture terrains (such as crease and jump fracture regions), its effect is still insufficient when dealing with fracture terrains with complex curvature changes such as slopes.

[0005] Therefore, it is necessary to propose a high-precision DEM modeling and feature line extraction method based on a second-order variational model to solve the above technical problems existing in the prior art. Summary of the Invention

[0006] The object of the present invention is to provide a high-precision DEM modeling and feature line extraction method based on a second-order variational model, so as to solve the problem that traditional surface modeling methods cannot effectively preserve fracture features for complex terrains and achieve accurate extraction of terrain feature lines.

[0007] To achieve the above object, the present invention provides the following technical solutions:

[0008] The high-precision DEM modeling and feature line extraction method based on a second-order variational model includes the following steps:

[0009] Step 1. Establish a second-order variational model including a fidelity term, a smoothness term, and a discontinuity detection term according to the characteristics of fractured terrains.

[0010] Step 2. Use a vector block coordinate descent algorithm based on finite differences to minimize the second-order variational model to achieve DEM modeling and identification of fracture feature points.

[0011] Step 3. Convert the fracture feature points identified by the second-order variational model into feature lines through clustering denoising, shrinkage refinement, feature point connection, and line smoothing to achieve extraction of terrain feature lines.

[0012] Compared with the prior art, the present invention has the following beneficial effects:

[0013] As described above, the high-precision DEM modeling and feature line extraction method based on a second-order variational model described in the present invention uses the multi-feature constraint property of the second-order variational model to preserve fracture terrain features such as cliffs, ridges, and slopes, thereby achieving high-precision modeling in complex terrain areas; in addition, the discontinuity detection term in the second-order variational model is used to identify feature points, and through operations such as clustering denoising and shrinkage refinement, they are converted into feature lines, achieving high-quality extraction of terrain feature lines. The present invention solves the problem that traditional surface modeling methods cannot effectively preserve fracture features for complex terrains. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required to be used in the embodiments.

[0015] Figure 1 It is a flowchart of the high-precision DEM modeling and feature line extraction method based on a second-order variational model of the present invention;

[0016] Figure 2 It is a schematic diagram of the geometric process of second-order variational minimization in an embodiment of the present invention;

[0017] Figure 3 It is a slope change diagram at a fracture feature;

[0018] Figure 4 Schematic diagram of the process for feature line extraction

[0019] Figure 5 Reference DEM of 6 groups of data in the embodiment of the present invention

[0020] Figure 6 Schematic diagram of comparison of RMSE and MAE when each method in the embodiment of the present invention processes 6 groups of data

[0021] Figure 7 DEM hill shade map constructed by each method in the embodiment of the present invention

[0022] Figure 8 Topographic feature line map extracted by each method in the embodiment of the present invention Detailed implementation manners

[0023] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0024] Embodiment

[0025] As Figures 1 to 8 shown, this embodiment describes a high-precision DEM modeling and feature line extraction method based on a second-order variational model. The method includes the following steps:

[0026] Step 1. Establish a second-order variational model including a fidelity term, a smoothness term, and a discontinuity detection term according to the characteristics of fractured terrain.

[0027] Let the topographic surface function be expressed as: z i = f(x i , y i ) + e i , where e i is the surface modeling error, and f(·) is the mapping function from plane coordinates to elevation. Therefore, the key to surface modeling is to accurately calculate this mapping function. Generally speaking, traditional smooth surface modeling methods achieve surface simulation through an effective compromise between fidelity and smoothness, and the objective function can be expressed as:

[0028]

[0029] In the formula, λ is the smoothing parameter, is the fidelity term, is the smoothness term.

[0030] As can be seen from the above formula, the smooth surface modeling method lacks the ability to recognize the geometric features of first-order and second-order discontinuities, resulting in the over-smoothing of fracture features and the distortion of the modeling results, as Figure 2 shown.

[0031] Therefore, in this embodiment, a second-order variational model with the ability to recognize discontinuity features is introduced in surface modeling to avoid the over-smoothing problem in complex terrain areas and thus maintain the fracture terrain features. Specifically, after adding the second-order variational Blake-Zisserman (BZ) model to the above formula, it is expressed as:

[0032]

[0033] In the formula, is the fidelity term, is the smooth term, is the discontinuity detection term, λ is the smoothing parameter, f is the piecewise smooth approximation of z, S f and respectively represent the discontinuity set of f and the discontinuity set of the gradient of f, H 1 represents the Hausdorf measure, Ω represents the region of sampled data, z represents the elevation of the point cloud data, and α, β represent the weight parameters.

[0034] As can be seen from Equation (1), due to the existence of the unknown measurement term H 1 , the functional F(f) is non-differentiable. To address this, this embodiment introduces two auxiliary functions m, n: Ω → [0, 1] (i.e., the indicator function of the discontinuity set), and represents the Γ-convergence approximation of F(f) through a uniformly elliptic functional. The formula is:

[0035]

[0036] In the formula, ε is the Γ-convergence parameter, and α, β, λ 1 , λ 2 are the weight parameters. After the improvement, the functional is differentiable.

[0037] As Figure 2 shown, by analyzing the roles of each function term in Equation (2), the geometric behaviors of the functions m, n, and f in the fracture terrain area can be predicted. Specifically:

[0038] (1) Due to the minimization of the overall functional, the approximation term ||f – z|| 2 forces the function f to approximate the sampled value z. To keep (m - 1) 2 / 4ε bounded, m needs to be 1 throughout the region Ω. However, when the sampling point is at a jump fracture, such as Figure 3 a, the constraint term The value will increase significantly, resulting in an increase in the value of the target functional. To maintain the minimization of the functional, m is changed from 1 to 0. At the same time, The constraint effect of is heterogeneous, causing the function f to perfectly fit the sampled value z, avoiding the smoothing of the discontinuous features in this area, and thus effectively retaining the jump fracture terrain features. In addition, during the minimization process, keeping

[0039] (2) Similarly, when the point cloud is located at the curvature and crease fracture points, such as Figure 3 b and Figure 3 c, the value will increase significantly. To ensure the minimization of the objective functional, the value of n is changed from 1 to 0. This process also eliminates the constraint effect, thereby retaining the curvature and crease fracture features.

[0040] As can be seen from the above, f is an approximation of z under Γ - convergence, while m and n are approximations of the discontinuous sets S f and respectively. In addition, the points where the auxiliary functions m and n are equal to 0 represent the first - order and second - order discontinuous points of f, that is, the fracture feature points.

[0041] Step 2. Since the non - convex property of the second - order variational model makes its minimization solution extremely difficult, in this embodiment, a vector block coordinate descent algorithm based on finite differences is used to solve the minimization of equation (2) to achieve DEM modeling and fracture feature point recognition.

[0042] Step 2.1. Discretize equation (2) using grid - based discretization techniques and approximate the first - order and second - order differential operators in the form of finite differences, and then represent equation (2) as a matrix and decompose it into a multi - variable block form.

[0043] Assume the vector is represented as v, then Rv represents a diagonal matrix with diagonal elements equal to v, and v 2 is the vector of the square coefficients of v, that is, [v 2 i =([v] i ) 2 . Therefore, the differential operator at each grid point can be approximated as:

[0044]

[0045] where v ij represents the vector at (i, j), w(i, j)=(j - 1)N + i, and i and j respectively represent the horizontal and vertical coordinates after the point cloud data is gridded, i = 1,…,N, j = 1,…,M; D x v, D y ​v, D xx v, D yy v, D xy v is a first-order or second-order differential operator in the finite difference format.

[0046] Therefore, Equation (2) can be expanded as follows:

[0047]

[0048] where A m = A m (f), A n = A n (f), A f = A f (m, n).

[0049] A m 、A n 、A f 、b m 、b n 、b f can be expressed as follows respectively:

[0050]

[0051] where e = (1, 1,..., 1) T , I is the identity matrix, and b = 1,..., B.

[0052] Step 2.2. Equation (4) is a structure composed of multi-variable blocks. Using the Gauss-Seidel method (GS method) to solve it variable by variable, the variable blocks m, n, and f can be expressed as follows respectively:

[0053]

[0054] Among them, for the function f(x), argmin f(x) represents the value of the variable x when the function f(x) reaches the minimum value; in this embodiment, m k+1 = argmin m F(m, n k , f k ) means that when F(m, n k , f k ) reaches the minimum value, the value of the variable m is equal to m k+1 .

[0055] The functional is quadratic with respect to each variable block, and the descent solution can be performed for each variable block respectively (Steps 2.2.1 to 2.2.4). In order to find a suitable gradient-related search direction d k m 、d k n and dk fb , the present invention uses the preconditioned conjugate gradient method to perform iterative solutions for linear and .

[0056] The process of performing a descent solution for each variable block is as follows:

[0057] Step 2.2.1. Define the initial values of m k and n k as m 0 = n 0 = (1, 1,..., 1) T , The initial value of 0 is f m = z, the step size parameter γ n = γ f = 1, the step size parameter γ

[0058] Step 2.2.2. Calculate the search directions d k m and d k n ;

[0059] Calculate

[0060] Iteratively calculate

[0061] Calculate

[0062] Iteratively calculate

[0063] Step 2.2.3. Calculate the search direction

[0064] Calculate

[0065] Iteratively calculate

[0066] Step 2.2.4. After k iterations, determine whether the relative change of the functional satisfies the convergence condition, and its formula is:

[0067]

[0068] In the formula, tol F is the preset threshold;

[0069] If the convergence condition is satisfied, stop the iterative calculation; otherwise, let k = k + 1, and return to Step 2.2.2 until convergence.

[0070] During the minimization process, functions m and n are 0 at the breaks and 1 in other regions; while f is considered as a smoothed constraint approximation of z, that is, the output result that avoids over - smoothing to retain fracture terrain features such as cliffs, ridges, and slopes. The method in this embodiment can stop after k iterations, and the threshold tol F is taken as 0.01m.

[0071] Step 3. To extract accurate and complete terrain feature lines, the fracture feature points identified by the second - order variational model are transformed into feature lines through clustering denoising, shrinking and thinning, feature point connection, and line smoothing, so as to achieve the accurate extraction of terrain feature lines, as Figure 4 shown.

[0072] First, fracture feature points are extracted by functions m and n in Equation (2) to form a data set( Figure 4 a);

[0073] Secondly, the DBSCAN algorithm is used to cluster the fracture feature points, specifying the search radius r and the minimum number of samples MINpts to determine the cluster to which each point belongs and generate cluster labels; after clustering, calculate the distance from each point to the center of its cluster, and regard the points whose distance from the cluster center exceeds the set threshold r (search radius r) as noise points and remove them from the data set( Figure 4 b);

[0074] Then, the Laplacian - based skeleton shrinking method is used to thin the data set( Figure 4 c);

[0075] Finally, the Minimum Spanning Tree (MST) is used to connect the remaining fracture feature points in the data set and form terrain feature lines. To make the lines smoother, a spline function is used to fit them to generate smooth curves( Figure 4 d).

[0076] In addition, to prove the effectiveness of the method of the present invention, the following experiments are also given.

[0077] First, 6 groups of ground sampling points are selected from the airborne LiDAR point cloud data set provided by the open - source website https: / / portal.opentopography.org as the research objects, which include landform types such as mountains, hills, and river networks, and there are various fracture terrain features, as Figure 5 shown.

[0078] Then, each group of point cloud data is randomly divided into training data and verification data at a ratio of 9:1. Among them, the training data is used to construct the DEM, and the verification data is used for the quantitative analysis of the modeling. Since there may be some outlier points in the acquired ground point cloud, in order to ensure the quality of the ground point cloud data, it is necessary to preprocess the data. During the data preprocessing, first use Terrascan software to perform initial filtering on the point cloud, and then manually edit to remove the misclassified points, and finally obtain accurate ground points. Table 1 lists the statistical information such as the elevation, average slope, number of training points, and number of verification points of 6 groups of point cloud data.

[0079] Table 1 Statistical Information of 6 Groups of Point Cloud Data

[0080]

[0081] To quantitatively analyze the efficiency of the modeling of the present invention, it is compared with IDW, RBF, OK, constrained TIN, and M-RBF (Multivariate Radial Basis Function). In order to ensure that the performance of each method reaches the optimal, the cross-validation technique is used to obtain the optimal parameters of each method. Two indicators, root mean square error (RMSE) and mean absolute error (MAE), are used to verify the accuracy, and the formulas of each accuracy index are as follows:

[0082]

[0083] Where M i * is the predicted value, M i is the reference value, and n is the number of verification points.

[0084] In addition, this embodiment uses three indicators, completeness, correctness, and quality, to verify the accuracy of extracting the terrain feature lines of the present invention, and compares it with D8 and TSAE (Two-Step Adaptive Extraction). The formulas of each accuracy index are as follows:

[0085]

[0086] Where TP represents the number of feature points that actually exist, that is, the number of feature points that actually exist and are extracted; FN represents the number of feature points that are not extracted, that is, the number of feature points that actually exist but are not extracted; FP represents the number of feature points that are wrongly extracted, that is, the number of feature points that do not actually exist but are extracted.

[0087] Next, a comparison of the DEM modeling is carried out.

[0088] Combined with Table 1 andFigure 6 It can be seen that the complexity of terrain features has a great impact on the accuracy of DEM. Specifically, the accuracy of DEM decreases as the terrain complexity increases. For example, Data2, Data5, etc. have low modeling accuracy due to the inclusion of complex terrain features such as steep banks, ridges, and slopes. On the contrary, Data1 and Data3 have flat surfaces or high point cloud densities, so their modeling accuracies are high. From Figure 6 It can also be known that the accuracy of IDW is always the worst under various terrains, because the predicted values of IDW cannot exceed the maximum and minimum values of the sampling point values, and thus large errors are likely to occur in the sparse point cloud areas selected in this embodiment.

[0089] In comparison, the method of the present invention always has the best modeling accuracy in all terrains. The average RMSE is reduced by 36.6%, 21.1%, and 19.6% respectively compared with the traditional methods IDW, RBF, and OK, and the average MAE is reduced by 39.8%, 26.1%, and 24.4% respectively; compared with the better-performing M-RBF and constrained TIN, the average RMSE is reduced by 9.1% and 10.9% respectively, and the average MAE is reduced by 15.6% and 14.5% respectively.

[0090] In order to further compare the modeling performances of the 6 methods in the fracture area, local point cloud areas containing crease fractures ( Figure 5 b), jump fractures ( Figure 5 d), and curvature fractures ( Figure 5 f) are intercepted from Data2, Data4, and Data6 respectively for modeling comparative analysis (Table 2, Figure 8 ).

[0091] Due to the complex and polymorphic terrain and sparse point clouds around the fracture features, the modeling accuracy of each method in this area drops significantly compared with that in the overall area. From Table 2 and Figure 7 it can be seen that the reduction of RBF is the largest, indicating that the modeling performance of this method in the fracture area is the worst. On the contrary, the reduction of the method of the present invention is the smallest.

[0092] Table 2 also shows the average RMSE and MAE of each modeling method when dealing with local fracture terrains. It can be known from this that the average RMSE of the method of the present invention is reduced by 38.3%, 40.7%, and 31.9% respectively compared with the traditional methods IDW, RBF, and OK, and the average MAE is reduced by 43.8%, 47%, and 40.3% respectively. Compared with the terrain feature-considering M-RBF and constrained TIN, the average RMSE is reduced by 17.3% and 11.8% respectively, and the average MAE is reduced by 25.9% and 18.5% respectively.

[0093] Table 2 RMSE and MAE (m) of each method when dealing with the local fracture area

[0094]

[0095]

[0096] To more intuitively compare the modeling effects of six methods, 3D mountain shadow maps of the local fracture area were generated using each method ( Figure 7 ).

[0097] As Figure 7 can be seen, the visualization effect of the method of the present invention is the best and the terrain features are maintained the best ( Figure 7 f), while the other methods exhibit different defects. Specifically, the visualization effect of IDW is the worst, and there are a large number of obvious discontinuities in the interpolation surface ( Figure 7 a); the RBF surface is extremely rough and accompanied by many small pits, and this phenomenon is mostly seen around the crease fracture features ( Figure 7 b); the image texture obtained by OK can well reflect the real ground surface, but at the crease and curvature fracture features, there are many fine cracks ( Figure 7 c); M-RBF is prone to pseudo-features at the curvature fracture terrain ( Figure 7 d); the constrained TIN can better maintain various fracture features, but a small amount of terrain details are lost ( Figure 7 e).

[0098] The following is a comparison of the feature line extraction.

[0099] Table 3 shows the accuracy comparison of the feature line extraction of Data2, Data4, and Data6 by three methods.

[0100] The results show that the extraction accuracy of the method of the present invention is the highest in all terrains. Among them, the accuracy and integrity of the feature line extraction by the TSAE and D8 algorithms in Data2 differ greatly. Since this data is mostly rugged mountain terrain, there is an easy phenomenon of incorrect connection of feature lines; for Data4, the extraction accuracy of the method of the present invention is as high as 85.2%, but the precision is 70.3%. This is mainly because the ground surface of Data4 has large undulations and the point cloud density at the fracture is small, resulting in incomplete extraction of the feature lines at some fractures; in Data6, the extraction integrity of TSAE is poor (45.3%), attributed to its principle based on plane intersection, resulting in the inability to identify the feature lines at the elevation continuous but curvature discontinuous places. In contrast, the extraction integrity of D8 has been greatly improved (66.3%), but the accuracy is low (42.1%), which is due to the extraction of a large number of pseudo-feature lines. The method of the present invention better balances the above two indicators, and the precision has been greatly improved.

[0101] Table 3 Results and Precision of Feature Line Extraction by Each Method

[0102]

[0103]

[0104] The terrain feature lines extracted by each method are as Figure 8 shown. The results show that the extraction effect of TSAE is the worst ( Figure 8 a), many feature lines are not extracted and the trend of the line body is chaotic; the feature lines extracted by the D8 algorithm are more comprehensive ( Figure 8 b), but there are many pseudo feature lines and the connection of the line body is disordered; in contrast, the feature lines extracted by the method of the present invention are significantly superior to the other two methods in terms of accuracy, integrity, and line body continuity ( Figure 8 c).

[0105] Aiming at the problem that the traditional surface modeling method cannot effectively maintain the fracture features of complex terrain, the present invention proposes a high-precision DEM modeling and feature line extraction method based on a second-order variational model, which maintains the fracture terrain features such as cliffs, ridges, and slopes through the multi-feature constraint property of the second-order variational model, so as to realize high-precision modeling in complex terrain areas. In addition, the discontinuity detection term in the second-order variational model is used to identify feature points, and they are converted into feature lines through operations such as clustering denoising and shrinking refinement to achieve high-quality extraction of terrain feature lines. With the help of 6 groups of data containing various fracture features, the method of the present invention is compared with 5 other interpolation methods (IDW, RBF, OK, M-RBF, constrained TIN) for modeling. The results show that the method of the present invention has the best accuracy, the average RMSE is at least 19.6% lower than that of the classical interpolation algorithm, and the average MAE is at least 24.4% lower, and it can better maintain the terrain details in the fracture area. In addition, taking Data2, Data4, and Data6 as examples, by comparing the feature lines extracted by the method of the present invention with those extracted by D8 and TSAE, it can be seen that the accuracy of the method of the present invention is at least increased by 36.5%, and the extracted terrain feature lines have a smooth trend and good integrity, which are more in line with the actual terrain features.

[0106] The embodiments of the present invention are only used to illustrate the technical solutions of the present invention rather than to limit them. For those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirits of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A high-precision DEM modeling and feature line extraction method based on a second-order variational model, characterized in that: The steps include: Step 1. Establish a second-order variational model including fidelity term, smooth term and discontinuity detection term according to the characteristics of the fault terrain; Step 2. Use the vector block coordinate descent algorithm based on finite difference to minimize the second-order variational model to achieve DEM modeling and fracture feature point identification; Step 3. The fracture feature points identified by the second-order variational model are converted into feature lines through clustering denoising, shrinkage and refinement, feature point connection and line smoothing to achieve the extraction of terrain feature lines.

2. The high-precision DEM modeling and feature line extraction method according to claim 1 is characterized in that: In the step 1, based on the traditional smooth surface modeling, a second-order variational model with discontinuity feature recognition capability is introduced to avoid the over-smoothing problem in complex terrain areas and to maintain the fracture terrain characteristics.

3. The high-precision DEM modeling and feature line extraction method according to claim 2 is characterized in that: In step 1, based on the smooth surface modeling, the second-order variational Blake-Zisserman model is introduced, and its objective function F(f) is expressed as: In the formula, For authenticity, is a smooth term, is the discontinuity detection term, λ is the smoothing parameter, f is the piecewise smoothing approximation of z, S f and They represent the discontinuity set of f and the discontinuity set of the gradient of f, respectively. 1 represents the Hausdorf measure, Ω represents the area of ​​the sampling data, z represents the elevation of the point cloud data, and α and β represent weight parameters.

4. The high-precision DEM modeling and feature line extraction method according to claim 3 is characterized in that: In step 1, for the unknown measurement item H 1 There is a problem that causes the functional F(f) to be non-differentiable. Two auxiliary functions m and n are introduced: Ω→[0,1], and the Γ-convergent approximation of F(f) is expressed by a uniform elliptic functional. The formula is: Wherein, ε is the Γ-convergence parameter, α, β, λ1, λ2 are weight parameters, and after improvement, the functional F(f) is differentiable.

5. The high-precision DEM modeling and feature line extraction method according to claim 4 is characterized in that: The step 2 is specifically as follows: Step 2.

1. First, discretize equation (2) based on the grid discretization technology, and use the finite difference form to approximate the first-order and second-order differential operators, then express equation (2) as a matrix and decompose it into a multivariable block form; Step 2.

2. Implement descent solution for each variable block based on Gauss-Seidel method.

6. The high-precision DEM modeling and feature line extraction method according to claim 5, characterized in that: The step 2.1 is specifically as follows: First, assuming that the vector is denoted by v, then Rv represents a diagonal matrix whose diagonal elements are equal to v, v 2 is the vector of v squared coefficients, i.e. [v 2 ] i =([v] i ) 2 , so the differential operator at each grid point is expressed as: In the formula, v ij represents the vector at (i, j), w(i, j) = (j-1)N+i, i and j represent the horizontal and vertical coordinates of the gridded point cloud data, i = 1, ..., N, j = 1, ..., M; D x v、D y v、D xx v、D yy v、D xy v is the first-order and second-order differential operators in the finite difference format; Then, equation (2) is expressed as a matrix and decomposed into a multivariate block form: In the formula, A m , b m , A n , b n , A f , b f Respectively expressed as: Where, e=(1,1,…,1) T , I is the unit matrix, b=1,…,B.

7. The high-precision DEM modeling and feature line extraction method according to claim 6 is characterized in that: The step 2.2 is specifically as follows: The Gauss-Seidel method is used to solve equation (4) variable by variable, and each variable block m, n, f is expressed as:

8. The high-precision DEM modeling and feature line extraction method according to claim 7, characterized in that: In step 2.2, the process of descending and solving each variable block is as follows: Step 2.2.

1. Define m k and n k The initial value is m 0 =n 0 =(1, 1, ..., 1) T , The initial value is f 0 =z, step size parameter γ m =γ n =1, step size parameter γ f =1.5, k=0; Step 2.2.

2. Calculate the search direction d k m and d k n ; calculate Iterative Calculation calculate Iterative Calculation Step 2.2.

3. Calculate the search direction d k fb ; calculate Iterative Calculation Step 2.2.

4. Determine whether the relative change of the functional satisfies the convergence condition. If so, stop the iterative calculation; otherwise, set k=k+1 and return to step 2.2.2 until convergence.

9. The high-precision DEM modeling and feature line extraction method according to claim 8, characterized in that: In step 2.2.4, after iterating k times, it is determined whether the relative change of the functional satisfies the convergence condition, and the formula is: In the formula, tol F is the preset threshold.

10. The high-precision DEM modeling and feature line extraction method according to claim 4, characterized in that: The step 3 is specifically as follows: First, the fracture feature points are extracted through the functions m and n in formula (2) to form a data set; Secondly, the DBSCAN algorithm is used to cluster the fracture feature points, specifying the search radius r and the minimum number of samples MINpts to determine the cluster to which each point belongs and generate cluster labels; After clustering is completed, the distance from each point to the center of its cluster is calculated, and the points whose distance from the cluster center exceeds the search radius r are regarded as noise points and removed from the data set; Then, the dataset is refined using Laplacian’s skeleton shrinkage method; Finally, the minimum spanning tree is used to connect the fracture feature points retained in the data set to form terrain feature lines, and then the spline function is used to fit them to generate a smooth curve.

Images