A highway road domain deformation field construction method based on point cloud data
Patent Information
- Application Number
- CN202610527424.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-04-21
- Publication Date
- 2026-08-04
Smart Images
Figure CN122510491A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of point cloud data processing and highway geological disaster monitoring technology, specifically a method for constructing a highway roadway deformation field based on point cloud data. Background Technology
[0002] The stability of infrastructure such as slopes and roadbeds within the highway area is directly related to driving safety. Geological disasters such as landslides, collapses, and flood damage are characterized by their sudden occurrence and high destructiveness, often manifesting in their early stages as minute deformations on the order of millimeters to centimeters. Traditional manual inspections and discrete point measurements are insufficient for large-scale, high-frequency, and accurate monitoring, and cannot effectively capture continuous deformation signals during the disaster's incubation process. With the development of UAV-borne lidar technology, high-precision three-dimensional point cloud data of the ground surface can be rapidly acquired, providing a new technical means for monitoring geological disasters in highway areas.
[0003] However, existing deformation analysis methods based on point cloud data mostly rely on interpolating single-period point clouds to generate digital elevation models for comparison. This method is susceptible to uneven point cloud density and interpolation algorithm errors, especially in areas with complex terrain such as highway slopes and high embankment subgrades, making it difficult to guarantee the detection accuracy and reliability of millimeter-level minute deformations from the data source. Furthermore, existing change detection methods are mostly designed for building walls or specific structures, lacking optimization for the complex terrain of highway areas, resulting in insufficient early identification capabilities for hazards such as slope instability.
[0004] Therefore, how to overcome the accuracy loss caused by the existing methods relying on digital elevation model interpolation and comparison, realize the automated and high-precision extraction of millimeter-level minute deformations in the road area, and optimize the deformation detection algorithm for the characteristics of geological disaster monitoring are the technical problems to be solved in this field. Summary of the Invention
[0005] The purpose of this invention is to provide a method for constructing a roadway deformation field based on point cloud data, so as to solve the problems raised in the prior art.
[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for constructing a highway roadway deformation field based on point cloud data. The methods include: Step S1: Collect several periods of raw point cloud data from LiDAR in the target area of the highway according to the time period, and perform noise reduction, point cloud filtering, ground point extraction and georeference processing on the raw point cloud data of each period in sequence to obtain the digital ground model point cloud corresponding to the raw point cloud data of each period. Step S2: Using the digital ground model point cloud acquired in the first cycle as the reference point cloud, the digital ground model point cloud in each subsequent cycle is registered with the reference point cloud with high precision, so that the point cloud data of all cycles are unified under the same spatial coordinate system, and the registered point cloud data of several cycles is obtained. Step S3: Using the point cloud registered in the first cycle as the reference point cloud, and the point clouds registered in each subsequent cycle as the target point cloud, the three-dimensional point cloud change detection algorithm is used to calculate the preliminary three-dimensional deformation field of each subsequent cycle relative to the reference cycle. Step S4: Perform confidence assessment and multi-scale fusion on the preliminary three-dimensional deformation field to obtain the final three-dimensional deformation field of each subsequent period relative to the reference period; Step S5: Project the final three-dimensional deformation field between each cycle to generate a two-dimensional deformation field map. The two-dimensional deformation field map includes the deformation, significance level, and confidence information of each pixel position.
[0007] Furthermore, step S1 includes: Step S1-1: Collect raw point cloud data of the LiDAR in the target area of the highway according to a preset cycle, and obtain the raw point cloud data for each cycle, forming a raw point cloud data set D={D1,D2,…,D…} t ,...,D T}, where D1, D2, ..., D t ,...,D T Let P1, P2, ..., Pt represent the original point cloud data for the 1st, 2nd, ..., tth, ..., Tth periods, respectively; the spatial point sets corresponding to the original point cloud data for each period constitute the point set set P = {P1, P2, ..., PtT}. t ,...,P T}, where P1, P2, ..., P t ,...,P T Let P represent the spatial point sets of the 1st, 2nd, ..., tth, ..., Tth periods, respectively; and the spatial point set P of the tth period. t ={p t,1 ,p t,2 ,…,p t,i ,...,p t,m}, where p t,1 ,p t,2 ,…,p t,i ,...,p t,I Let P represent the spatial point set of the t-th period respectively. t The first spatial point, the second spatial point, ..., the i-th spatial point, ..., the I-th spatial point, each spatial point contains three-dimensional spatial coordinate information; Step S1-2: Use the statistical outlier removal algorithm to process the spatial point set P in the t-th period.t Denoising is performed on the spatial point set P of the t-th period. t The i-th spatial point p t,i Calculate the spatial point set P of the t-th period. t The i-th spatial point p t,i With the spatial point set P in the t-th period t The average spatial distance between the k nearest spatial points in the middle Euclidean distance is denoted as d. t,i The average distance {d} between all m spatial points in the spatial point set of the t-th period. t,1 ,d t,2 ,…,d t,i ,...,d t,I} Calculate the global average distance μ of the spatial point set in the t-th period. t and global standard deviation σ t ,in , ; will satisfy condition |d t,i -μ t |>n·σ t Spatial points are identified as noise points, and the spatial point set P in the t-th period is selected. t Remove from the middle to obtain the denoised spatial point set Q in the t-th cycle. t ={q t,1 ,q t,2 ,…,q t,j ,…,q t,J}, where q t,1 ,q t,2 ,…,q t,j ,…,q t,J Let Q represent the set of spatial points after denoising in the t-th period. t The first spatial point, the second spatial point, ..., the j-th spatial point, ..., the J-th spatial point, and J≤I; where n is the standard deviation multiple parameter; Step S1-3: Use the progressive triangulation encryption algorithm to denoise the spatial point set Q after the t-th period. t Filtering and ground point extraction are performed to construct an initial triangulation, setting a distance threshold L and an angle threshold θ. In each iteration, based on the currently constructed triangulation, the denoised spatial point set Q in the t-th cycle is processed. t The j-th spatial point q t,j Determine the triangular face into which the vertical projection of the point falls, and calculate the denoised spatial point set Q in the t-th period. t The j-th spatial point q t,j The perpendicular distance l to the triangular face t,j And the spatial point set Q after denoising in the t-th cycle t The j-th spatial point qt,j The included angle α between the normal vector and the triangular face normal vector t,j ; Determine the spatial points that satisfy the condition l t,j < L and α t,j < θ as ground points, and add them to the triangular network to encrypt and update the triangular network structure; Iteratively perform the above determination and encryption process until no new ground points are added, and obtain the ground point cloud point set G t ={g t,1 , g t,2 ,…, g t,k ,…, g t,K} separated from the point cloud in the t-th cycle, where g t,1 , g t,2 ,…, g t,k ,…, g t,K respectively represent the first spatial point, the second spatial point, …, the k-th spatial point, …, the K-th spatial point in the ground point cloud point set G t in the t-th cycle, and K ≤ J; where, L is a distance threshold parameter used to control the maximum allowable distance from the spatial point to the triangular face; θ is an angle threshold parameter used to control the maximum allowable included angle between the spatial point normal vector and the triangular face normal vector; Step S1-4: Perform georeferencing processing on the ground point cloud point set G t in the t-th cycle using the positioning data and attitude data recorded by the UAV system, and read the POS position data and IMU attitude data synchronously recorded with the original point cloud data D t in the t-th cycle; Through integrated navigation solution, convert each spatial point in the ground point cloud point set G t from the vehicle coordinate system to the absolute geographic coordinate system, and obtain the digital terrain model point cloud point set Dtm t ={dtm t,1 , dtm t,2 ,…, dtm t,k ,…, dtm t,K} with accurate geographic coordinates in the t-th cycle, where dtm t,1 , dtm t,2 ,…, dtm t,k ,…, dtm t,K respectively represent the first spatial point, the second spatial point, …, the k-th spatial point, …, the K-th spatial point in the digital terrain model point cloud point set Dtm t in the t-th cycle; Step S1-5: Repeat steps S1-2 to step S1-4 for the original point cloud data of each cycle, obtain the digital terrain model point cloud point sets corresponding to the original point cloud data of each cycle, and form a digital terrain model point cloud data set Dtm={Dtm1, Dtm2, …, Dtm t,…,Dtm T}, where Dtm1, Dtm2, ..., Dtm t ,…,Dtm T Let represent the point cloud set of the digital ground model for the 1st cycle, the 2nd cycle, ..., the tth cycle, ..., the Tth cycle, respectively.
[0008] Furthermore, step S2 includes: Step S2-1: Using the point cloud set Dtm1 of the digital ground model in the first cycle as the reference point cloud, and using the point cloud set Dtm1 of the digital ground model in the t-th cycle... t As the target point cloud, a stable control point set F is intelligently selected from both the reference point cloud and the target point cloud. This is achieved by extracting local features from stable and invariant feature regions in the two periodic point clouds and robustly selecting the optimal matching point pairs using a random sample consensus algorithm. t ={f t,1 ,f t,2 ,…,f t,a ,…,f t,A}, where f t,1 ,f t,2 ,…,f t,a ,…,f t,A These represent the 1st control point, the 2nd control point, ..., the ath control point, ..., the Ath control point in the stable control point pair matched with the reference point cloud during the t-th period; the stable control point set is selected from stable ground features outside the monitoring area to establish a unified deformation analysis benchmark. Step S2-2: Based on the stable control point set F t The weighted iterative nearest point algorithm is used to analyze the target point cloud Dtm. t Perform fine registration; in each iteration, for the target point cloud Dtm t For each spatial point in the target point cloud Dtm, find the nearest matching point in the reference point cloud Dtm1 using Euclidean distance to form a set of matching point pairs; t The k-th spatial point dtm t,k Where k = 1, 2, ..., K, K is the target point cloud Dtm t The total number of spatial points; in the reference point cloud Dtm1, find the spatial point with the closest Euclidean distance as the matching point, denoted as u. t,k Forming matching point pairs (dtm) t,k ,u t,k In the error function, each matching point pair (dtm) is considered. t,k ,u t,k Set a weight w t,k The weight w t,kThe weight is determined based on the Euclidean distance between matching point pairs; the larger the matching distance, the smaller the weight. The optimal transformation matrix is obtained iteratively by minimizing the weighted error function, transforming the target point cloud Dtm. t The k-th spatial point dtm t,k Transform the coordinates of the reference point cloud Dtm1 to obtain the transformed spatial point, denoted as r. t,k All transformed spatial points constitute the point cloud set Regt={r} after registration in the t-th period. t,1 ,r t,2 ,…,r t,k ,…,r t,K}, where r t,1 ,r t,2 ,…,r t,k ,…,r t,K Let Dtm represent the original target point cloud in the t-th period. t The first spatial point dtm in t,1 The second spatial point dtm t,2 ..., the kth spatial point dtm t,k ..., the Kth spatial point dtm t,K The new spatial points obtained after registration transformation; Step S2-3: Perform registration on the point cloud set Reg t Registration accuracy assessment is performed based on the stable control point set F. t Calculate the root mean square error of the registration residuals. , where e t,a The residual of the a-th stable control point pair in the t-th period after registration is the target point cloud control point f in the stable control point pair. t,a The Euclidean distance between the spatial points obtained after registration transformation and the corresponding control points in the reference point cloud Dtm1 is required to achieve RMSE. t Less than the preset threshold ε; Step S2-4: Repeat steps S2-1 to S2-3 for the target point cloud in each period, that is, process t=2,3,…,T sequentially to obtain the registered point cloud set for each period, forming a registered point cloud data set Reg={Reg1,Reg2,…,Reg...} for several periods. t ,…,Reg T}, where Reg1 represents the reference point cloud set Dtm1 itself for the first period, Reg2, ..., Reg t ,…,Reg T Let represent the point cloud sets after registration in the 2nd, ..., tth, ..., Tth periods, respectively.
[0009] Furthermore, step S3 includes: Step S3-1: Using the point cloud set Reg1 registered in the first period as the reference point cloud, and using the point cloud set Reg1 registered in the t-th period as the reference point cloud... t As the target point cloud; analysis points of a regular grid are generated on the reference point cloud Reg1 to form the analysis point set H. t ={h t,1 ,h t,2 ,…,h t,v ,…,h t,V}, where h t,1 ,h t,2 ,…,h t,v ,…,h t,V Let represent the 1st, 2nd, ..., vth, ..., Vth analysis points generated for the target point cloud in the t-th period, respectively; Step S3-2: For the analysis point set H t The v-th analysis point h in t,v Search in the baseline point cloud Reg1 to analyze point set H t The v-th analysis point h in t,v For all spatial points within a spherical neighborhood centered at r1, calculate and analyze the point set H based on the spatial points within the neighborhood. t The v-th analysis point h in t,v The normal vector is denoted as m. t,v The displacement direction along the monitoring area determined by prior geological information includes one of the main sliding direction of a landslide, the collapse direction of a collapse, the descent direction of a ground collapse, and the expansion direction of a ground fissure; the prior geological information includes at least one of geological survey data, historical deformation records, and rock strata occurrence data, used to identify the type of geological hazard and its corresponding deformation direction; the normal vector is adjusted; the normal vector m is... t,v Rotate the monitored area by a preset angle δ in the direction of displacement to obtain the adjusted normal vector, denoted as n. t,v ; Step S3-3: Along the adjusted normal vector n t,v Direction, to analyze point set H t The v-th analysis point h in t,v Construct a cylinder with radius r2 and height h centered at the reference point cloud Reg1. Extract the set of spatial points within the cylinder, denoted as B. t,v ={b t,v,1 ,b t,v,2 ,…,b t,v,x}, where b t,v,1 ,b t,v,2 ,…,b t,v,x These represent the first, second, ..., xth spatial points in the reference point cloud Reg1 that fall into the cylinder; and the target point cloud Reg...t Extract the set of spatial points within the same cylinder, denoted as C. t,v ={c t,v,1 ,c t,v,2 ,…,c t,v,y}, where c t,v,1 ,c t,v,2 ,…,c t,v,y Representing the target point cloud Reg t The first spatial point, the second spatial point, ..., the yth spatial point fall into the cylinder; Step S3-4: Calculate the set of spatial points B in the reference point cloud Reg1 that fall within the cylinder. t,v All spatial points along the adjusted normal vector n t,v The average value of the projected coordinates of the direction, denoted as b. t,v,mean ; Calculate the target point cloud Reg t The set of spatial points C falling into the cylinder t,v All spatial points along the adjusted normal vector n t,v The average value of the projected coordinates of the direction, denoted as c. t,v,mean Then analyze the point set H t The v-th analysis point h in t,v The shape variable d at the center t,v =c t,v,mean -b t,v,mean ; Step S3-5: Analyze the point set H t Perform steps S3-2 to S3-4 at all analysis points to obtain the preliminary three-dimensional deformation field of the t-th period relative to the 1st period, denoted as Δ. t,pre ={(h t,v ,d t,v )|v=1,2,…,V}, where each element contains the analysis point position and deformation.
[0010] Furthermore, step S4 includes: Step S4-1: For the initial three-dimensional deformation field Δ of the t-th period relative to the 1st period t,pre The v-th analysis point h in t,v Calculate the deformation d t,v The corresponding standard deviation σ t,v Take the set of spatial points B in the reference point cloud Reg1 that fall within the cylinder. t,v and target point cloud Reg t The set of spatial points C falling into the cylinder t,v The average point cloud density is used as the analysis point set H. t The v-th analysis point h in t,v The local point cloud density ρ at that location t,vThe analysis point set H is obtained through interpolation. t The v-th analysis point h in t,v Registration residual f at the location t,v ; will σ t,v ρ t,v f t,v Normalization is performed separately to obtain the normalized value σ. t,v,norm ρ t,v,norm f t,v,norm Construct a comprehensive confidence index z t,v =ω1·(1-σ t,v,norm )+ω2·ρ t,v,norm +ω3·(1-f t,v,norm ), where ω1, ω2, and ω3 are preset weight coefficients and satisfy ω1 + ω2 + ω3 = 1; let z t,v The deformation d corresponding to analysis points below the preset confidence threshold γ t,v Treat it as noise and filter it out; Step S4-2: Using multiple sets of preset parameter combinations (r1, r2, h), execute steps S3-2 to S3-5 respectively to obtain the deformation at different scales, denoted as d. t,v,1 ,d t,v,2 ,…,d t,v,L , where d t,v,1 ,d t,v,2 ,…,d t,v,L These represent the analysis point set H obtained by calculating using the parameter combinations of group 1, group 2, ..., group L, respectively. t The v-th analysis point h in t,v The deformation at each point is L, where L is the total number of parameter combinations; the deformation at different scales is weighted and fused to obtain the analysis point set H. t The v-th analysis point h in t,v The shape variable after fusion o t,v ; Step S4-3: Analyze the point set H t The v-th analysis point h in t,v Repeat steps S3-2 to S3-6 to obtain the final three-dimensional deformation field of the t-th period relative to the 1st period, denoted as Δt={(h t,v ,o t,v ,σ t,v ,z t,v )|v=1,2,…,V}, where each element contains the analysis point location, the fused deformation, the standard deviation, and the confidence level information; Step S4-4: Repeat steps S3-1 to S3-7 for each period's target point cloud, that is, process t=2,3,…,T sequentially to obtain the final three-dimensional deformation field of each period relative to the reference period, forming the final deformation field set Δ={Δ2,Δ3,…,Δ t ,…,Δ T}, where Δ2, Δ3, ..., Δ t ,…,Δ T Let represent the final three-dimensional deformation fields of the 2nd period, 3rd period, ..., t-th period, and T-th period relative to the 1st period, respectively.
[0011] Furthermore, step S5 includes: Step S5-1: For the final three-dimensional deformation field Δ of the t-th period relative to the 1st period t Extract the three-dimensional deformation field Δ of the t-th period relative to the 1st period. t The vth analysis point h t,v The plane coordinates are denoted as (x t,v ,y t,v ), and with the v-th analysis point h t,v The fused shape variable o t,v Standard deviation σ t,v and confidence level z t,v As the vth analysis point h t,v The deformation property value; Step S5-2: Arrange all analysis points according to planar coordinates to generate a two-dimensional deformation field spectrum for the t-th period. The two-dimensional deformation field spectrum includes a deformation variable layer, a significance level layer, and a confidence level layer. The deformation variable layer is composed of the fused deformation variable o from all analysis points. t,v The significance level layer is composed of the standard deviation σ of each analysis point. t,v The confidence layer is composed of the confidence scores z for each analysis point. t,v constitute; Step S5-3: Repeat steps S4-1 to S4-2 for t=2,3,…,T to obtain the two-dimensional deformation field spectrum corresponding to each period, forming a set of two-dimensional deformation field spectra.
[0012] Compared with the prior art, the beneficial effects of the present invention are: 1. This invention directly processes multiple phases of original point cloud data, avoiding the accuracy loss caused by traditional methods that rely on digital elevation model interpolation comparisons. By combining statistical outlier removal algorithms, progressive triangulation encryption algorithms, and weighted iterative nearest-point algorithms, it achieves high-precision denoising, ground point extraction, and registration of original point clouds, providing a clean and reliable data foundation for subsequent deformation analysis, thereby enabling accurate extraction of millimeter-level minute deformations in highway roadways; 2. This invention optimizes and improves the M3C2 algorithm to address the characteristics of geological disaster monitoring. By introducing normal vector fine-tuning based on prior geological information, the deformation calculation is made more consistent with the actual displacement direction of disasters such as landslides. By constructing a comprehensive confidence index that includes standard deviation, local point cloud density, and registration residuals, the interference of noise points on the deformation results is effectively filtered out, significantly improving the reliability and accuracy of deformation detection.
[0013] 3. This invention employs a multi-scale analysis strategy, using different parameter combinations to calculate multiple sets of deformation variables and then weighted and fused them. This allows for the simultaneous capture of deformation characteristics at different spatial scales, ranging from micro-cracks to overall slippage, thus meeting the needs of different development stages of geological hazards. Furthermore, by generating a two-dimensional deformation field map containing deformation variables, significance levels, and confidence information, it provides intuitive and quantitative data support for the early identification and risk classification of geological hazards in highway areas. Attached Figure Description
[0014] Figure 1 This is a flowchart illustrating a method for constructing a roadway deformation field based on point cloud data according to the present invention. Detailed Implementation
[0015] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0016] Example: Figure 1 As shown, this invention provides a technical solution: a method for constructing a roadway deformation field based on point cloud data. The right-side cutting slope of a highway section from K35+200 to K35+400 was used as the monitoring target. Data was collected over three cycles using an unmanned aerial vehicle (UAV)-borne lidar system: Cycle 1 was March 1, XXXX (baseline period); Cycle 2 was June 1, XXXX; and Cycle 3 was September 1, XXXX. Only 11 spatial points on the slope surface were collected as example data in each cycle, and all coordinates were based on a local coordinate system (unit: meters). Two stable control points, (0,0) and (2,2), were set around the perimeter of the slope monitoring area, and their coordinates remained unchanged throughout the cycles.
[0017] The methods include: Step S1: Collect several periods of raw point cloud data from LiDAR in the target area of the highway according to the time period, and perform noise reduction, point cloud filtering, ground point extraction and georeference processing on the raw point cloud data of each period in sequence to obtain the digital ground model point cloud corresponding to the raw point cloud data of each period. Step S1-1: Collect raw point cloud data for three cycles, with each cycle containing 10 spatial points; The first periodic point set P1 is: p 1,1 =(0,0,10), p 1,2 =(1,0,10.05), p 1,3 =(2,0,10.1), p 1,4 =(0,1,10.02), p 1,5 =(1,1,10.07), p 1,6 =(2,1,10.12), p 1,7 =(0,2,10.03), p 1,8 =(1,2,10.08), p 1,9 =(2.00,2.00,10.13), p 1,10 =(1.5,1.5,10.5), p 1,11 =(3.00,3.00,10.5); The second periodic point set P2 is: p 2,1 =(0,0,10), p 2,2 =(1,0,10.07), p 2,3 =(2,0,10.13), p 2,4 =(0,1,10.03), p 2,5 =(1,1,10.09), p 2,6 =(2,1,10.15), p 2,7 =(0,2,10.04), p 2,8 =(1,2,10.1), p 2,9 =(2,2,10.13), p 2,10 =(1.5,1.5,10.52), p 2,11 =(3,3,11); The third periodic point set P3 is: p 3,1 =(0,0,10), p 3,2 =(1,0,10.09), p 3,3 =(2,0,10.16), p 3,4 =(0,1,10.04), p 3,5 =(1,1,10.11), p 3,6 =(2,1,10.18), p 3,7=(0,2,10.05), p 3,8 =(1,2,10.12), p 3,9 =(2,2,10.13), p 3,10 =(1.5,1.5,10.55), p 3,11 =(3,3,11.50); Step S1-2: Use the statistical outlier removal algorithm for noise reduction, where k=3 is the number of nearest neighbors and n=2 is the standard deviation factor; in the first period, calculate the average Euclidean distance d between each spatial point and its three nearest neighbors. 1,1 =1.138, d 1,2 =1.001, d 1,3 =1.334, d 1,4 =1,d 1,5 =1.001, d 1,6 =1,d 1,7 =1.138, d 1,8 =1.001, d 1,9 =1.334, d 1,10 =0.818, d 1,11 =2.002; Global mean μ1=1.161, standard deviation σ1=0.304, threshold T=n·σ1=0.608; Check |d 1,i -μ1|:p 1,11 The difference between the points is 1 > 0.788, so they are identified as noise points and removed. The remaining 10 spatial points are retained, resulting in the denoised point set Q1 = {p 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,p 1,5 ,p 1,6 ,p 1,7 ,p 1,8 ,p 1,9 ,p 1,10 Similarly, p is calculated for the second period. 2,11 For noise points, Q2={p 2,1 ,p 2,2 ,p 2,3 ,p 2,4 ,p 2,5 ,p 2,6 ,p 2,7 ,p 2,8 ,p 2,9 ,p 2,10}; For the third period, we get p 3,11 For noise points, Q3={p 3,1 ,p 3,2 ,p 3,3 ,p 3,4 ,p 3,5 ,p3,6 , p 3,7 , p 3,8 , p 3,9 , p 3,10}。
[0018] Step S1-3: Filter and extract ground points using the progressive triangulation encryption algorithm, set the distance threshold L = 0.2 m and the angle threshold θ = 10°. For the first cycle Q1, first select the lowest point to construct the initial triangulation, and select p 1,1 (0, 0, 10), p 1,2 (1, 0, 10.05), p 1,4 (0, 1, 10.02) to form triangle T0; perform iterative determination on each unprocessed point: point p 1,3 (2, 0, 10.1): The projection falls into T0, calculate the vertical distance l 1,3 ≈0.077 m < L, the included angle α of the normal vector 1,3 ≈5° < θ, meets the conditions, is determined as a ground point, and added to the triangulation; point p1,5(1, 1, 10.07): The projection falls into a triangular face of the updated triangulation, calculate l 1,5 ≈0.05 m < L, α 1,5 ≈4° < θ, is determined as a ground point, and added to the triangulation; similarly, p 1,6 , p 1,7 , p 1,8 , p 1,9 , etc. all meet the conditions and are added in sequence; point p 1,10 (1.5, 1.5, 10.5): The projection falls into a triangular face, calculate the elevation at the projection point of this triangular face ≈10.09 m, the vertical distance l 1,10 =10.50 - 10.09 = 0.41 m > L, and the included angle α of the normal vector 1,10 ≈15° > θ, does not meet the conditions, and is excluded. Iterate until no new points are added, and obtain the ground point cloud point set G1 = {p 1,1 , p 1,2 , p 1,3 , p 1,4 , p 1,5 , p 1,6 , p 1,7 , p 1,8 , p 1,9}}, a total of 9 spatial points; similarly, for the second cycle Q2, exclude the vegetation point p 2,10 , and get G2 = {p 2,1 , p 2,2 , p 2,3 , p 2,4 , p 2,5 , p 2,6 , p 2,7 , p2,8 ,p 2,9 For period 3 Q3, remove p. 3，10 , thus G3={p 3,1 ,p 3,2 ,p 3,3 ,p 3,4 ,p 3,5 ,p 3,6 ,p 3,7 ,p 3,8 ,p 3,9 ,p 3,10}; Step S1-4: Use UAV POS / IMU data to convert the point cloud to the absolute geographic coordinate system. The coordinates remain unchanged after conversion, and the point cloud set of the digital ground model is obtained: Dtm1=G1, Dtm2=G2, Dtm3=G3. Step S1-5: Construct a digital ground model point cloud dataset Dtm={Dtm1,Dtm2,Dtm3}.
[0019] Step S2: Using the digital ground model point cloud acquired in the first cycle as the reference point cloud, the digital ground model point cloud in each subsequent cycle is registered with the reference point cloud with high precision, so that the point cloud data of all cycles are unified under the same spatial coordinate system, and the registered point cloud data of several cycles is obtained. Step S2-1: Using Dtm1 as the reference point cloud and Dtm2 and Dtm3 as the target point clouds; by extracting local features of stable and invariant feature regions from the point clouds of the two periods, and using the random sampling consensus algorithm to robustly select the optimal matching point pairs, a stable control point set is obtained; stable ground features on the periphery of the monitoring area are selected as control points, with reference points (0,0,10) and (2,2,10,13); after feature matching, the target point positions corresponding to the second period are p 2,1 (0,0,10) and p 2,9 (2.001, 2.001, 10.132); the target point positions corresponding to the 3rd cycle are p 3,1 (0,0,10.00) and p 3,9 (2.002, 1.998, 10.135); The stable control point set is obtained: F2 = {f 2,1 :(0,0,10.00) (0,0,10.00),f 2,2 :(2,2,10.13) (2.001, 2.001, 10.132)};F3={f 3,1 :(0,0,10.00) (0,0,10.00),f 3,2 :(2,2,10.13) (2.002, 1.998, 10.135)}; Step S2-2: Based on the stable control point set, a weighted iterative nearest-point algorithm is used to perform fine registration of the target point cloud. In each iteration, for each spatial point in the target point cloud Dtm2, the nearest spatial point in the reference point cloud Dtm1 is found, and weights are assigned according to the distance. After iterative solution, the optimal transformation matrix (containing small rotations and translations) is obtained, transforming the target point cloud to the reference coordinate system. After the transformation, the coordinates of the points in the original target point cloud undergo slight changes, and the control point p... 2,1 (0,0,10.00) transforms into (0.001,-0.002,10.003), p 2,9 (2.001, 2.001, 10.132) is transformed into (2.002, 1.998, 10.131); all the transformed points constitute the registered point cloud Reg2; similarly, registering Dtm3 yields Reg3, and the transformed p 3,1 It becomes (0.002, -0.001, 10.004), p 3,9 It becomes (2.003, 1.997, 10.136); Step S2-3: Based on the stable control point set, calculate the residuals of the registered control point pairs: For t=2, e 2,1 ≈0.0037m; e 2,2 ≈0.003m; RMSE2≈0.0034m; For t=3, e 3,1 ≈0.0046m; e 3,2 ≈0.0061m; all are less than the preset threshold ε=0.05m, meeting the requirements; Step S2-4: Form the registered point cloud set Reg={Reg1=Dtm1,Reg2,Reg3}.
[0020] Step S3: Using the point cloud registered in the first cycle as the reference point cloud, and the point clouds registered in each subsequent cycle as the target point cloud, the three-dimensional point cloud change detection algorithm is used to calculate the preliminary three-dimensional deformation field of each subsequent cycle relative to the reference cycle. Step S3-1: Generate regular grid analysis points on the reference point cloud Reg1, with a spacing of 1m, for a total of 9 spatial points. For the 2nd and 3rd cycles, construct analysis point sets H2 and H3 respectively. The spatial positions of all analysis points are the same. The interpolation value from the reference point cloud is: h 2,1 =(0,0,10),h 2,2 =(1,0,10.05),h 2,3 =(2,0,10.1),h 2,4 =(0,1,10.02),h2,5 =(1,1,10.07), h 2,6 =(2,1,10.12),h 2,7 =(0,2,10.03),h 2,8 =(1,2,10.08),h 2,9 =(2,2,10.13); h 3,1 =(0,0,10),h 3,2 =(1,0,10.05),h 3,3 =(2,0,10.1),h 3,4 =(0,1,10.02),h 3,5 =(1,1,10.07), h 3,6 =(2,1,10.12),h 3,7 =(0,2,10.03),h 3,8 =(1,2,10.08),h 3,9 =(2,2,10.13).
[0021] Step S3-2: Using h 2,1 Let r1 = 1.5m, and search for neighborhood points in the baseline point cloud Reg1: (0,0,10), (1,0,10.05), (0,1,10.02), (1,1,10.07); fit the plane to obtain the normal vector m. 2,1 =(0.05,0.02,-1), after normalization it becomes (0.0499,0.0200,-0.9985); based on the prior sliding direction southeast and the rotation angle δ=0°, the adjusted normal vector n is obtained. 2,1 =m 2,1 Similar calculations were performed on other analysis points to obtain their respective normal vectors.
[0022] Step S3-3: Use the first set of parameters (r1=1.5, r2=0.5, h=0.2); with h 2,1 Along n 2,1 Construct a cylinder (radius 0.5m, height 0.2m) and extract B from the baseline point cloud Reg1. 2,1 ={(0,0,10)}, extract C from the target point cloud Reg2. 2,1 ={(0.001,-0.002,10.003)}; Similarly, for t=3, h 3,1 Corresponding B 3,1 ={(0,0,10)},C 3,1 ={(0.002,-0.001,10.004)}; Step S3-4: Calculate the average value of the projected coordinates to obtain the deformation d.t,v =c t,v,mean -b t,v,mean For h 2,1 d 2,1 =10.003-10=0.003m; similarly, calculate all points to obtain the preliminary deformation d of the second period. 2,v :h 2,1 =0.003,h 2,2 =0.02,h 2,3 =0.03,h 2,4 =0.01,h 2,5 =0.02,h 2,6 =0.03,h 2,7 =0.01,h 2,8 =0.02, h 2,9 =0.03; Preliminary deformation of the 3rd period: d 3,1 =0.004,d 3,2 =0.04,d 3,3 =0.06,d 3,4 =0.02,d 3,5 =0.04,d 3,6 =0.06,d 3,7 =0.02,d 3,8 =0.04,d 3,9 =0.06; Step S3-5: Constructing the preliminary three-dimensional deformation field: Δ 2,pre ={(h 2,v ,d 2,v )},Δ 3,pre ={(h 3,v ,d 3,v )}.
[0023] Step S4: Perform confidence assessment and multi-scale fusion on the preliminary three-dimensional deformation field to obtain the final three-dimensional deformation field of each subsequent period relative to the reference period; Step S4-1: Calculate the confidence index and filter out low confidence points; using h 2,1 Standard deviation σ 2,1 =0 (single point on a cylinder), local density ρ 2,1 =1 / (π×0.5 2 ×0.2)≈6.37, registration residual f 2,1 The control point f is obtained by interpolation using its planar position. 2,1 and f 2,2 The plane coordinates are (0,0) and (2,2), with corresponding registration residuals of 0.0037m and 0.003m; using inverse distance weighted interpolation, h 2,1 Located at (0,0), therefore f 3,1≈0.0037m; Let σ max =0.1, ρ max =10, f max =0.05, then the normalized value is: σ 2,1,norm =σ 2,1 / σ max =0, ρ 2,1,norm =ρ 2,1 / ρ max =0.637, f 2,1,norm =f 2,1 / f max =0.074; Taking weights ω1=0.3, ω2=0.4, ω3=0.3, the overall confidence level z 2,1 =ω1·(1-σ 2,1,norm )+ω2·ρ 2,1,norm +ω3·(1-f 2,1,norm )=0.8326; Preset reliability threshold γ=0.5, z 2,1 If the confidence level is greater than 0.5, the analysis point is retained. Similar calculations are performed on all analysis points, and the confidence level of each analysis point is greater than 0.5, so no points are filtered out. Step S4-2: Use L=2 preset parameter combinations: Group 1 is (r1=1.5, r2=0.5, h=0.2), and Group 2 is (r1=2, r2=1, h=0.3); the deformation d of Group 1 2,v,1 This refers to the initial deformation d calculated in step S3-4. 2,v The shape variable d of group 2 2,v,2 Recalculation is required; using h 2,1 By re-executing steps S3-2 to S3-4 using parameters from group 2, we obtain d. 2,1,2 =0.0125m; Similarly, the deformation d of all analysis points in group 2 is obtained. 2,v,2 Using equal-weight fusion, the fused deformation variable o is obtained. 2,v =(d 2,v,1 +d 2,v,2 ) / 2; The specific results are as follows: o 2,1 =0.00775, o 2,2 =0.0225, o 2,3 =0.03375, o 2,4 =0.01125, o 2,5 =0.0225, o 2,6 =0.03375, o 2,7 =0.01125, o 2,8 =0.0225, o 2,9 =0.03375; For the 3rd period, the deformation d of group 2 3,v,2 For: d 3,1,2=0.025, d 3,2,2 =0.05, d 3,3,2 =0.075, d 3,4,2 =0.025, d 3,5,2 =0.05, d 3,6,2 =0.075, d 3,7,2 =0.025, d 3,8,2 =0.05, d 3,9,2 =0.075; fusion yields: o 3,1 =0.0145, o 3,2 =0.045, o 3,3 =0.0675, o 3,4 =0.0225, o 3,5 =0.045, o 3,6 =0.0675, o 3,7 =0.0225, o 3,8 =0.045, o 3,9 =0.0675; Step S4-3: Integrate the fused deformation, standard deviation, and confidence level to obtain the final three-dimensional deformation Δ2={(h 2,v ,o 2,v ,σ 2,v ,z 2,v The final three-dimensional deformation of the third period is Δ3 = {(h}. 3,v ,o 3,v ,σ 3,v ,z 3,v ); Step S4-4: Form the final deformation field set Δ={Δ2,Δ3}.
[0024] Step S5: Project the final three-dimensional deformation field between each cycle to generate a two-dimensional deformation field map, which includes the deformation, significance level and confidence information of each pixel position; Step S5-1: For each analysis point in the final three-dimensional deformation field, extract its planar coordinates (x, y, z). t,v ,y t,v And deformation attribute values: the fused deformation value o t,v Standard deviation σ t,v Confidence level z t,v ; for h in Δ2 2,1 The plane coordinates are (0,0), o 2,1 =0.00775, σ 2,1 =0, z 2,1 =0.8326; Step S5-2: Arrange all analysis points according to planar coordinates to generate the deformation layer, significance level layer, and confidence layer for the second and third periods; Step S5-3: Obtain the two-dimensional deformation field spectra corresponding to the second and third periods, forming a set of two-dimensional deformation field spectra.
[0025] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.
Claims
1. A method for constructing a roadway deformation field based on point cloud data, characterized in that: Includes the following steps: Raw point cloud data of lidar in the target area of the highway is collected periodically at preset time intervals. The raw point cloud data of each period is preprocessed to obtain the digital ground model point cloud corresponding to each period. Taking the digital ground model point cloud of the first period as the reference point cloud, the digital ground model point cloud of each subsequent period is registered with the reference point cloud with high precision to unify the point cloud data of all periods under the same spatial coordinate system. Taking the registered point cloud of the first period as the reference point cloud, and taking the registered point cloud of each subsequent period as the target point cloud, the preliminary three-dimensional deformation field of each subsequent period relative to the reference period is calculated sequentially using a three-dimensional point cloud change detection algorithm. The initial three-dimensional deformation field is evaluated for confidence and fused at multiple scales to obtain the final three-dimensional deformation field of each subsequent period relative to the reference period. The final three-dimensional deformation field of each period is projected to generate a two-dimensional deformation field map containing the deformation, significance level and confidence information of each pixel position.
2. The method for constructing a highway roadway deformation field based on point cloud data according to claim 1, characterized in that: The preprocessing includes: A statistical outlier removal algorithm is used to denoise the original point cloud data for each period; a progressive triangulation encryption algorithm is used to filter the denoised point cloud and extract ground points to obtain the ground point cloud; georeferencing is performed on the ground point cloud using synchronously recorded positioning and attitude data, and it is converted to an absolute geographic coordinate system to obtain the digital ground model point cloud.
3. The method for constructing a highway roadway deformation field based on point cloud data according to claim 1, characterized in that: The high-precision registration includes: A stable control point set is intelligently selected for both the reference point cloud and the target point cloud. Local features of points in the two periodic point clouds are extracted and matched. A random sampling consensus algorithm is used to select the optimal matching point pairs from the matching results, forming a stable control point set. Based on this stable control point set, a weighted iterative nearest point algorithm is employed to perform fine registration of the target point cloud. In each iteration, the nearest matching point in the reference point cloud is found for each point in the target point cloud, and weights are set according to the matching distance. The optimal transformation matrix is solved by minimizing the weighted error function. The registration accuracy of the registered point cloud is evaluated. The root mean square error of the registration residual is calculated based on the stable control point set, and this root mean square error is less than a preset threshold.
4. The method for constructing a highway roadway deformation field based on point cloud data according to claim 1, characterized in that, The 3D point cloud change detection algorithm includes: deformation calculation, confidence assessment, and multi-scale fusion.
5. The method for constructing a highway roadway deformation field based on point cloud data according to claim 4, characterized in that: The deformation calculation includes: Regular grid analysis points are generated on the reference point cloud. For each analysis point, the normal vector of the analysis point is calculated by searching its neighborhood point set in the reference point cloud. The normal vector is adjusted along the displacement direction of the monitoring area determined by prior geological information. A cylinder with a set radius and a set height is constructed with the analysis point as the center along the adjusted normal vector direction. Point sets falling into the cylinder from the reference point cloud and the target point cloud are extracted respectively. The difference between the average values of the projected coordinates of the reference point cloud point set and the target point cloud point set along the adjusted normal vector direction is calculated as the deformation at the analysis point. The deformation of all analysis points constitutes a preliminary three-dimensional deformation field.
6. The method for constructing a highway roadway deformation field based on point cloud data according to claim 4, characterized in that: The confidence assessment includes: For each analysis point in the preliminary three-dimensional deformation field, the standard deviation of the deformation at each analysis point and the local point cloud density at the analysis point are calculated, and the registration residual at the analysis point is obtained by interpolation. After normalizing the standard deviation, local point cloud density and registration residual of each analysis point, a comprehensive confidence index is constructed, and analysis points and the deformation at the analysis points whose comprehensive confidence index is lower than the confidence threshold are filtered out.
7. The method for constructing a highway roadway deformation field based on point cloud data according to claim 4, characterized in that: The multi-scale fusion includes: For each analysis point retained after filtering, several sets of parameter combinations with different set radii and set heights are used to perform the deformation calculation process, resulting in several deformation variables for each analysis point under different parameter combinations. These deformation variables are then weighted and fused to obtain the fused deformation variable at the analysis point. All retained analysis points are traversed to obtain the fused deformation variable, standard deviation, and confidence level at each analysis point. Using the positions of all retained analysis points as spatial locations and the fused deformation variable, standard deviation, and confidence level at each location as attribute values, the final three-dimensional deformation field of each period relative to the reference period is constructed.
8. The method for constructing a highway roadway deformation field based on point cloud data according to claim 1, characterized in that: The step of projecting the final three-dimensional deformation field between each cycle to generate a two-dimensional deformation field map includes: Extract the planar coordinates of each analysis point in the final three-dimensional deformation field of each period, as well as the corresponding fused deformation, standard deviation, and confidence level, as deformation attribute values; arrange all analysis points according to planar coordinates to generate a two-dimensional deformation field map containing a deformation layer, a significance level layer, and a confidence level layer.