TIN_DDM terrain feature line automatic extraction method based on rolling ball transformation

Abstract

The application provides a TIN_DDM terrain feature line automatic extraction method based on a rolling ball transformation, and belongs to the technical field of terrain feature line extraction and application of seabed terrain. The application can overcome the problems in traditional terrain feature line extraction, such as difficulty in quantitative regulation and control of a threshold, incapability of self-adaptive adjustment of a connection mode and incompleteness of terrain feature line types, realizes complete, continuous and subdivided automatic extraction and multi-scale expression of a TIN-DDM terrain feature line, and the extracted terrain feature line has higher terrain reconstruction precision. The application can quickly realize construction of a buffer surface under any buffer radius, and the overall buffer surface construction efficiency is far greater than that of the same type of invention, and with the increase of the number of sampling points, the efficiency can be improved by nearly 40 times.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of terrain feature line extraction and application of seabed terrain, and relates to a TIN_DDM terrain feature line automatic extraction method based on rolling ball transformation. BACKGROUND

[0002] Terrain feature line is a structured framework of Digital Depth Model (DDM) and is a link between local terrain feature analysis and overall terrain form control. Compared with DDM which directly expresses terrain undulation form, terrain feature line contains relatively abstract and implicit terrain structure information and is usually indirectly obtained from DDM. As a basic framework and form reference of DDM terrain feature analysis, terrain feature line is widely used in DDM data compression, terrain reconstruction, geomorphology identification, hydrological analysis and the like.

[0003] At present, the terrain feature line extraction research method is relatively in-depth for GRID_DEM and mainly includes image processing, curved surface geometry analysis, and surface water simulation. Although the image processing method can realize rapid determination of GRID_DEM terrain feature points, it cannot effectively provide a connection mode between sampling points. The curved surface geometry analysis method also has the characteristics of disordered connection of terrain feature points, and the terrain feature line extraction is incomplete. The surface water simulation method has the characteristics of expansion from local to overall, realizes automatic determination and adaptive connection of terrain feature points according to the accumulation threshold of catchment, and is widely used in DEM watershed and confluence line extraction engineering practice. However, the method is affected by the accumulation threshold of catchment, and the terrain feature line (ridge line and valley line) extraction result is not unique, and the terrain feature line extracted by the method lacks the extraction of slope break line, and it is difficult to construct a complete DEM terrain structured framework. Therefore, the method is relatively limited in the application fields of DEM data compression, terrain reconstruction, and geomorphology identification. In addition, compared with GRID_DEM which is processed by interpolation, TIN_DDM uses original sampling points to construct seabed terrain, the terrain feature points are more accurate, and the terrain feature line is more real, and is more suitable for terrain feature line extraction. SUMMARY

[0004] In order to overcome the limitations that the threshold value cannot be adaptively adjusted and the terrain feature line connection cannot be automatically adjusted in GRID_DEM terrain feature line extraction, and further solve the difficulty of complete TIN_DDM terrain feature line extraction, the application provides a TIN_DDM terrain feature line adaptive extraction method based on rolling ball transformation, which can realize complete, subdivided and continuous terrain feature line extraction and multi-scale expression.

[0005] The technical scheme adopted by the application to achieve the above purpose is:

[0006] An automatic extraction method for TIN_DDM terrain feature lines based on rolling ball transform includes the following steps:

[0007] a. Read the sampling point set Λ={P1,P2…P} within TIN_DDM on the computer. i …P n The coordinate data of the sampling points in |i∈(1,n)} are obtained, and the P of each sampling point is calculated. i The positive critical rolling ball radius r i1 And the negative critical rolling ball radius r i2 ;

[0008] b. According to the negative critical rolling ball radius r i2 The sampling points are arranged in descending order of their values ​​to form a new sampling point set Λ. NEW ={P1',P'2…P i '…P' n |i∈(1,n)} and its corresponding new positive critical rolling ball radius r i1 ', New negative critical rolling ball radius r i2 ';

[0009] c. Based on the sampling point set Λ NEW Construct a Delaunay triangulation and solve for P at each sampling point. i 'Natural neighborhood point set Ω0(P i ') and store;

[0010] d. Set the sampling point set Λ NEW Each sampling point P i 'The set of index numbers of the subdivision unit Index(P) i '), whose initial state is an empty set;

[0011] e. Combining sampling point P i 'Natural neighborhood point set Ω0(P i ') and the set of index numbers Index(P i Based on the terrain boundary point determination criteria, the new sampling point set Λ is determined sequentially. NEW Construct a set of terrain boundary points Ψ = {P1”, P2”…P4”} for the inner terrain boundary points. k …P” m |k∈(1,m)}, and set the index of the sampling points within the terrain boundary point set Ψ to the corresponding index set Index(P”). k Updated to {k}; where m represents the new sampling point set Λ NEW The number of inland terrain boundary points;

[0012] f. Combining sampling point P i 'Natural neighborhood point set Ω0(P i') and the set of index numbers Index(P i Based on the characteristic attributes of the sampling points, and according to the new sampling point set Λ NEW Sequentially determine sampling points P, excluding terrain boundary points. i 'Belonging to the subdivision unit Ω(P)' k And update its index set Index(P) i ');

[0013] g. Divide the same subdivision unit Ω(P”) k ) boundary line Ω B (P) k The connection of these lines forms the boundary line network Γ, which is the positive terrain feature line network extracted based on the critical rolling sphere radius. Pos-Ball ;

[0014] h. In the positive terrain feature line network Γ Pos-Ball Based on this, and according to the multi-scale expression conditions of spatial scale, a positive terrain feature line Γ extracted based on the critical rolling sphere radius is constructed under the spatial scale threshold R condition. Pos-Ball-R ;

[0015] i. Change the order of steps b from descending to ascending, then replace the negative critical rolling sphere radius of step bh with the water depth value, repeat step bh, and use the result as the ridgeline Γ under the spatial scale threshold R condition. Ridge-R ;

[0016] j. Positive terrain feature lines Γ extracted based on the critical rolling sphere radius under the spatial scale threshold R. Pos-Ball-R Ridge lines extracted from water depth values ​​Γ Ridge-R Based on this, and according to the subdivision conditions of terrain feature lines, the ridgeline Γ is calculated under the spatial scale threshold R. Ridge-R With the positive slope broken line Γ Pos-Slope-R .

[0017] k. Replace the negative critical rolling sphere radius in step bh with the positive critical rolling sphere radius, run step bh again, and use the result as the negative terrain feature line Γ under the spatial scale threshold R condition. Neg-Ball-R Replace the negative critical rolling sphere radius from step bh with the water depth value, repeat step bh, and use the result as the valley line Γ under the spatial scale threshold R condition. Valley-R ;

[0018] l. Negative terrain feature lines Γ extracted based on the critical rolling sphere radius under the spatial scale threshold R. Neg-Ball-R Valley lines Γ extracted from water depth numerical values Valley-R Based on the terrain feature line subdivision conditions in step j, the valley line Γ under the spatial scale threshold R is calculated. Valley-RWith negative slope broken line Γ Neg-Slope-R .

[0019] In step a, the positive critical rolling ball radius r i1 Solution process (negative critical rolling ball radius r) i2 The solution is similar:

[0020] 1) Set the maximum buffer distance r max Select Λ = {P1, P2, ..., P} in sequence. i …P n Sampling point P within |i∈(1,n)} i Selected sampling point P i As the initial reference sampling point P ij (j=1).

[0021] 2) Based on sampling point P i The formula for the distance difference between the intersection points along the z-axis is used to solve the point set Ω. i The specific process is as follows:

[0022] Select the sampling point set Λ={P1,P2…P} in sequence. i …P n Sampling point P within |i∈(1,n)} t Calculate point P t With point P ij An equidistant sphere at point P i The distance difference between the intersection points along the z-axis is calculated using the following formula:

[0023]

[0024] In the formula: Sampling point P t Reference sampling point P ij z-coordinate, r max To preset the maximum rolling ball radius, Sampling point P t Reference sampling point P i1 To sampling point P i The horizontal distance. If Δ t If >0, then the sampling point P will be... t Place the point set Ω i .

[0025] 3) Calculate the point set Ω using the formula for the radius of the intersection point. i Internal sampling point P t With sampling point P i The radius of intersection when the intersection points coincide t,ij The formula for calculating the radius of the intersection point is:

[0026]

[0027] 4) Select the minimum intersection overlap radius, denoted as r min , and its corresponding sampling point, denoted as P min .

[0028] 5) If there exists a minimum intersection overlap radius r min and its corresponding sampling point P min , and the intersection overlap radius r min is greater than the last critical radius (the initial critical radius r i1 is 0), then record it in the data chain of the selected sampling point in TIN_DDM, that is, the next key sampling point P ij+1 = P min , and its corresponding critical radius r ij+1 = r min .

[0029] 6) Take P ij+1 as the reference sampling point of the next stage, repeat steps 2)-5), until the minimum intersection overlap radius of this stage is less than the critical radius of the last stage or there is no minimum intersection overlap radius and corresponding sampling point, then select the next sampling point in the point set Λ in turn, and when all the sampling points are solved, the forward index data chain as shown in Figure 2 can be constructed, wherein P1 to P n are each sampling point in TIN_DDM, j1 to j n are the number of reference sampling points corresponding to each sampling point, is a series of key sampling points of sampling point P i , is the critical radius corresponding to the series of key sampling points of sampling point P i .

[0030] Based on the forward index data chain, take the maximum value of the critical radius r ij corresponding to the reference sampling point P i when P ij as the forward critical rolling ball radius r i of the sampling point P i1 , and the forward critical rolling ball radius of other sampling points is solved in the same way.

[0031] In step e, the natural neighborhood point set Ω0(P i ') of each sampling point is known, and the terrain boundary point judgment condition is:

[0032]

[0033] In the formula, P' y is the natural neighborhood Ω0(P i ') of the sampling point Pi Sampling points within ') These represent sampling points P respectively. i '、P' y The positive critical rolling sphere radius. This formula indicates that if the sampling point P i The positive critical rolling sphere radius of ' in the natural neighborhood Ω0(P) i If the value within ') is the largest, then P i ' is the positive terrain boundary point.

[0034] In step f, the natural neighborhood set Ω0(P) of each sampling point is known. i If '), then the condition for determining the feature attributes of the sampling points is:

[0035]

[0036] In the formula: point P' y For sampling point P i 'Natural Neighborhood Ω0(P i The formula represents: if P is a sampling point within a certain range, then P is a sampling point within a certain range. i 'Natural Neighborhood Ω0(P i All sampling points P' within ') y Both belong to the same subdivision unit Ω(P”) k ), that is, each sampling point P' y Index set Index(P' y If the intersection of the two sets of indexes is {k}, then the set of indexes Index(P) is... i Let ') be {k}, and determine P. i ' is a non-terrain feature point; if sampling point P i 'Natural Neighborhood Ω0(P i Sampling point P' within ') y Not belonging to the same subdivision unit Ω(P) k ), that is, each sampling point P' y Index set Index(P' y If the intersection of the two samples is empty, then the sampling point P is empty. i ' represents the set of index numbers for terrain feature points, Index(P) i ') is Ω0(P i Sampling point P' within ') y Index set Index(P' y The union of ).

[0037] In step h, the condition for multi-scale spatial scale expression is:

[0038]

[0039] Where: Φ positive (Φ negativerepresents the convex (concave) terrain range constructed by the convex (concave) terrain point and the convex (concave) terrain nested point under the spatial scale threshold R; Γ Pos-Ball-R represents the positive terrain feature line network Γ positive within the convex terrain range Φ Pos-Ball ; Γ Neg-Ball-R represents the negative terrain feature line network Γ negative within the concave terrain range Φ Neg-Ball . Wherein, the convex (concave) terrain point and the convex (concave) terrain nested point determination condition is:

[0040]

[0041] In the formula: r i1 ' represents the positive critical rolling ball radius; r i2 ' represents the negative critical rolling ball radius; Q (P i ') represents the terrain type attribute of the sampling point P i ', wherein: 1 represents the convex terrain point; -1 represents the concave terrain point; 0 represents the flat terrain point; 2 represents the concave terrain point within the convex nested; -2 represents the convex terrain point within the concave nested. By the adjacency clustering of the convex terrain point and the convex terrain nested point in TIN_DDM, the convex terrain range Φ positive is constructed; the adjacency clustering of the negative terrain point and the negative terrain nested point in TIN_DDM, the concave terrain range Φ negative is constructed.

[0042] In the step j, the positive slope break line Γ Pos-Slope-R under the spatial scale threshold R determination condition is:

[0043]

[0044] In the formula: Γ Pos-Slope-R , Γ Neg-Slope-R respectively represent the positive terrain slope break line and the negative terrain slope break line under the spatial scale threshold R. Formula (7) represents: under the spatial scale threshold R, the terrain feature line Γ Feature-R is composed of four lines of ridge line Γ Ridge-R , valley line Γ Valley-R , positive terrain slope break line Γ Pos-Slope-R and negative terrain slope break line Γ Neg-Slope-R .

[0045] The present application can overcome the problems of traditional terrain feature line extraction, such as the difficulty in quantitative regulation of threshold value, the incapability of self-adaptive adjustment of connection mode, and the incompleteness of terrain feature line types, realize complete, continuous and subdivided TIN-DDM terrain feature line automatic extraction and multi-scale expression, and the extracted terrain feature line has higher terrain reconstruction precision. BRIEF DESCRIPTION OF DRAWINGS

[0046] Figure 1 is a flow chart of the TIN_DDM terrain feature line automatic extraction method based on rolling ball transformation.

[0047] Figure 2 is a schematic diagram of forward index data.

[0048] Figure 3 is a schematic diagram of the natural neighborhood point set of a sampling point.

[0049] Figure 4 is a schematic diagram of terrain boundary point determination.

[0050] Figure 5 is a schematic diagram of sampling point feature attribute determination.

[0051] Figure 6 is a schematic diagram of the forward terrain feature line network extracted according to the critical rolling ball radius.

[0052] Figure 7(a) is a sampling point terrain type attribute diagram.

[0053] Figure 7(b) is a schematic diagram of the forward terrain feature line extracted according to the critical rolling ball radius under the condition of a spatial scale threshold R = 600 m.

[0054] Figure 8(a) is a schematic diagram of the ridge line and the forward slope break line in a high area under the condition of a spatial scale threshold R = 600 m.

[0055] Figure 8(b) is a schematic diagram of the ridge line and the forward slope break line in a gully area under the condition of a spatial scale threshold R = 600 m. DETAILED DESCRIPTION

[0056] The present application proposes a TIN_DDM terrain feature line automatic extraction method based on rolling ball transformation, which realizes the determination of the subdivision unit to which the sampling point belongs and the connection of the subdivision unit boundary by using a computer. For example, Figure 3 As shown in Figure 8, assuming that the existing TIN_DDM sampling point data, taking the extraction of forward terrain feature lines (including forward slope break lines and ridge lines) as an example, includes the following steps (as shown in Figure 1 ​

[0057] Step a. Read the set of sampling points Λ = {P1, P2, ..., P3} within TIN_DDM on the computer. i …P n The coordinate data of the sampling points in |i∈(1,n)} are used to calculate the positive critical rolling sphere radius r of each sampling point. i1 And the negative critical rolling ball radius r i2 ;

[0058] Step b. According to the negative critical rolling ball radius r i2 The sampling points are arranged in descending order of their values ​​to form a new sampling point set Λ. NEW ={P1',P'2…P i '…P' n |i∈(1,n)} and its corresponding new positive critical rolling ball radius r i1 ', New negative critical rolling ball radius r i2 ';

[0059] Step c. Based on the sampling point set Λ NEW Construct a Delaunay triangulation and solve for P at each sampling point. i 'Natural neighborhood point set Ω0(P i And store, specifically as follows Figure 3 As shown, the dashed triangular network represents the Delaunay triangular network, and the black circle represents any sampling point in the sampling point set Λ, denoted by P. Point P and the sampling points A, B, C, D, E, F, and G represented by the black dashed circles sharing the Delaunay triangle sides constitute the natural neighborhood set Ω0(P) of point P.

[0060] Step d. Set the new sampling point set Λ NEW Each sampling point P i 'The set of index numbers of the subdivision unit Index(P) i '), whose initial state is an empty set;

[0061] Step e. Combine the natural neighborhood point set Ω0(P) of the sampling point i ') and the set of index numbers Index(P i Based on the terrain boundary point determination criteria, determine the new sampling point set Λ. NEW Construct a set of terrain boundary points Ψ = {P1”, P2”…P4”} for the inner terrain boundary points. k …P” m |k∈(1,m)}, for the set of index numbers corresponding to the sampling points within the terrain boundary point set Ψ, Index(P) i The ') is updated to {k}, and the specific judgment result is as follows: Figure 4As shown, the black hexagon represents the selected terrain boundary point, and its negative critical rolling sphere radius is the maximum value of its natural neighborhood point set;

[0062] Step f. Combining the natural neighborhood point set of the sampling points, and based on the characteristic attribute judgment conditions of the sampling points, according to the new sampling point set Λ NEW Sequentially determine the sampling point P outside the terrain boundary point i 'Belonging to the subdivision unit Ω(P)' k And update the index set Index(P) i The judgment result is as follows: Figure 5 As shown, black squares represent non-positive terrain feature points, and black dots represent positive terrain feature points.

[0063] Step g. Divide the same subdivision unit Ω(P”) k The boundary line Ω B (P) k The connection results in a positive terrain feature line network Γ. Pos-Ball ;like Figure 6 As shown, the black solid lines represent the positive terrain feature line network, and each black solid line grid is a subdivision unit.

[0064] Step h. In the forward terrain feature line network Γ Pos-Ball Based on this, and according to the conditions for multi-scale spatial scale expression, a positive terrain feature line Γ is constructed under the condition that the spatial scale threshold is R. Pos-Ball-R Figure 7(a) shows the terrain type attribute map of the sampling points, where convex terrain points and nested convex terrain points are represented by black triangles, and negative terrain points and nested negative terrain points are represented by black rhombuses; Figure 7(b) shows the positive terrain feature line Γ extracted based on the critical rolling sphere radius under the spatial scale threshold R = 600m. Pos-Ball-R (This includes positive slope break lines and some ridge lines).

[0065] Step i. Change the order of steps b from descending to ascending, and replace the negative critical rolling sphere radius of step bh with the water depth value. Repeat step bh above, and use the result as the ridgeline Γ under the spatial scale threshold R condition. Ridge-R ;

[0066] Step j. Positive terrain feature line Γ under the condition of spatial scale threshold R. Pos-Ball-R With the ridgeline Γ Ridge-R Based on this, and according to the terrain feature line subdivision conditions, the ridgeline Γ is calculated under the condition of spatial scale threshold R. Ridge-R With the positive slope broken line Γ Pos-Slope-RFigures 8(a) and 8(b) are schematic diagrams of ridge lines and positive slope break lines for different regional spatial scale thresholds R = 600 m, in which the ridge lines are represented by solid black lines and the positive slope break lines are represented by dashed black lines.

[0067] The present application is described by way of example in the specification. Those skilled in the art will recognize that changes and modifications can be made to the described features and examples without departing from the spirit and scope of the application. Further, those skilled in the art will recognize that modifications can be made to adapt a particular situation or material to the teachings of the application without departing from the scope of the application. Therefore, other embodiments are intended to be within the scope of the application.

Claims

1. A TIN_DDM terrain feature line automatic extraction method based on rolling ball transformation, characterized in that, The method comprises the following steps: a. Read the sampling point set Λ = {P1, P2…Pn} in TIN_DDM on the computer, and calculate the forward critical rolling ball radius r i …P n n of each sampling point P i i∈(1,n)} and the negative critical rolling ball radius r i1 i∈(1,n)} of each sampling point P i2 ​ b. According to the negative critical rolling ball radius r i2 The numerical values are arranged in descending order to form a new set of sampling points Λ NEW = {P1', P2'…P i '…P n '|i∈(1,n)} and its corresponding new positive critical rolling ball radius r i1 ' and the new negative critical rolling ball radius r i2 ' c. According to the sampling point set Λ NEW Construct Delaunay triangulation, solve each sampling point P i Natural neighborhood point set Ω0(P i ) and store; d. Set Λ of sampling points NEW Each sampling point P i Index(P i ), which is initially empty. e.Combining the sampling points P i The natural neighborhood point set Ω0(P i ) and the index number set Index(P i ), according to the terrain boundary point judgment condition, sequentially judge the new sampling point set Λ NEW The terrain boundary point, construct the terrain boundary point set Ψ = {P1”, P2”…P k ”…P m ”|k ∈ (1, m)}, and update the index number set Index(P k ”) corresponding to the sampling point in the terrain boundary point set Ψ to {k}; wherein m represents the number of terrain boundary points in the new sampling point set Λ NEW ​ f. Combining sampling points P i The natural neighborhood point set Ω0(P i ) and the index number set Index(P i ), according to the sampling point characteristic attribute judgment condition, sequentially judge the new sampling point set Λ NEW , except for the terrain boundary point, the sampling point P i belongs to the subdivision unit Ω(P k ) and updates its index number set Index(P i ). g. Connect the boundary lines Ω k (P B ”) of the same dissected unit Ω k ”, and the connected boundary lines are the positive terrain feature line network Γ Pos-Ball extracted according to the critical rolling ball radius. h.In the positive terrain feature line network Γ Pos-Ball On this basis, according to the spatial scale multi-scale expression condition, the positive terrain feature line Γ Pos-Ball-R extracted by the critical rolling ball radius under the spatial scale threshold R condition is constructed i. The order of the arrangement of step b is changed from large to small to small to large, and the negative critical rolling ball radius of step b-h is replaced by the water depth value. Repeat the above steps b-h, and the result is the ridge line Γ of the mountain under the spatial scale threshold R condition Ridge-R ; j. positive terrain feature line Γ extracted by critical rolling ball radius under spatial scale threshold R condition Pos-Ball-R ridge line Γ extracted by water depth value Ridge-R On this basis, according to the terrain feature line subdivision condition, the ridge line Γ under the spatial scale threshold R condition is calculated Ridge-R positive slope break line Γ Pos-Slope-R ; k. Replace the negative critical rolling ball radius of steps b-h with a positive critical rolling ball radius, run steps b-h, and plot the results as a negative topographic feature line Γ for a spatial scale threshold R Neg-Ball-R ; Replace the negative critical rolling ball radius of steps b-h with a water depth value, repeat steps b-h above, and plot the results as a valley line Γ for a spatial scale threshold R Valley-R ; l. Negative topographic feature line Γ extracted with critical rolling ball radius under spatial scale threshold R condition Neg-Ball-R Valley line Γ extracted with water depth value Valley-R Based on the valley line Γ, according to the topographic feature line subdivision condition of step j, the valley line Γ under the spatial scale threshold R condition is calculated Valley-R Negative slope break line Γ Neg-Slope-R .

2. The method of claim 1, wherein, In step a, the forward critical rolling ball radius r i1 The solution process is: 1) Set the maximum buffer distance r max , sequentially select sampling points P i …P n | i ∈ (1, n)} in Λ = {P1, P2…P i n} ; select the sampling point P i as the initial reference sampling point P ij , j = 1; 2) According to the distance difference formula between the intersection points in the z-axis direction of the sampling point P i , the point set Ω i is solved; the specific process is as follows: Select the sampling point set Λ={P1,P2…P} in sequence. i …P n Sampling point P within |i∈(1,n)} t Calculate point P t With point P ij An equidistant sphere at point P i The distance difference between the intersection points along the z-axis is calculated using the following formula: In the formula: respectively are the z coordinates of the sampling point P t , the reference sampling point P ij , and r max is a preset maximum rolling ball radius, respectively are the horizontal distances from the sampling point P t , the reference sampling point P i1 to the sampling point P i ; if Δ t > 0, the sampling point P t is put into the point set Ω i ; 3) Calculate the intersection point set Ω according to the intersection radius solving formula i The inner sampling point P t The intersection point P i The intersection radius r when the intersection point is coincident t,ij The intersection radius calculation formula is: 4) Select the smallest intersection overlap radius, denoted as r min , and its corresponding sampling point, denoted as P min ; 5) if the minimum intersection coincidence radius r min and the corresponding sampling point P min and the intersection coincidence radius r min is greater than the last critical radius, where the initial critical radius r i1 is 0, then it is recorded in the data chain of the selected sampling point in the TIN_DDM, i.e. the next key sampling point P ij+1 = P min with the corresponding critical radius r ij+1 = r min ; 6) P ij+1 As the next stage of the reference sampling point, repeat steps 2) -5), until the minimum intersection point coincidence radius of the stage is less than the critical radius of the previous stage or there is no minimum intersection point coincidence radius and the corresponding sampling point, then in turn select the next sampling point in the point set Λ, when all the sampling points are calculated, the forward index data chain can be constructed. Based on this forward index data chain, a reference sampling point P is selected in the forward index data chain. ij For P i The critical radius r corresponding to time ij The maximum value is taken as the sampling point P. i The positive critical rolling ball radius r i1 The calculation of the positive critical rolling sphere radius for the remaining sampling points is similar.

3. The method according to claim 1 or 2, characterized in that, In step e, the natural neighborhood point set Ω0(P i ) of each sampling point P is known, and the terrain boundary point determination condition is: In the formula, P y 'For sampling point P' i 'Natural Neighborhood Ω0(P i Sampling points within ') These represent sampling points P respectively. i '、P y The positive critical rolling ball radius; This equation means that if the sampling point P i ' has the maximum positive critical rolling ball radius within the natural neighborhood Ω0(P i '), then P i ' is a positive terrain boundary point.

4. The method according to claim 1 or 2, characterized in that, In step f, the natural neighborhood point set Ω0(P i ) of each sampling point is known, and the sampling point feature attribute determination condition is: where P is a sampling point y Ω0(P i ) is the natural neighborhood of P i Ω0(P i ) is the natural neighborhood of P i Ω0(P y ) is the natural neighborhood of P k Ω0(P y ) is the natural neighborhood of P y Ω0(P i ) is the natural neighborhood of P i Ω0(P i ) is the natural neighborhood of P i Ω0(P y ) is the natural neighborhood of P k Ω0(P y ) is the natural neighborhood of P y Ω0(P i ) is the natural neighborhood of P i Ω0(P i ) is the natural neighborhood of P y Ω0(P y ) is the natural neighborhood of P 5. The method of claim 3, wherein, In step f, the natural neighborhood point set Ω0(P i ) of each sampling point is known, and the sampling point feature attribute determination condition is: In the formula: point P y 'For sampling point P i 'Natural Neighborhood Ω0(P i The formula represents: if P is a sampling point within a certain range, then P is a sampling point within a certain range. i 'Natural Neighborhood Ω0(P i All sampling points P within ') y 'Belongs to the same subdivision unit Ω(P) k “), that is, each sampling point P y 'Index set Index(P y If the intersection of the two sets of indexes is {k}, then the set of indexes Index(P) is... i Let ') be {k}, and determine P. i ' is a non-terrain feature point; if sampling point P i 'Natural Neighborhood Ω0(P i Sampling point P within ') y 'Not belonging to the same subdivision unit Ω(P) k “), that is, each sampling point P y 'Index set Index(P y If the intersection of the two samples is empty, then the sampling point P is empty. i ' represents the set of index numbers for terrain feature points, Index(P) i ') is Ω0(P i Sampling point P within ') y 'Index set Index(P y The union of ').

6. The method of claim 1, 2, or 5, wherein, In the step h, the spatial scale multi-scale expression condition is: In the formula: Φ positive represents the convex terrain range constructed by the convex terrain point and the convex terrain nested point adjacent clustering under the spatial scale threshold R condition; Γ Pos-Ball-R represents the positive terrain feature line network Γ positive located in the convex terrain range Φ Pos-Ball ; Γ Neg-Ball-R represents the negative terrain feature line network Γ negative located in the concave terrain range Φ Neg-Ball ; wherein the convex terrain point and the convex terrain nested point and the concave terrain point and the concave terrain nested point determination condition is: wherein: r i1 represents the positive critical rolling ball radius; r i2 represents the negative critical rolling ball radius; Q(P i ) represents the terrain type attribute of the sampling point P i , wherein: 1 represents a convex terrain point; -1 represents a concave terrain point; 0 represents a flat terrain point; 2 represents a convexly nested concave terrain point; -2 represents a concavely nested convex terrain point. Constructing convex terrain range Φ by clustering the adjacent points of convex terrain points and convex terrain nest points in TIN_DDM positive ; constructing concave terrain range Φ by clustering the adjacent points of negative terrain points and negative terrain nest points in TIN_DDM negative .

7. The method of claim 3, wherein, In the step h, the spatial scale multi-scale expression condition is: In the formula: Φ positive represents the convex terrain range constructed by the convex terrain point and the convex terrain nested point adjacent clustering under the spatial scale threshold R condition; Γ Pos-Ball-R represents the positive terrain feature line network Γ positive located in the convex terrain range Φ Pos-Ball ; Γ Neg-Ball-R represents the negative terrain feature line network Γ negative located in the concave terrain range Φ Neg-Ball ; wherein the convex terrain point and the convex terrain nested point and the concave terrain point and the concave terrain nested point determination condition is: wherein: r i1 represents the positive critical rolling ball radius; r i2 represents the negative critical rolling ball radius; Q(P i ) represents the terrain type attribute of the sampling point P i , wherein: 1 represents a convex terrain point; -1 represents a concave terrain point; 0 represents a flat terrain point; 2 represents a convex-type nested inner concave terrain point; -2 represents a concave-type nested inner convex terrain point. Constructing convex terrain range Φ by clustering the adjacent points of convex terrain points and convex terrain nest points in TIN_DDM positive ; constructing concave terrain range Φ by clustering the adjacent points of negative terrain points and negative terrain nest points in TIN_DDM negative .

8. The method of claim 4, wherein, In the step h, the spatial scale multi-scale expression condition is: In the step h, the spatial scale multi-scale expression condition is: In the formula: Φ positive represents the convex terrain range constructed by the convex terrain point and the convex terrain nested point adjacent clustering under the spatial scale threshold R condition; Γ Pos-Ball-R represents the positive terrain feature line network Γ positive located in the convex terrain range Φ Pos-Ball ; Γ Neg-Ball-R represents the negative terrain feature line network Γ negative located in the concave terrain range Φ Neg-Ball ; wherein the convex terrain point and the convex terrain nested point and the concave terrain point and the concave terrain nested point determination condition is: wherein: r i1 represents the positive critical rolling ball radius; r i2 represents the negative critical rolling ball radius; Q(P i ) represents the terrain type attribute of the sampling point P i , wherein: 1 represents a convex terrain point; -1 represents a concave terrain point; 0 represents a flat terrain point; 2 represents a convex nested-in-concave terrain point; -2 represents a concave nested-in-convex terrain point. Constructing convex terrain range Φ by clustering the adjacent points of convex terrain points and convex terrain nest points in TIN_DDM positive ; constructing concave terrain range Φ by clustering the adjacent points of negative terrain points and negative terrain nest points in TIN_DDM negative .

9. The method of claim 1, 2, 5, 7, or 8, wherein, In step j, the forward slope fold line Γ under the spatial scale threshold R condition Pos-Slope-R The determination condition is: In the formula, Γ Pos-Slope-R , Γ Neg-Slope-R respectively represent the positive terrain slope break line, the negative terrain slope break line under the condition of the spatial scale threshold R; formula (7) represents: under the condition of the spatial scale threshold R, the terrain feature line Γ Feature-R is composed of four lines of the ridge line Γ Ridge-R , the valley line Γ Valley-R , the positive terrain slope break line Γ Pos-Slope-R and the negative terrain slope break line Γ Neg-Slope-R .

Images