An airborne LIDAR power line extraction method based on target space distribution characteristics

Abstract

The application provides an airborne LIDAR power line extraction method based on target space distribution characteristics, and relates to the technical field of remote sensing data processing. Original airborne LIDAR point cloud data is read, near-ground data is removed from the original airborne LIDAR data, and a non-near-ground data set is obtained; the non-near-ground data set is regularized into a binary 3D virtual voxel data set; based on the spatial geometric characteristics of the power line, linear virtual voxels and laser points in the virtual voxels are extracted, and power line extraction based on virtual voxel spatial geometric characteristics is completed. The method of the application makes good use of the neighborhood relationship between each virtual voxel in the 3D virtual voxel data and the characteristics of the power line, and is helpful to the development of airborne LIDAR data processing and application based on virtual voxel modeling theory.

Description

Technical Field

[0001] This invention relates to the field of remote sensing data processing technology, and in particular to an airborne LIDAR power line extraction method based on the spatial distribution characteristics of targets. Background Technology

[0002] Power lines are an indispensable component of 3D geographic information products, making automatic, high-precision, and rapid power line target extraction a research hotspot. Airborne LiDAR (Light Detection and Ranging) technology can quickly and accurately acquire high-precision, high-density three-dimensional (3D) spatial information of surface features, making rapid extraction of 3D power line information possible. Traditional power line extraction methods can be categorized by data structure into raster-based and pixel-based methods. In raster-based methods, the raw data is organized into a raster structure, and the raster is used as the data processing unit. Power line data within the raster is extracted based on the characteristics of high height and low density, combined with empirical or dynamic threshold comparisons. Raster-based methods are simple to implement and efficient, but the completeness of power line data extraction is low due to processing on a raster-by-raster basis. Pixel-based methods typically require projecting the data onto a horizontal plane and using mature image processing techniques (such as line detection algorithms and least squares methods) for power line extraction. However, this method may lead to data loss during projection, reducing extraction accuracy. Therefore, traditional power line extraction methods are not conducive to leveraging the true 3D technical advantages of airborne LiDAR. Voxel structures are a true 3D data structure that can be used to represent LiDAR data without information loss. Furthermore, the voxels within this structure implicitly contain geometric topological relationships, making the design of data processing algorithms based on this data structure relatively easy. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to address the shortcomings of the prior art by providing an airborne LIDAR power line extraction method based on the spatial distribution characteristics of the target. The method constructs virtual voxels and combines the voxel structure with power line target extraction. The final power line extraction result is expressed as discrete laser points.

[0004] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0005] A method for extracting power lines using airborne LiDAR based on target spatial distribution characteristics includes the following steps:

[0006] Step 1: Read the raw airborne LiDAR data to form the raw airborne LiDAR dataset;

[0007] Step 2: Remove near-ground data from the original airborne LiDAR data to obtain the non-near-ground dataset;

[0008] Step 3: Regularize the non-near-ground data into a virtual 3D voxel structure, and assign a value of 1 to the virtual voxel containing the laser point and a value of 0 to the other virtual voxels to obtain a virtual binary 3D voxel.

[0009] Step 4: Based on the spatial distribution characteristics of the target, use principal component analysis to determine the distribution shape of each 1-valued virtual voxel and the target points in its spatial neighborhood. Mark the virtual voxels that are distributed in a linear pattern, and identify the laser points in the linear virtual voxels as electric field lines, thus completing the extraction of electric field lines based on the spatial distribution characteristics of the target.

[0010] Further, step 2 includes the following steps:

[0011] Step 2.1: Separate ground points and non-ground points in the original airborne LiDAR data based on the cloth simulation filtering algorithm;

[0012] Step 2.2: Interpolate the original airborne LiDAR data and the separated ground data respectively to obtain the digital surface model (DSM) and digital elevation model (DEM), and further subtract the DSM from the DEM to generate a regularized digital surface model (nDSM).

[0013] Step 2.3: Based on the spatial distribution characteristics of elevation values ​​in nDSM, the elevation threshold H0 is dynamically obtained using the "elevation-frequency" histogram analysis method. nDSMs below the elevation threshold H0 are classified as near-ground nDSMs, while those above are classified as non-near-ground nDSMs. Based on the mapping relationship between laser points and grids, the laser points corresponding to non-near-ground nDSMs in each grid are read to obtain the non-near-ground dataset.

[0014] Furthermore, step 2.3 specifically includes the following steps:

[0015] Step 2.3.1: Count the frequency of elevation values ​​in nDSM and visualize the statistical results in the form of a histogram. Record the elevation value mutation point as the elevation threshold H0 of nDSM.

[0016] Step 2.3.2: For the elevation values ​​in nDSM, determine whether the data in nDSM that is below the elevation threshold H0 is near-ground nDSM data, and vice versa.

[0017] Step 2.3.3: Based on the mapping relationship between laser points and grids in nDSM, read the laser points of the grids containing non-near-ground data in each grid to obtain the non-near-ground dataset.

[0018] Furthermore, step 3 specifically includes the following steps:

[0019] Step 3.1: Represent the 3D spatial extent using directed bounding boxes from the non-near-ground dataset;

[0020] Step 3.2: Determine the resolution (Δx, Δy, Δz) of the virtual voxel in the x, y, and z directions based on the average point spacing of the laser points in the non-near-ground data set. Δx, Δy, and Δz are the sizes of the virtual voxels.

[0021] Step 3.3: Based on the virtual voxel resolution (Δx, Δy, Δz), the directed bounding box is virtually divided to obtain a virtual 3D voxel grid. Each virtual 3D voxel grid cell is a virtual 3D voxel.

[0022] Step 3.4: Map each laser point in the non-near-ground data set to a virtual 3D voxel, and then assign the value 1 to the virtual 3D voxel containing the laser point and the value 0 to the other virtual 3D voxels, as shown in Equation (9).

[0023]

[0024] Among them, v s (r s ,c s ,l s ) represents the voxels in the 3D voxel array, and V is the set of voxels in the 3D voxel array, V = {v s (r s ,c s ,l s ), s=1,...,k}, where s is the volume index; k is the number of volume elements; v s It is the voxel value of the s-th voxel; (r s ,c s ,l s ) represents the coordinates of the s-th voxel in the voxel array, and r s c s l s These represent the row number, column number, and layer number, respectively. The floor operator; (x j ,y j ,z j ) represents the coordinates of the j-th point in the corresponding non-near-ground dataset; x min =min{x j ,j=1,…,t},y min =min{y j ,j=1,…,t},z min =min{z j ,j=1,...,t}, where t is the number of laser points in the non-near-ground data set;

[0025] Thus, a virtual binary 3D voxel is obtained, completing the regularization of non-near-ground data.

[0026] Furthermore, step 4 includes the following steps:

[0027] Step 4.1: Initialize the laser point closest to the center of each 1-value virtual voxel as the target point of each 1-value virtual voxel;

[0028] Step 4.2: Based on principal component analysis, construct linear, planar, and spherical distribution indices L to describe the spatial distribution characteristics of the target points of each 1-valued virtual volume element. λ P λ and S λ And the direction vector describing the directional features of the target point of each 1-valued virtual volume element.

[0029] Step 4.3: Determine the linear distribution index L λ Greater than the areal distribution index P λ and spherical distribution index S λ And direction vector The target point of the 1-value virtual volume element that is approximately parallel to the horizontal plane is the seed point of the electric field line. The 1-value virtual volume element containing the seed point of the electric field line is marked as the electric field line virtual volume element. The laser points in the electric field line virtual volume element are read as the electric field line data to complete the extraction of the electric field line.

[0030] Furthermore, step 4.1 specifically includes the following steps:

[0031] Step 4.1.1: Calculate the distance d between all laser points within each 1-value virtual voxel and the center of each 1-value virtual voxel. c ;

[0032] Step 4.1.2: Select the distance between each 1-value virtual voxel and its center as the nearest distance d. c,min The laser point is used as the target point of each 1-value virtual volume element.

[0033] Furthermore, step 4.2 specifically includes the following steps:

[0034] Step 4.2.1: Scan each 1-value virtual volume element, perform a spatial neighborhood search on it, and obtain its neighborhood set;

[0035] Step 4.2.2: Construct the covariance matrix Cov based on each 1-valued virtual voxel target point and its neighborhood target point set, then construct a homogeneous equation, and solve for the eigenvalues ​​λ1, λ2, λ3 of the covariance matrix Cov and the corresponding eigenvectors. Where λ1≥λ2≥λ3, the eigenvector corresponding to the largest eigenvalue λ1 For each 1-valued virtual volume element target point, there is a direction vector.

[0036] Step 4.2.3: Based on the feature values ​​λ1, λ2, and λ3 of each of the above-mentioned one-valued virtual volumetric target points, construct a linear distribution index L describing the distribution characteristics of each one-valued virtual volumetric target point. λ The areal distribution index P λ and spherical distribution index S λ As shown in equation (14):

[0037]

[0038] Among them, L λ S λ P λ ∈[0,1],L λ +S λ +P λ =1.

[0039] The beneficial effects of adopting the above technical solution are as follows: The airborne LiDAR power line extraction method based on target spatial distribution characteristics provided by this invention first removes near-ground data by setting an elevation threshold based on the characteristic that there is a certain distance between power lines and the ground, obtaining a non-near-ground dataset. Then, the non-near-ground data is regularized into virtual voxels, and the power lines within the virtual voxels are extracted according to the linear spatial distribution characteristics of the power lines. Compared with power line extraction methods based on traditional grid, pixel, and other data structures, voxels are a simple spatial structure, making data processing easier. The method proposed in this invention virtually divides the spatial data and transforms target extraction based on discrete laser points into target extraction based on virtual voxels. It effectively utilizes the implicit neighborhood relationships between virtual voxels in 3D voxels and the characteristics of power lines, which contributes to the development of airborne LiDAR data processing and applications based on virtual voxel modeling theory. Attached Figure Description

[0040] Figure 1 A flowchart illustrating the airborne LIDAR power line extraction method based on target spatial distribution characteristics, provided for a specific embodiment of the present invention;

[0041] Figure 2 Raw airborne LiDAR data provided for specific embodiments of the present invention;

[0042] Figure 3 A flowchart for removing near-ground data is provided as a specific embodiment of the present invention;

[0043] Figure 4 The nDSM generation result diagram provided for a specific embodiment of the present invention;

[0044] Figure 5The "elevation-frequency" histogram in nDSM provided for specific embodiments of the present invention;

[0045] Figure 6 The image showing the results of near-ground data removal is provided for a specific embodiment of the present invention.

[0046] Figure 7 A flowchart for regularizing non-near-ground data into binary virtual 3D voxels is provided for a specific embodiment of the present invention;

[0047] Figure 8 A flowchart of the airborne LIDAR power line extraction steps provided for a specific embodiment of the present invention;

[0048] Figure 9 The present invention provides the results of airborne LIDAR power line extraction based on the spatial distribution characteristics of the target, which is provided for a specific embodiment of the present invention. Detailed Implementation

[0049] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0050] In this embodiment, the method of the present invention is implemented using MATLAB 2019a on an Intel(R) Core(TM) i5-7300HQ CPU@2.50GHz, 8GB of memory, and Windows 10 Ultimate system. The effectiveness of the method is further verified by evaluating the accuracy of the method.

[0051] like Figure 1 As shown, the method of this embodiment is described below.

[0052] Step 1: Read the raw airborne LiDAR data to form the raw airborne LiDAR dataset.

[0053] In this embodiment, the overhead power lines located in Sichuan Province were obtained using the NS-300RT PRO airborne laser scanner (at an altitude of 100 meters) from Chengdu Hangwei Zhixin Co., Ltd. as experimental data. Figure 2 As shown, the effectiveness and feasibility of the method were verified. The experimental data included low-lying vegetation, a small amount of tall vegetation, residential areas surrounded by crops, and overhead power lines. The data density was approximately 50 points / m². 2 .

[0054] In this embodiment, the original airborne LIDAR dataset P = {p i (x i ,y i ,z ip), i = 1, ..., n}, where i is the index of the original airborne LiDAR data, n is the number of original airborne LiDAR data, and p i This is the i-th raw airborne LiDAR data, with coordinates (x, y). i y i , z i ).

[0055] Step 2: Remove near-ground data from the raw airborne LiDAR data to obtain the non-near-ground dataset. The process for removing near-ground data is as follows: Figure 3 As shown.

[0056] Step 2.1: Separate ground points and non-ground points in the original data based on the cloth simulation filtering algorithm.

[0057] Step 2.2: Perform elevation interpolation on the original airborne LiDAR data and the separated ground data respectively to obtain a Digital Surface Model (DSM) and a Digital Elevation Model (DEM). Further, subtract the DEM from the DSM to generate a normalized Digital Surface Model (nDSM). The nDSM generation result is shown below. Figure 4 As shown.

[0058] In this embodiment, the DSM elevation value obtained based on the original airborne LiDAR data is denoted as {z}. i' ,i'=1,...,n'}, where i' is the index of the elevation value in the DSM; z i' Let be the i'-th elevation value in the DSM; n' be the number of elevation values ​​in the DSM. The DEM elevation value obtained based on ground data is denoted as {z}. o The elevation values ​​of the nDEM grid are denoted as {Δz, 0, ..., h}, where o is the index of the DEM grid elevation value; zo is the lowest elevation value of the o-th grid in the DEM; and h is the number of DEM grid elevation values. i' ,i'=1,...,n'}, where i' is the index of the elevation value in nDSM; Δz i' Let be the i'th elevation value in nDSM; n' is the number of elevation values ​​in nDSM.

[0059] The elevation value of nDSM is calculated by equation (1). Let the i'th elevation value in DSM be located in the oth DEM grid of DEM, and the lowest elevation value of the oth DEM grid be z. o Then the i'th elevation value Δz in nDSM i' It can be determined by equation (2).

[0060] nDSM = DSM - DEM (1)

[0061] Δz i' =z i' -z o (2)

[0062] Step 2.3: Based on the spatial distribution characteristics of elevation values ​​in the nDSM, an elevation threshold H0 is dynamically obtained using the "elevation-frequency" histogram analysis method. nDSMs below the elevation threshold H0 are classified as near-ground nDSMs, while those above are classified as non-near-ground nDSMs. Based on the mapping relationship between laser points and grids, the laser points corresponding to the non-near-ground nDSMs in each grid are read to obtain the non-near-ground dataset. The specific method is as follows:

[0063] Step 2.3.1: Count the frequency of elevation values ​​in nDSM and visualize the statistical results in the form of a histogram. Record the elevation value mutation point as the elevation threshold H0 of nDSM.

[0064] In this embodiment, the elevation threshold H0 is a constant, and its value needs to be determined based on the spatial distribution of elevation values ​​in nDSM. The optimal elevation threshold is determined by the following scheme (taking experimental data as an example):

[0065] The elevation value {Δz} in nDSM i' ,i'=1,…,n'} according to the interval width d Δz Within the elevation range Δz min ~Δz max If divided within, then the maximum interval level b max As shown in equation (3):

[0066]

[0067] Where, Δz min Δz max These are the minimum and maximum elevation values ​​in nDSM, respectively.

[0068] For the interval level b after the i'th elevation value in nDSM is divided i' As shown in equation (4):

[0069]

[0070] The frequency of elevations within each interval level after elevation division in nDSM is statistically analyzed, and the statistical results are displayed using a histogram, such as... Figure 5 As shown, the elevation value corresponding to the marked elevation change point is the elevation threshold H0 = 9.63m.

[0071] Step 2.3.2: Determine whether an nDSM below the elevation threshold H0 is a near-ground nDSM, and otherwise a non-near-ground nDSM.

[0072] Step 2.3.3: Based on the mapping relationship between laser points and grids in nDSM, read the laser points of the grids containing non-near-ground data in each grid, thus obtaining the non-near-ground dataset, such as... Figure 6 As shown.

[0073] In this embodiment, the non-near-ground dataset is denoted as G = {g j (x j ,y j ,z j ),j=1,...,t}, where j is the index of the laser point in the non-near-ground dataset G; g j For the j-th point in the non-near-ground dataset, the coordinates are (x... j ,y j ,z j ); t is the number of laser points in the non-near-ground data set.

[0074] Step 3: Regularize the non-near-ground data into a virtual 3D voxel structure, assign a value of 1 to the virtual voxel containing the laser point, and assign a value of 0 to other virtual voxels, to obtain a virtual binary 3D voxel. The specific process is as follows: Figure 7 As shown.

[0075] Step 3.1: Represent the 3D spatial extent using directed bounding boxes from the non-near-ground dataset.

[0076] In this embodiment, the directed bounding box of the non-near-ground dataset G is a cuboid. The size of its base is determined by finding the minimum bounding rectangle of the projection of the non-near-ground dataset G onto the XY plane, and its height is determined by (z...). max -z min ) Determined. Wherein, z max It is the maximum value of the z-coordinate of the laser point in the non-near-ground dataset G, z min It is the minimum z-coordinate of the laser point in the non-near-ground dataset G, z max =max{z j ,j=1,...,t},z min =min{z j ,j=1,...,t}.

[0077] Step 3.2: Determine the resolution (Δx, Δy, Δz) of the virtual voxel in the x, y, and z directions based on the average point spacing of the laser points in the non-near-ground data set. Δx, Δy, and Δz are the sizes of the virtual voxels.

[0078] In this embodiment, the formulas for calculating the resolutions Δx, Δy, and Δz of the voxel in the x, y, and z directions are shown in equation (5):

[0079]

[0080] Among them, S xy ={(x j ,y j Let {j = 1, ..., t} be a two-dimensional dataset obtained by projecting the non-near-ground dataset Q onto the XOY plane, and let C(S) be the dataset. xy ) is the dataset S xy A(C(S) convex shell xy )) is a convex shell C(S) xy )area.

[0081] Step 3.3: Based on the virtual voxel resolution (Δx, Δy, Δz), the directed bounding box is virtually divided to obtain a virtual 3D voxel grid. Each virtual 3D voxel grid unit is a virtual 3D voxel.

[0082] In this embodiment, the directed bounding box can be divided into a 3D voxel grid based on the voxel resolution (Δx, Δy, Δz), which is represented by a 3D voxel array. Let V be the set of voxels in the 3D voxel array, as shown in equation (6):

[0083] V = {v} s (r s ,c s ,l s ),s=1,...,k} (6)

[0084] Where s is the volume element index; k is the number of volumes; v s It is the voxel value of the s-th voxel; (r s ,c s ,l s ) represents the coordinates (row, column, and layer number) of the s-th voxel in the voxel array. The number of voxels in the X direction is R, the number of voxels in the Y direction is C, and the number of voxels in the Z direction is L. R, C, and L are given by formula (7):

[0085]

[0086] in, The rounding up operator, x max =max{x j ,j=1,...,t},x min =min{x j ,j=1,...,t},y max =max{y j ,j=1,...,t},ymin =min{y j ,j=1,...,t}.

[0087] Therefore, the number of volume elements k is as shown in equation (8):

[0088] k=R*C*L (8)

[0089] Step 3.4: Map each laser point in the non-near-ground dataset to a virtual 3D voxel, then assign the value 1 to the virtual 3D voxel containing the laser point and assign the value 0 to the other virtual 3D voxels to obtain a virtual binary 3D voxel.

[0090] In this embodiment, each laser point in the non-near-ground dataset G is mapped to a 3D voxel grid, and then each voxel is assigned a value based on whether it contains a laser point. Specifically, virtual voxels containing laser points are assigned a value of 1, and other virtual voxels are assigned a value of 0, as shown in equation (9):

[0091]

[0092] in, This is the floor operation. Thus, a virtual binary 3D voxel is obtained, completing the regularization of non-near-ground data.

[0093] Step 4: Based on the spatial distribution characteristics of the target, principal component analysis is used to determine the distribution shape of each 1-valued virtual voxel and the target points in its spatial neighborhood. Virtual voxels with linear distributions are marked, and laser points within these linear virtual voxels are identified as electric field lines. This completes the extraction of electric field lines based on the spatial distribution characteristics of the target. The specific process is as follows: Figure 8 As shown.

[0094] Step 4.1: Initialize the laser point closest to the center within each 1-value virtual voxel as the target point for each 1-value virtual voxel. The specific method is as follows:

[0095] Step 4.1.1: Calculate the distance d between all laser points within each 1-value virtual voxel and the center of each 1-value virtual voxel. c .

[0096] In this embodiment, based on the row, column, and layer coordinates (r) of the s'th 1-valued virtual voxel in array V... s’ ,c s’ ,l s’ The center coordinates (x, y) of the virtual volume element with a value of 1. s’c ,y s’c ,z s’c As shown in equation (10):

[0097]

[0098] Step 4.1.2: Select the distance from the center of each virtual voxel with a value of 1 as the nearest distance d. c,min The laser point is used as the target point of each 1-value virtual volume element.

[0099] In this embodiment, based on the number of laser points within each 1-value virtual voxel, the 1-value virtual voxels can be divided into two categories: those containing only one laser point and those containing multiple laser points. For a 1-value virtual voxel containing only one laser point, that laser point is the target point within the 1-value virtual voxel. For a 1-value virtual voxel containing multiple laser points, it is necessary to find the laser point closest to the center of the 1-value virtual voxel as the target point of that voxel. Let the point set within the s'-th 1-value virtual voxel be denoted as G'={g j' (x j' ,y j' ,z j' ),j'=1,...,t'}, where j' is the index of the laser point in the point set G'; g j' Let j' be the j-th point in the corresponding point set, with coordinates (x, y, y). j' ,y j' ,z j' ); t' is the number of laser points in the point set. Then the distance between the laser point in the s'th virtual volume element with a value of 1 and the center of the volume element is calculated as shown in equation (11):

[0100]

[0101] Step 4.2: Based on principal component analysis, construct linear, planar, and spherical distribution indices L to describe the spatial distribution characteristics of target points within each 1-valued virtual volume element. λ P λ and S λ And a direction vector describing the directional features of the target point of each 1-valued virtual volume element. The specific method is as follows:

[0102] Step 4.2.1: Scan each 1-value virtual element and perform a spatial neighborhood search to obtain its neighborhood set.

[0103] In this embodiment, applying different spatial neighborhood scales will result in different neighborhood sets, which in turn will affect the subsequent extraction results of electric field lines. The optimal neighborhood scale will be determined experimentally.

[0104] Step 4.2.2: Construct the covariance matrix Cov based on each 1-valued virtual voxel target point and its neighborhood target point set, then construct a homogeneous equation, and solve for the eigenvalues ​​λ1, λ2, λ3 (λ1≥λ2≥λ3) of the covariance matrix Cov and the corresponding eigenvectors. Among them, the eigenvector corresponding to the largest eigenvalue λ1 Let be the direction vector of each 1-valued virtual volume element target point.

[0105] In this embodiment, the target point within any 1-valued virtual voxel and the set of target points within its m-neighborhood are denoted as Q = {q}. u (x u ,y u ,z u ), u=1,2,...,m}, where u is the index of the target point within the 1-valued virtual volume element; q u Let the target point be the u-th virtual volume element with a value of 1, and its coordinates be (x, y, y). u ,y u ,z u ); m is the number of 1-valued virtual elements in the neighborhood set. The covariance matrix Cov is shown in equation (12):

[0106]

[0107] in, Let Q be the 3D centroid of the neighborhood set.

[0108] A homogeneous equation is constructed and solved for the covariance matrix Cov, as shown in equation (13):

[0109]

[0110] Step 4.2.3: Based on the feature values ​​λ1, λ2, and λ3 of each of the above-mentioned one-valued virtual volumetric target points, construct a linear distribution index L describing the distribution characteristics of each one-valued virtual volumetric target point. λ The areal distribution index P λ and spherical distribution index S λ .

[0111] In this embodiment, the linear, planar, and spherical distribution indices are respectively represented by L. λ S λ P λ This is represented as shown in equation (14):

[0112]

[0113] Among them, L λ S λ P λ ∈[0,1],L λ +S λ +P λ =1.

[0114] Step 4.3: Determine the linear distribution index L λ Greater than the areal distribution index P λ and spherical distribution index S λ And direction vector The target point of the 1-value virtual volume element that is approximately parallel to the horizontal plane is the seed point of the electric field line. The 1-value virtual volume element containing the seed point of the electric field line is marked as a linear virtual volume element. The laser points in the linear virtual volume element are read as electric field line data to complete the extraction of the electric field line.

[0115] Step 4.3.1: For each target point within a 1-value virtual voxel, determine the linear distribution index L. λ At the same time, it is greater than the areal distribution index P. λ and spherical distribution index S λ And direction vector The target point within a 1-valued virtual element that is approximately parallel to the horizontal plane is the seed point for the electric field line.

[0116] In this embodiment, the target point type is determined based on the target point feature distribution index within each 1-value virtual volume element, as shown in Table 1:

[0117] Table 1 Target Point Type Identification

[0118] <![CDATA[λ1>>λ2≈λ3]]> <![CDATA[L λ >P λ And L λ >S λ ]]> linear <![CDATA[λ1≈λ2>>λ3]]> <![CDATA[P λ >L λ And P λ >S λ ]]> planar <![CDATA[λ1≈λ2≈λ3]]> <![CDATA[S λ >L λ And S λ >P λ ]]> spherical

[0119] For the direction vector of each 1-valued virtual volume element target point The criterion for determining whether a direction vector is approximately parallel to the horizontal plane is as follows: if the direction vector is approximately parallel to the horizontal plane... unit vector with x-axis The closer the included angle θ is to 0°, the more likely the principal direction of the virtual volume element target point with a value of 1 is to be approximately parallel to the horizontal plane. Considering the catenary-like characteristics of electric field lines, the direction vector... unit vector with x-axis The included angle θ should be greater than or equal to 0° and less than the threshold T, as shown in equation (15):

[0120]

[0121] in, These are the direction vector and the x-axis unit vector of the target point of the 1-valued virtual volume element, respectively, where θ is... and The included angle. The threshold T is a constant; in this embodiment, T = 10°.

[0122] Step 4.3.2: Mark the 1-value virtual element where the power line seed point is located as the power line virtual element, and read the laser points in the power line virtual element as power line data to complete the extraction of the power line.

[0123] The following uses manually extracted power lines obtained from the domestic LiDAR power line inspection data processing software "Cloud Compare v2.12" as reference data to quantitatively evaluate the power line extraction accuracy of the method of this invention. The number of power line data points extracted by the method of this invention is compared with the reference data. The effectiveness of the proposed power line extraction method is then quantitatively evaluated using Type I error (the proportion of power line laser points misclassified as non-power line laser points), Type II error (the proportion of non-power line laser points misclassified as power line laser points), total error (the proportion of misclassified power line laser points), accuracy (the proportion of correctly extracted power line laser points to the total number of power line laser points in the extraction result), completeness (the proportion of correctly extracted power line laser points to the total number of power line laser points in the reference data), and quality.

[0124] Table 2 shows the accuracy of electric field line extraction results for experimental data at neighborhood scales of 18, 26, 56, 80, and 124. The data in this table aim to examine the impact of different neighborhood scales on the electric field line extraction results and thereby determine the optimal neighborhood scale.

[0125] Table 2. Accuracy index of power line extraction results at different neighborhood scales

[0126]

[0127] Table 2 shows that the quality of the method of this invention in 18, 26, 56, 82, and 124 neighborhoods is 23.23%, 57.46%, 88.17%, 38.90%, and 98.65%, respectively. This indicates that the electric field line extraction has high quality in the 80-neighborhood, making it the optimal neighborhood scale for electric field line extraction. However, increasing the neighborhood scale does not necessarily improve extraction accuracy. Taking 18 and 26 neighborhoods as examples, their smaller neighborhood scale makes them more susceptible to uneven data distribution. If there are too many consecutive blank voxels in the virtual voxels, the spatial neighborhood voxel set cannot be obtained, meaning spatial shape feature calculation is impossible. This leads to incomplete electric field line extraction results, as shown in Table 2 where the completeness of the 18-neighborhood and 26-neighborhood algorithms is significantly lower. As the neighborhood scale increases, more virtual voxels are incorporated into the spatial neighborhood of the target voxel, reducing the probability of Type I errors. This explains why 56 and 80 neighborhoods produce better extraction results than 18 and 26 neighborhoods. However, if the neighborhood scale is too large, too many neighboring voxels are incorporated into the spatial neighborhood of the target voxel, leading to an increased probability of Type II errors. As shown in Table 2, the accuracy of 124 neighborhoods significantly decreases compared to 80 neighborhoods, resulting in a decrease in algorithm accuracy. This explains why the accuracy of 124 neighborhoods is actually lower than that of 80 neighborhoods.

[0128] The electric field line extraction results obtained by applying the method of this invention are as follows: Figure 9 As shown in the figure. The size of a single virtual volume element is 0.4m×0.4m×0.3m. The number of electric field laser points and non-electric field laser points extracted using the optimal neighborhood scale of 80 neighborhoods are 64926 and 20566, respectively.

[0129] Table 3 presents a quantitative evaluation of the power line extraction accuracy at the 80-neighborhood scale of the experimental data, using the reference data as a standard.

[0130] Table 3. Accuracy of power line extraction results

[0131] 99.08 99.82 98.90

[0132] As shown in Table 3, the completeness, accuracy, and quality of the extracted power lines were 99.08%, 99.82%, and 98.90%, respectively. This verifies the effectiveness of the method proposed in this invention.

[0133] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. A method for extracting power lines using airborne LiDAR based on the spatial distribution characteristics of a target, characterized in that: The method includes the following steps: Step 1: Read the raw airborne LiDAR data to form the raw airborne LiDAR dataset; Step 2: Remove near-ground data from the original airborne LiDAR data to obtain the non-near-ground dataset; Step 3: Regularize the non-near-ground data into a virtual 3D voxel structure, assign a value of 1 to the virtual voxel containing the laser point, and assign a value of 0 to other virtual voxels, to obtain a virtual binary 3D voxel; specifically including the following steps: Step 3.1: Represent the 3D spatial extent using directed bounding boxes from the non-near-ground dataset; Step 3.2: Determine the virtual volume element based on the average point spacing of the laser points in the non-near-ground data set. x , y , z Resolution in direction , That is, the size of the virtual volume element; Step 3.3: Based on virtual voxel resolution The directed bounding box is virtually divided to obtain a virtual 3D volumetric grid, and each virtual 3D volumetric grid cell is a virtual 3D volumetric element. Step 3.4: Map each laser point in the non-near-ground data set to a virtual 3D voxel, and then assign the value 1 to the virtual 3D voxel containing the laser point and the value 0 to the other virtual 3D voxels, as shown in equation (9). （9） in, Voxels in a 3D voxel array. V It is a collection of voxels in a 3D voxel array. , s It is a volume element index; k It is the number of volume elements; It is the first s The voxel value of an individual element; r s , c s , l s ) is the first s The coordinates of an individual element in the voxel array r s 、c s 、l s These represent the row number, column number, and layer number, respectively. This is the floor operator; For the corresponding non-near-ground dataset, the first j The coordinates of the points; , , , t This represents the number of laser points concentrated in the non-near-ground data set. Thus, virtual binary 3D voxels are obtained, completing the regularization of non-near-ground data; Step 4: Based on the spatial distribution characteristics of the target, use principal component analysis to determine the distribution shape of each 1-valued virtual voxel and the target points in its spatial neighborhood. Mark the virtual voxels that are distributed in a linear pattern, and identify the laser points in the linear virtual voxels as electric field lines, thus completing the extraction of electric field lines based on the spatial distribution characteristics of the target.

2. The airborne LIDAR power line extraction method based on target spatial distribution characteristics according to claim 1, characterized in that: Step 2 includes the following steps: Step 2.1: Separate ground points and non-ground points in the original airborne LiDAR data based on the cloth simulation filtering algorithm; Step 2.2: Interpolate the original airborne LiDAR data and the separated ground data respectively to obtain the digital surface model (DSM) and digital elevation model (DEM), and further subtract the DSM from the DEM to generate a regularized digital surface model (nDSM). Step 2.3: Based on the spatial distribution characteristics of elevation values ​​in nDSM, dynamically obtain the elevation threshold using the "elevation-frequency" histogram analysis method. H 0, and will be below the elevation threshold. H An nDSM with a value of 0 is identified as a near-ground nDSM, while one with a value of 0 is identified as a non-near-ground nDSM. Based on the mapping relationship between laser points and grids, the laser points corresponding to the non-near-ground nDSMs in each grid are read to obtain the non-near-ground dataset.

3. The airborne LIDAR power line extraction method based on target spatial distribution characteristics according to claim 2, characterized in that: Step 2.3 specifically includes the following steps: Step 2.3.1: Count the frequency of elevation values ​​in nDSM and visualize the statistical results in the form of a histogram. Record the elevation value abrupt change points as the elevation threshold of nDSM. H 0; Step 2.3.2: For the elevation values ​​in nDSM, determine those below the elevation threshold in nDSM. H Data with a value of 0 represents near-surface nDSM data, while data with a value of 0 represents non-near-surface nDSM data. Step 2.3.3: Based on the mapping relationship between laser points and grids in nDSM, read the laser points of the grids containing non-near-ground data in each grid to obtain the non-near-ground dataset.

4. The airborne LIDAR power line extraction method based on target spatial distribution characteristics according to claim 1, characterized in that: Step 4 includes the following steps: Step 4.1: Initialize the laser point closest to the center of each 1-value virtual voxel as the target point of each 1-value virtual voxel; Step 4.2: Based on principal component analysis, construct a linear distribution index describing the spatial distribution characteristics of the target points of each 1-valued virtual volume element. areal distribution index spherical distribution index And the direction vector describing the directional features of the target point of each 1-valued virtual volume element. ; Step 4.3: Determine the linear distribution index Greater than the areal distribution index and spherical distribution index And direction vector The target point of the 1-value virtual volume element that is approximately parallel to the horizontal plane is the seed point of the electric field line. The 1-value virtual volume element containing the seed point of the electric field line is marked as the electric field line virtual volume element. The laser points in the electric field line virtual volume element are read as the electric field line data to complete the extraction of the electric field line.

5. The airborne LIDAR power line extraction method based on target spatial distribution characteristics according to claim 4, characterized in that: Step 4.1 specifically includes the following steps: Step 4.1.1: Calculate the distances between all laser points within each 1-value virtual voxel and the center of each 1-value virtual voxel. d c ; Step 4.1.2: Select the distance between each 1-value virtual voxel and its center as the shortest distance. d c,min The laser point is used as the target point of each 1-value virtual volume element.

6. The airborne LIDAR power line extraction method based on target spatial distribution characteristics according to claim 4, characterized in that: Step 4.2 specifically includes the following steps: Step 4.2.1: Scan each 1-value virtual volume element, perform a spatial neighborhood search on it, and obtain its neighborhood set; Step 4.2.2: Construct the covariance matrix based on each 1-valued virtual volume element target point and its neighborhood target point set. Cov Then, a homogeneous equation is constructed, and the covariance matrix is ​​solved. Cov eigenvalues λ 1. λ 2. λ 3 and the eigenvectors corresponding to the eigenvalues ​​respectively , , ,in, λ 1≥ λ 2≥ λ 3. Largest eigenvalue λ eigenvector corresponding to 1 For each 1-valued virtual volume element target point, there is a direction vector. Step 4.2.3: Based on the feature values ​​of each of the above-mentioned 1-value virtual volume element target points λ 1. λ 2. λ 3. Construct a linear distribution index to describe the distribution characteristics of target points in each 1-valued virtual volume element. areal distribution index and spherical distribution index As shown in equation (14): （14） in, , .

Images