AUTOMATIC DETECTION AND CHARACTERISATION OF POWER LINES AND THEIR SURROUNDINGS USING LIDAR DATA

Abstract. Light Detection and Ranging (LiDAR) is nowadays one of the most used tools to obtain geospatial data. In this paper, a method to detect and characterise power lines of both high and low voltage and their surroundings from 3D LiDAR point clouds exclusively is proposed. First, to identify points of the power lines a global search of candidate points is carried out based on the height of each point compared to its neighbours. Then, the Hough Transform (HT) is applied on the set of candidate points to extract the catenaries that belong to each power line, allowing the identification of each conductor individually. Finally, conductors located on the same power line are grouped, their geometric characteristics analysed, and the quantitative features of the surroundings are computed. A very high accuracy of power line classification is reached with these methods, while the computational time is optimised by efficient memory usage and parallel implementation of the code.



INTRODUCTION
Today, most of the activities of modern-day societies depend on electricity. Thus, it is necessary to ensure the efficient operation of the electric transmission networks, which typically include nationwide, regional and distribution networks. In forested countries, most parts of the networks are located inside forests (Matikainen et al., 2016), and they must be checked and maintained periodically. The absence of obstacles in the transmission line corridor, i.e. Right-of-Ways (ROWS), is a security requirement that must be fulfilled. The process of checking this requirement is traditionally performed through visual inspection of the elements of the power line, which can be expensive, error prone, and even risky, because the observation areas are not always easily accessible. Thus, methods based on airborne LiDAR data are proposed in order to characterise the power lines and its surroundings with high precision, and therefore to detect dangerous obstacles that may interfere with the proper power lines operation.
In the literature a number of proposal can be found. (McLaughlin, 2006) proposes the segmentation of the LiDAR point cloud in elliptic neighbourhoods and the analysis of the covariance matrix eigenvalues. This allows to distinguish among three types of classifications: power lines, vegetation and buildings. Furthermore, the reconstruction of the power lines in order to obtain the parameters of each transmission line is proposed.
In (Zhu , Hyyppä, 2014) a method to power line extraction in forest areas using a rasterized point cloud is proposed. This method is based, on the one hand, on the statistical analysis to identify power line candidates, using height and point density criteria. On the other hand, candidate points are converted into a binary image, where image processing techniques are applied to find power lines.
In (Clode , Rottensteiner, 2005) a method for tree and power line extraction in urban areas is shown, the height difference between the first and the last pulse is used and finally the intensity of LiDAR points is taken into account. The distance between trees and power lines is computed by applying a classification method based on the theory of Demper-Shafer.
A different method for power line extraction is shown in (Wang et al., 2017). A DTM is built to select candidate points. Then, multi-scalar spherical neighbourhoods are used to group the candidates into clusters. Inside the clusters, feature extraction is performed in order to distinguish between power lines and other objects. Finally, the authors classify line points using support vector machine, which needs a training data set to calibrate the classifier. (Liu et al., 2009) carry out an analysis of the skewness and kurtosis in order to classify LiDAR points as ground and non-ground. They employ an improved Hough Transform over a rasterization of the points to detect power lines in the grayscale 2D image.
Finally, (Sohn et al., 2012) perform the segmentation of power lines using Markov Random Fields. Once the segmentation is complete, a voxel based rasterization is carried out. For each voxel, they apply the Hough Transform. The linear candidate points are then converted into line segments by 3D line fitting based on Random Sample Consensus (RANSAC). To detect the pylons, line segments are converted into a binary image and Random Forests are used to perform a binary classification of pylon vs no pylon.
In this paper we present a novel power line classification method from airborne LiDAR that relies only on the spatial coordinates of the points. Neither intensity data nor number of returns are needed to perform a successful power line classification. The algorithms presented in this paper work at point level, so neither a 2D binary image nor a rasterization is necessary. Furthermore, the methods to be exposed in this paper can be immediately applied to any data set, since a large data set to train any supervised classification method is not needed. The method is robust in the sense that it works for non specific flights along the power lines. In fact, all the use cases in this paper are obtained from general purpose LiDAR data. Additionally, we consider that the computational costs needed to perform the classification are very low, despite only a few authors provide this information. We deal with two techniques that improve computational performance by using efficient 3D data structures that take advantage from the memory locality of data, and by using OpenMP implementations to parallelize the code for multi core systems. These solutions decrease significantly the execution times of our proposal. Experiments show that an average correctness of 99, 24% was achieved, with an execution time less than two minutes for a 8 · 10 6 points data set.

POWER LINE EXTRACTION
The two main elements of the power lines are wires and pylons. According to our proposal, the first part of the power line detection consists in the extraction of the individual conductors. Detection of pylons is performed once the wires are detected, since the position of them depends on the position and orientation of the wires.

Iterative Candidate Search
The objective of the Iterative Candidate Search (ICS) is to label every LiDAR point that could be part of a conductor. This stage also allows us to reduce the computational costs of the subsequent steps, since the algorithms will be run on a reduced amount of points compared to the original size of the point cloud.
The ICS is based on two premises: Firstly, transmission lines are elevated over the ground and secondly, there are no structures above and near the transmission lines. These premises can be checked by just using the spatial information of the LiDAR points, so a previous classification of the data is not needed. The proposed algorithm computes the percentage of neighbour points (C) that are located at a vertical distance of the current analysed point. If the percentage is big enough, the target point is labelled as a power line candidate.
Let pi = (pix, piy, piz) be the current LiDAR point defined only by its spatial coordinates. The point cloud is defined as the set of all points obtained in a specific flight, i.e. Ω = {p1, p2, . . . , pn}. Thus, for each LiDAR point pi, the subset Ni = {n1, n2, . . . , nm} ⊂ Ω / (pix − Rs) ≤ njx ≤ (pix + Rs), (piy − Rs) ≤ njy ≤ (piy + Rs) is built. This set consists of the neighbours of the currently analysed LiDAR point, where Rs is the size of the search distance around each LiDAR point.
Once the subset Ni is built, the algorithm performs two calculations: On the one hand it obtains the number of neighbours whose vertical distance to the currently analysed point is large enough to not belong to the same wire. To achieve this, the Wire Thickness Threshold (W th ) is defined as the estimated thickness of a conductor. The number of points Nnwn meeting this condition is given by: where Ni = set of neighbours of pi niz = height of each neighbour piz = height of the currently analysed point W th = wire thickness | · | = the cardinality of a set of points On the other hand, the algorithm computes the number of neighbours whose vertical distance is greater than H th , where H th is the Height Threshold parameter.
where Npos = number of points that meet the condition Ni = set of neighbours of pi niz = height of each neighbour piz = height of the currently analysed point H th = height threshold The points that eventually can influence these computations, are depicted in Figure 1. The blue point is the current analysed LiDAR point (pi). Npos is the number of green points and Nnwn is the sum of green and red points. Once the values Nnwn and Npos are computed for each point, the descriptor Ci is defined as It may be that for some point, Nnwn = 0. This will eventually happen on points in the ground. In that case, we assume that Nnwn = |Ni| where |Ni| is the cardinality of the set Ni, i.e. the number of neighbours no matter its height.
Finally, we label the points as members of the wire if they meet the following condition: where C th is the minimum ratio of points that must meet the conditions imposed by Equations 1 and 2. This algorithm was firstly designed assuming the absence of points above the wires, so in the case of several wires located on top of one another, the algorithm will struggle to correctly label the points belonging to the lower wires. This issue is solved by executing the candidate search iteratively, ignoring on each iteration the points previously labelled as wire members. The execution finishes when there are no more candidates to consider. The results of each iteration of the ICS are shown in Figure 2.

Power Line Detection
After the ICS stage, the 2D Hough Transform (HT) proposed in (Duda , Hart, 1972) is applied. A projection of the LiDAR points over the XOY plane were performed in order to detect straight lines.
On the Hough Space, each straight line is represented by where x, y = spatial coordinates of each LiDAR point ρ = distance from the origin of coordinates θ = angle between the x-axis and the normal to the line The sensitivity of the HT can be adjusted through the modification of three parameters: The angle step (As) that is the angle separation between two consecutive lines; The grid spacing (Gs) that is the minimum distance between lines that HT can distinguish. And (Nc) that is the minimum number of points that a line must have to be detected.
On each iteration, HT finds the straight line with the greatest number of points, and gives the polar coordinates of that line, i.e., (ρ, θ).
In order to identify points belonging to each line, at the end of each iteration a reversed voting is carried out. Knowing the pair (ρ l , θ l ) of the detected line, points whose spatial coordinates (x, y) fits into Equation 5 are searched. As a result of this step, each conductor of each power line is grouped in one specific cluster of points.
To identify every straight line in the LiDAR point cloud, HT is executed iteratively until the last line found contains at least Nc points. At the end of each iteration, the votes to the detected line are checked and the corresponding points are grouped in clusters.

Splitting Vertical Conductors
Due to the bi-dimensional nature of the HT, conductors located in the same vertical plane are detected as an unique line. A method to split the line into individual clusters, each one corresponding to one conductor was developed.
The first step is to represent the different wires contained in one plane, therefore a 2D representation of the lines is introduced. Each point is represented by its distance to a reference point and its z-coordinate, where the reference point is the one with the minimum x-coordinate. Thus, the 2D plot will show the height (H-axis) versus the distance to the first point (D-axis). An example of a triple wire configuration is shown in Figure 3.  Once the representation is computed, the next step is to identify the number of separated catenaries co-existing in the same vertical plane. The plots of each line are rasterized in order to turn each scatter plot into a 2D histogram, as shown in Figure 4. As the number of points in each case can be different, the number of bins in the histogram is chosen to be n bin × n bin , where n bin = √ N , according to the square-root choice method and N is the number of points of the analysed conductors.
When the histogram is built, the point density is computed. In order to avoid the effects of spurious points, we obtain two point densities: one using the total number of bins, and another using the number of non empty bins, as shown in Equation 6: where ρ = point density ρne = point density of the non-empty bins N = number of points of the power line n bin = number of total bins n nebin = number of non empty bins For each column in the 2D histogram, the number of adjacent bins with a difference in the number of points greater than ρ is computed, as shown in Equation 7: (7) where nj = number of detected gaps i = position of the bin bin = array of each column values of the 2D histogram The results of applying the Equation 7 to two columns of the 2D-histogram previously shown are displayed in Figure 5.
When the value of nj is obtained for each column in the histogram, we assume that the number of catenaries placed in the same vertical plane is equal to the mode of all the values computed of nj.
Due to the lack of power lines with enough spurious points in the available examples, we employ a synthetic test to test the robustness of the algorithm. A set of random catenaries with spurious points are generated, as shown in Figure 6. Each catenary consists of 300 points, and a dispersion both in distance and height is randomly generated in order to approximate the behaviour of real data sets. The number of generated catenaries is 3, and their extremes are 6 meters apart. For a given number of spurious points, the algorithm is run one hundred times, counting the number of successes and failures in the identification of the number of catenaries. The results are shown in Figure 7. The three catenaries were detected in every case study from 0 to 500 generated spurious points. In real data sets, find such amount of spurious points is not expected, so the algorithm can be considered robust.

Pylons Detection
Pylons are basic elements of power lines, as they are the structures that support and guide them. In this paper we propose a method for the detection of the areas where the pylons are located.
A wire hanging between two pylons under the action of its weight, adopts a well known shape named catenary. The proposed algorithm for pylon identification uses the first and second spatial derivatives of the wire points to locate discontinuities. At those discontinuities is where pylons are expected to be found. Furthermore, the existence of a pylon at the start and the end of a transmission line is assumed.
Due to the non uniform and scattered distribution of the LiDAR points, the use of finite difference methods produce inaccurate results. The proposed method is based on the computation of linear fits around each point. In order to perform the linear fits, it is necessary to transform each catenary into a 2D scatter plot, using the same strategy introduced in Section 2.3. The The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W13, 2019 ISPRS Geospatial Week 2019, 10-14 June 2019, Enschede, The Netherlands   Figure 9b). The points where the first derivative is zero correspond to two areas: The local minimums of the catenaries and the discontinuities. Theoretically, points corresponding to the discontinuities of the catenaries will be those whose value of the slope is nearest to zero. Eventually, using these points to locate the pylons yields to inaccurate results, due to the unbalance of adjacent catenaries: the end of one catenary and the start of the adjacent one do not have to be symmetrical; they may have different density of points or slope magnitude. This causes that the zero slope positions in the first derivative are located around the peaks of the catenary.
In order to discriminate the local minimums from the discontinuities, the second derivative is analysed, as shown in Figure 9c. The position of the discontinuities on the catenaries, and the position of the pylons, matches with the mid point of the valleys of the second spatial derivative.

Power line characterisation
Once the location of the pylons is known and the catenaries are individually identified, it is possible to extract the geometric characterisation of the power lines.
The catenary curve can be modelled as shown in Equation 8: where a = vertical displacement from the X-axis (m) b = horizontal displacement from the Y-axis (m) c = vertical position of the catenary vertex (m) The parameter c also corresponds to the ratio between the horizontal tension of the wire at its lowest point and the linear weight density: In order to compute the parameters a, b and c, a non linear fit of the data to Equation 8 is carried out.
Furthermore, using the catenary curve equation it is possible to compute the length of each conductor. The length of any given function between two points is given by: where L = computed length of the catenary p1 = initial horizontal position of the catenary p2 = end horizontal position of the catenary Substituting Equation 8 into 10 yields to Since we have individually clustered the LiDAR points of each conductor, it is possible to compute important distances in a very accurate way, by computing point-to-point: distances between conductors inside the same power line, width of the power line and distance between pylons. Furthermore, using the classification described in (Martínez et al., 2016), the distance to the closest vegetation, ground and building can be easily computed.

Experiments and Results
To validate the results obtained in the classification of power lines using the proposed algorithms, we conduct experiments over two point clouds data sets obtained from different sources: the point cloud corresponding to Rozas area was provided privately, while the point cloud corresponding to Diablo Canyon area was obtained from OpenTopography (OpenTopography, 2010). The values of the necessary parameters to run the algorithms are listed in Table 1 and the characteristics of the points clouds are listed in Table 2. The ground truth was acquired manually using the OLIVIA software (Martínez et al., 2018). The results of the classification are analysed attending to three parameters: correctness, completeness and quality, defined as shown in Equation 12:    The average completeness given in this paper is similar to those found in (Chen et al., 2018)  The correctness error is due to isolated points located at the top of high vegetation, and aligned with any power line. Completeness error is caused during the ICS stage, and it is due to those points belonging to transmission lines that do not pass this filter.
To increase the computational efficiency of the algorithms, the point clouds are organised into an octree structure, as shown in (Martínez et al., 2016). By using the octree, a higher neighbour points locality in memory is achieved, which improves the overall performance significantly. Also, both the ICS and the Hough Transform are coded to be executed in parallel using OpenMP. The following performance tests are carried out using the multi core system Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz. A common measure to test the quality of a parallel implementation is the speedup. The speedup is defined as the ratio between the execution time with n cores and with 1 core. The achieved speedup is shown in Figure 10. It presents an almost optimal linear behaviour when using the 4 physical cores. In the range of hyper-threading, from 4 to 8 threads, the speedup decreases, due to bottlenecks related with memory accesses.  Finally, Figure 11 displays the relationships between the number of points and the execution times. In addition, the relation between the number of detected points and the execution times of the ICS and HT are shown in Figure 12.
The tests were performed over reduced versions of the Diablo Canyon data sets. Note a linear growing of the execution time until 5 000 000 points, when it starts to grow faster. This behaviour is caused by the Hough transform and its high computational cost when the number of detected points is very high, as shown in Figure 12. The data sets used in the validation of the power line classification were acquired for general purposes, so we can consider that the obtained results in this paper are very competitive against other methods, since some authors like (Chen et al., 2018), (Guo et al., 2016) or (Liu et al., 2009) used airborne LiDAR data sets obtained specifically to detect transmission lines.

CONCLUSIONS
In this paper a general method to automatically detect and characterise power lines is shown. A fast algorithm to select power line candidates based on parameters was developed. The proposed iterative candidate search yielded very good results discriminating power line points from ground, buildings and vegetation. After the selection of candidates, the power lines were identified by applying the Hough transform, and a method to split different conductors in each detected line was developed. Finally, the position of the pylons was determined using geometric tools.
A very high degree of correctness and precision was achieved while keeping the computational costs very low. By applying our algorithms, an average quality of 93.84% was achieved in the power line detection. Our method has several advantages: specific flights over power lines are no needed, since the algorithms can be applied over any data set; not only for classifying power lines, but for detecting its presence. Furthermore, the algorithms are not restricted to a specific type of environment: they yield good results in open areas, urban zones and forest areas. Finally, it's possible to compute distances directly from point to point, which increases significantly the accuracy of the measures.

ACKNOWLEDGES
This work has received financial support from Babcock International Group, in the frame of the Civil UAVs Initiative of Xunta de Galicia. Also, it received support from the Consellería de Cultura, Educación e Ordenación Universitaria of Xunta de Galicia (accreditation 2016-2019, ED431G/08 and reference competitive group 2019-2021, ED431C 2018/19) and the European Regional Development Fund (ERDF). In addition, this research was also funded by the Ministerio de Economía, Industria y Competitividad within the project TIN2016-76373-P.