Natural gas micro-leakage stress vegetation detection method based on LiDAR data

Abstract

The invention relates to a natural gas micro-leakage stress vegetation detection method based on LiDAR data, and the method comprises the steps: carrying out the rasterization of vegetation point cloud data, so as to extract the height, homogeneity, contrast ratio and other stereo structure features, carrying out the clustering analysis of the dimension-reduced LiDAR data through employing a Kmeans clustering algorithm, carrying out the classification of the clustering results under the constraint of a JM distance criterion, and carrying out the detection of the natural gas micro-leakage stress vegetation. And determining a natural gas stress vegetation stress range by using a morphological threshold segmentation method, and finally determining the position of the micro-leakage point of the underground gas storage according to the stress range, so that the micro-leakage point of the natural gas can be effectively detected.

Description

Technical Field

[0001] The present invention belongs to the technical field of natural gas leakage detection, and in particular relates to a natural gas micro-leakage stress vegetation detection method based on LiDAR data. Background Art

[0002] Natural gas is an environmentally friendly, low-carbon, clean fuel with a wide range of uses, and its demand and consumption continue to grow worldwide (Faramawy et al., 2016). Underground storage is the main method for storing natural gas, but there is a risk of leakage due to geological faults, human activities or other natural factors (Tarkowski et al., 2021; Zhang et al., 2021). Once a leak occurs in an underground storage, it will not only cause economic losses, but also cause safety and ecological environmental problems (Lukonge and Cao, 2020; Ren et al., 2018). Timely and accurate detection of micro-leakage points in natural gas storage can avoid serious economic and environmental consequences.

[0003] Traditional natural gas leak detection methods include manual inspection (Boaz et al., 2014), sound monitoring (Jin Het al., 2014), gas sampling (Sperl, 1991) and flow monitoring (Abdulshaheed et al., 2017). These methods consume a lot of manpower and time, are destructive, and cannot meet the requirements of large-scale rapid inspection. In addition, when a trace leak occurs in a natural gas storage reservoir, the natural gas concentration in the surface atmosphere is often difficult to reach the threshold of direct detection, and traditional methods are not suitable for natural gas micro-leakage (Y. Pan et al., 2022; Vodnik et al., 2006). However, leakage from underground natural gas storage facilities will change the soil environment in the area, displace oxygen in the soil, interfere with the respiration of vegetation roots, and cause surface vegetation to respond to natural gas stress. The aboveground part shows stress phenomena such as slow growth, sparse canopy, and reduced plant height (X. Pan et al., 2022; Y. Pan et al., 2022; Ran et al., 2020; Xiong et al., 2022). Remote sensing methods can be used to indirectly detect micro-leakage points of underground natural gas storage facilities in vegetation areas.

[0004] Hyperspectral technology can obtain continuous and detailed spectral information of the target object within a certain range. In recent years, it has been widely used to indirectly detect micro-leakage points of underground natural gas storage facilities in vegetation areas (Du et al., 2022; Jiang et al., 2023, 2020; Y. Pan et al., 2022; Ran et al., 2022). However, it can only reflect the spectral information of the vegetation surface and cannot capture the response of the three-dimensional structural characteristics of vegetation to natural gas stress. Due to its "same object, different spectrum" characteristics, the detection model that only uses spectral information is prone to confusion during the identification process. In addition, as a passive remote sensing technology, hyperspectral data acquisition needs to be carried out in a clear and cloudless state, which has high requirements for the acquisition environment. Light detection and ranging (LiDAR) is a remote sensing technology that uses pulsed lasers to measure distance. These pulses are usually combined with information recorded by airborne systems to generate a high-precision three-dimensional model of the target. It is highly adaptable to the acquisition environment and has been widely used in the field of agronomy to indicate the three-dimensional structural characteristics of vegetation in recent years (Hu et al., 2023; Jin et al., 2021; Li et al., 2023; Luo et al., 2021; Rivera et al., 2023; Sandonís-Pozo et al., 2022; Su et al., 2019). At present, the three-dimensional structural characteristics of vegetation under natural gas stress are still unclear, and no research has been conducted to identify micro-leakage points in underground natural gas storage reservoirs based on the three-dimensional structural characteristics of vegetation. Summary of the invention

[0005] In view of this, the present invention provides a natural gas micro-leakage stress vegetation detection method based on LiDAR data, which can effectively detect the micro-leakage points of natural gas.

[0006] The technical solution adopted by the present invention is:

[0007] A method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data, comprising the following steps:

[0008] S100, scanning a target plot of land by a 3D laser scanner to obtain point cloud data of the target plot of land;

[0009] S200, pre-processing the point cloud data to obtain CHM data;

[0010] S300, constructing a gray level co-occurrence matrix based on the CHM data to obtain homogeneity data and contrast data of each pixel;

[0011] S400, assigning homogeneity data and contrast data to CHM data respectively, so that each pixel object x i There are first data of height, second data of homogeneity and third data of contrast respectively;

[0012] S500, randomly selecting k pixels from all pixels, where the k pixels are called first central pixels, and the other pixels are called classification pixels;

[0013] S600, obtaining the Euclidean distance between each classified pixel and each central pixel according to the Euclidean model, and binding each classified pixel and the central pixel with the smallest Euclidean distance into one cluster, thereby obtaining k clusters, namely C 1 , C 2 ……C k ;

[0014] S700, re-determine the second central pixel of each cluster, and use the k second central pixels as the first central pixel, and re-execute step S600 until the second central pixel remains unchanged. The second central pixel at this time is called the third central pixel;

[0015] S900, comparing k third central pixels, selecting the third central pixel with the largest first data, the cluster where the third central pixel belongs is called the standard cluster, and the other clusters are called comparison clusters;

[0016] S1000 , calculating the JM distance between each comparison cluster and the standard cluster, and selecting comparison clusters with a JM distance greater than 1.8.

[0017] Preferably, in step S200, the preprocessing includes filtering noise points, generating DEM data and DSM data based on the point cloud data, and then obtaining CHM data by subtracting the DEM data from the DSM data.

[0018] Preferably, in step S300, homogeneity is expressed as:

[0019]

[0020] The contrast ratio is expressed as:

[0021]

[0022] Among them, P(i,j) represents the gray level co-occurrence matrix, i and j represent the gray level values, which are taken from the gray level range of the image, and their values ​​range from 1 to N, where N is the number of all possible gray levels in the image.

[0023] Preferably, in step S400, all pixels constitute a data set X, X = {x 1 ,x 2 ,…,x n}, each pixel x i is a three-dimensional vector (x i1 、x i2 、xi3 ), where x i1 is the height of the i-th pixel, x i2 is the homogeneity of the i-th pixel, x i3 is the contrast of the i-th pixel;

[0024] In step S600, the Euclidean model is:

[0025]

[0026] Among them, d(x i ,μ j ) represents the classified pixel x i With the first central pixel μ j The Euclidean distance between il and μ jl Classify pixels x separately i and the first central pixel μ j The value in the first dimension.

[0027] Preferably, in step S700, the second central pixel of each cluster is determined according to the central model. j The second central pixel is μ' j ,

[0028]

[0029] Among them, μ' j Cluster C j The second center pixel of |C j | is cluster C j The number of pixels in .

[0030] Preferably, the method further comprises the steps of:

[0031] S800, with k=2, 3, ..., n, respectively executing steps S500-S700, thereby obtaining multiple cluster groups, obtaining the sum of square errors of each cluster group through a sum of square error model, and determining the best cluster group;

[0032] The sum of squared error model is:

[0033]

[0034] Where SSE represents the within-cluster sum of squares error, k is the number of clusters in the cluster group, and x i Belongs to cluster C j Pixels, μ” j Cluster C j The third center pixel of ||x i -μ" j || 2 For xi With μ" j The square of the Euclidean distance between them;

[0035] S900 , comparing k third central pixels in the best cluster group.

[0036] Preferably, in step S800, a function curve of the sum of square errors and k is constructed, and the k value corresponding to the elbow of the function curve is the k value of the optimal cluster group.

[0037] Preferably, in step S1000, the JM distance is

[0038] JMD ij =2(1-e -α )

[0039]

[0040] Among them, α is the Bhattacharyya distance, μ” i is the comparison cluster mean vector, μ” j is the mean vector of the standard cluster, σ i is the corresponding matrix covariance of the contrast cluster, σ j is the covariance of the corresponding matrix of the standard cluster;

[0041] Filter out 1.8<JMD ij These clusters are called suspected stress clusters.

[0042] Preferably, the method further comprises the steps of:

[0043] S1100: among all suspected stress clusters, interconnected suspected stress clusters are merged to form connected candidate regions, circular connected candidate regions are screened out, and the circular connected candidate regions are target regions. The stress center and stress range are determined according to the target regions and the objective function model.

[0044] The objective function model is:

[0045]

[0046] Among them, (x i ,y i ) is the coordinate of the point at the edge of the target area, (x c ,y c ) are the coordinates of the stress center, and r is the radius of the stress area.

[0047] Beneficial effects of the present invention:

[0048] This algorithm is based on vegetation LiDAR data. It uses the Kmeans clustering algorithm to perform cluster analysis on the reduced-dimensional LiDAR data and classifies the clustering results under the constraint of the JM distance criterion. It can monitor the micro-leakage of underground natural gas storage tanks in vegetation areas without prior knowledge. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The above and other objects, features and advantages of the present invention will become more apparent through the following description of the embodiments of the present invention with reference to the accompanying drawings, in which:

[0050] Figure 1 This is a flow chart of a natural gas micro-leakage stress vegetation detection method based on LiDAR data;

[0051] Figure 2 It is a result display diagram;

[0052] Figure 3 It is the function curve of the sum of squared errors and k. DETAILED DESCRIPTION

[0053] The present invention is described below based on embodiments, but the present invention is not limited to these embodiments. In the following detailed description of the present invention, some specific details are described in detail. In order to avoid confusing the essence of the present invention, known methods, processes, procedures, and components are not described in detail.

[0054] In addition, persons of ordinary skill in the art will appreciate that the drawings provided herein are for illustration purposes and are not necessarily drawn to scale.

[0055] Unless the context clearly requires otherwise, throughout the specification and claims, the words "include", "comprising" and similar words should be interpreted in an inclusive sense rather than an exclusive or exhaustive sense; that is, in the sense of "including but not limited to".

[0056] In the description of the present invention, it should be understood that the terms "first", "second", etc. are only used for descriptive purposes and cannot be understood as indicating or implying relative importance. In addition, in the description of the present invention, unless otherwise specified, "plurality" means two or more.

[0057] In the present invention, the direction is defined according to the surgical process, the side away from the operator during the surgery is called "far", and the side close to the operator is called "near".

[0058] See also Figure 1-Figure 3The present invention provides a natural gas micro-leakage coercion vegetation detection method based on LiDAR data, which is used to detect the micro-leakage point of natural gas. When a natural gas micro-leakage occurs, it is not easy to be detected due to the small leakage amount. If it is not discovered in time, the leakage point will be corroded and enlarged, and serious safety accidents will occur. However, natural gas will inhibit the growth of plants around the surface of the leakage point. Therefore, by judging the growth of plants, it is determined whether there is a natural gas leak. When there is a leak, the leak point can be judged quickly, so that the micro-leakage point can be handled in time to avoid the expansion of the leakage accident.

[0059] The method comprises the steps of:

[0060] S100, scanning a target plot of land by a 3D laser scanner to obtain point cloud data of the target plot of land;

[0061] S200, pre-processing the point cloud data to obtain CHM data;

[0062] S300, constructing a gray level co-occurrence matrix based on the CHM data to obtain homogeneity data and contrast data of each pixel;

[0063] S400, assigning homogeneity data and contrast data to CHM data respectively, so that each pixel object x i There are first data of height, second data of homogeneity and third data of contrast respectively;

[0064] S500, randomly selecting k pixels from all pixels, where the k pixels are called first central pixels, and the other pixels are called classification pixels;

[0065] S600, obtaining the Euclidean distance between each classified pixel and each central pixel according to the Euclidean model, and binding each classified pixel and the central pixel with the smallest Euclidean distance into one cluster, thereby obtaining k clusters, namely C 1 , C 2 ……C k ;

[0066] S700, re-determine the second central pixel of each cluster, and use the k second central pixels as the first central pixel, and re-execute step S600 until the second central pixel remains unchanged. The second central pixel at this time is called the third central pixel;

[0067] S900, comparing k third central pixels, selecting the third central pixel with the largest first data, the cluster where the third central pixel belongs is called the standard cluster, and the other clusters are called comparison clusters;

[0068] S1000 , calculating the JM distance between each comparison cluster and the standard cluster, and selecting comparison clusters with a JM distance greater than 1.8.

[0069] For example, a field experiment on simulating the response of grassland to natural gas leakage was conducted in Changziying, Daxing District, Beijing (coordinates are 39°39′2.56″N, 116°34′33.10″E), aiming to study their response to natural gas stress under controlled conditions. On June 9, 2019, the grassland was planted in 8 2.5×2.5 m small plots, of which the control plots and experimental plots were alternately distributed, with 4 plots each. When the vegetation grew 2-3 layers of leaves and the growth was relatively stable, natural gas was continuously injected at a leakage rate of 1 liter / minute at a depth of about 60 cm in the center of each experimental plot through a gas control device and pipeline until the end of the experiment. The GeoSLAM ZEB HORIZON mobile laser scanner (PLS) developed by GeoSLAM Co., Ltd. (Nottingham, UK) was used to collect point cloud data for each plot. The device has a field of view (FOV) of 360° (H) × 270° (V), an angular resolution of 0.2° (H) × 2° (V), and acquires 300,000 scan points per second. The ZEB Horizon has a range of 100m and an overall accuracy of ±3mm.

[0070] Therefore, in step S100, the target plot can be scanned by a laser scanner to obtain point cloud data of the target plot. The point cloud data here are scanned points, and each point is represented by a three-dimensional coordinate, that is, a plane position and a height.

[0071] In addition, ground point classification, DEM (digital ground model) generation and CHM (vegetation canopy height model) generation are also required. The ground point filtering algorithm based on cloth simulation separates the point cloud into ground points and vegetation points. The specific parameters are that the scene is plane, the cloth resolution is 0.025m, the maximum number of iterations is 500, and the classification threshold is 0.1. Based on the ground points and vegetation points, the inverse distance weighted interpolation algorithm is applied to generate DEM and DSM grids with a resolution of 0.025 meters respectively (Munyati and Sinthumule, 2021). Finally, CHM is constructed by subtracting DEM from DSM. This process essentially realizes the normalization of vegetation height, eliminates the influence of terrain elevation, and enables the analysis to focus on the relative height changes of vegetation or objects of interest.

[0072] In step S300, three three-dimensional structural features are extracted through CHM: height, homogeneity and contrast, which is specifically achieved by using gray level co-occurrence matrix.

[0073] Homogeneity is expressed as:

[0074]

[0075] The contrast ratio is expressed as:

[0076]

[0077] Where P(i,j) represents the gray-level co-occurrence matrix, i and j represent the gray-level values, which are taken from the gray-level range of the image, and their values ​​range from 1 to N, where N is the number of all possible gray-levels in the image. i and j represent the gray-level values, which are taken from the gray-level range of the image, and their values ​​range from 1 to N, where N is the number of all possible gray-levels in the image. Specifically, the gray-level co-occurrence matrix is ​​calculated in a sliding window of size 25x25 and in four directions of 0°, 45°, 90° and 135°, and the mean of each feature in each direction is taken as the final feature parameter result.

[0078] At this time, CHM is pixelated, that is, each point (pixel) has coordinates on the plane to express its position on the plane. At the same time, each pixel is also assigned three values: height, homogeneity and contrast.

[0079] Therefore, all points can be expressed by data. For example, the position of any pixel on the plane can be expressed by x. i To express, pixels also have height, homogeneity, and contrast, which can be expressed by vectors, that is, (x i1 、x i2 、x i3 ).

[0080] In this way, all pixel sets are formed into a data set X, X = {x 1 ,x 2 ,…,x n}.

[0081] The purpose of the present invention is to find vegetation under natural gas stress, specifically by clustering.

[0082] In step S500, k pixels are randomly selected as first central pixels, for example, k=5, and the five first central pixels are expressed as μ 1 , μ 2 ,…,μ 5 , all pixels other than the five first center pixels are called classified pixels.

[0083] In step S600, the Euclidean distance between each classified pixel and the five first central pixels is calculated in turn to obtain the first central pixel with the smallest Euclidean distance, and the classified pixel is bound to the first central pixel.

[0084] For example, to classify a pixel x 3 Respectively with μ 1 , μ 2 ,…,μ 5Calculate the Euclidean distance. According to the calculation result, x 3 With μ 2 The Euclidean distance between them is the shortest, so x 3 With μ 2 Binding; similarly, each classified pixel is bound in the above manner, thus forming 5 clusters C1, C2...C5.

[0085] The above Euclidean calculation is implemented through the Euclidean model, which is:

[0086]

[0087] Among them, d(x i ,μ j ) represents the classified pixel x i With the first central pixel μ j The Euclidean distance between il and μ jl Classify pixels x separately i and the first central pixel μ j The value in the first dimension.

[0088] Since the initial first central pixel is randomly generated, it is not the center of each cluster in a real sense, and thus the center of each cluster is re-determined through S700.

[0089] In step S700, the second central pixel of each cluster is determined according to the central model. j The second central pixel is μ' j ,

[0090]

[0091] Among them, μ' j Cluster C j The second center pixel of |C j | is cluster C j The number of pixels in x i Cluster C j Pixels in .

[0092] Since there are k clusters, k second central pixels can be obtained. Then, these k second central pixels are used as first central pixels, and step S600 is executed again to redetermine the k clusters. Then step S700 is executed again to redetermine the second central pixels of the k clusters. This process is repeated until the k second central pixels remain unchanged (third central pixels). In this way, all pixels can be divided into the final k clusters, achieving the purpose of clustering.

[0093] The most appropriate number of clusters to divide the pixels into is the most appropriate value of k.

[0094] The present invention determines the optimal number of clusters by using the elbow method.

[0095] Specifically, with k=2, 3, ...n, steps S500-S700 are respectively performed to obtain multiple cluster groups, for example, one cluster group when k=2, one cluster group when k=3, and so on, and the sum of square errors are calculated for each cluster group.

[0096] The sum of squared error model is:

[0097]

[0098] Where SSE represents the within-cluster sum of squares error, k is the number of clusters in the cluster group, and x i Belongs to cluster C j Pixels, μ” j Cluster C j The third center pixel of ||x i -μ" j || 2 For x i With μ" j The square of the Euclidean distance between them;

[0099] Therefore, the change line of SSE can be proposed according to the change of k from small to large. The k value of the point corresponding to the elbow of the change line is the optimal cluster group, thereby determining the value of k. Figure 3 In the illustrated embodiment, k=5 is the optimal cluster group.

[0100] In step S900, since the K value is determined and the k clusters are also determined, the third central pixel with the largest first data (height) is picked out, and the cluster corresponding to it is the standard cluster, and the other clusters are comparison clusters.

[0101] The standard cluster has the highest height and the corresponding vegetation has the highest height, so it can be defined as a healthy vegetation cluster. By comparing other comparison clusters with the standard cluster, it can be determined whether these comparison clusters are under stress.

[0102] Specifically, in step S1000, whether the standard cluster is under duress is determined by using the JM distance.

[0103] The JM distance is

[0104] JMD ij =2(1-e -α )

[0105]

[0106] Among them, α is the Bhattacharyya distance, μ” i is the comparison cluster mean vector, μ” jis the mean vector of the standard cluster, σ i is the corresponding matrix covariance of the contrast cluster, σ j is the covariance of the corresponding matrix of the standard cluster;

[0107] The JM distance has convergence, and the criteria for judging separability are as follows: when 0 < JMD ij ≤1, the two categories are not separable; when 1<JMD ij ≤1.8, the two categories are somewhat separable, but overlap to a certain extent; when 1.8<JMD ij When ≤2, the two categories are well separable, and 1.8<JMD is selected. ij These clusters are called suspected stress clusters.

[0108] Clusters that are more separable from healthy vegetation clusters, that is, clusters with JM distance values ​​greater than 1.8, indicate that the corresponding vegetation structure characteristics have changed significantly. The pixels belonging to this cluster are classified as "suspected stress pixels", while the pixels belonging to clusters that are less separable from healthy vegetation clusters are classified as "non-stress pixels". Finally, the "suspected stress pixels" are merged to form a candidate area for stress pixel connectivity.

[0109] In step S1100, circular connected candidate areas (the impact of natural gas leakage is close to a circle) are screened out. The screening method can be manual screening, or thresholds are set for the three elements of connected area area, body ratio and rectangularity to eliminate non-circular connected candidate areas that do not meet the requirements.

[0110] The area of ​​the connected region is defined as the number of pixels contained in the connected region. The aspect ratio is defined as the ratio of the length and width of the object's smallest circumscribed rectangle. This parameter reflects the degree to which an object is close to a circle. The calculation formula is as follows:

[0111] P=W / L

[0112] Where W and L are the width and length of the object's minimum enclosing rectangle, respectively. Rectangularity is defined as the ratio of the area of ​​an object to the area of ​​its minimum enclosing rectangle. This parameter reflects the degree of "fullness" of an object to its minimum enclosing rectangle. The calculation formula is as follows:

[0113] R=A 0 / A R

[0114] In the formula, A 0 is the area of ​​the connected region, A R is the area of ​​its minimum circumscribed rectangle. The thresholds are set as follows: the connected area is 0.075 square meters, the P body ratio is 0.55, and the R rectangularity is above 0.50. Areas that meet these requirements are considered as coercive areas, and non-circular connected candidate areas will be eliminated.

[0115] The circular connected candidate area is the target area. The coercion center and coercion range are determined according to the target area and the objective function model.

[0116] The objective function model is:

[0117]

[0118] Among them, (x i ,y i ) is the coordinate of the point at the edge of the target area, (x c ,y c ) is the coordinate of the coercion center, and r is the radius of the coercion area. In this way, the location of the natural gas micro-leakage can be obtained, which is convenient for subsequent treatment and prevents the accident from expanding.

[0119] It should be understood that the above-mentioned embodiments are merely illustrative and not restrictive. Without departing from the basic principles of the present invention, various obvious or equivalent modifications or substitutions that can be made by those skilled in the art to the above-mentioned details will all be included in the scope of the claims of the present invention.

Claims

1. A method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data, characterized in that: Includes steps: S100, scanning a target plot of land by a 3D laser scanner to obtain point cloud data of the target plot of land; S200, pre-processing the point cloud data to obtain CHM data; S300, constructing a gray level co-occurrence matrix based on the CHM data to obtain homogeneity data and contrast data of each pixel; S400, assigning homogeneity data and contrast data to CHM data respectively, so that each pixel object x i There are first data of height, second data of homogeneity and third data of contrast respectively; S500, randomly selecting k pixels from all pixels, where the k pixels are called first central pixels, and the other pixels are called classification pixels; S600, obtaining the Euclidean distance between each classified pixel and each central pixel according to the Euclidean model, and binding each classified pixel and the central pixel with the smallest Euclidean distance into a cluster, thereby obtaining k clusters, namely C1, C2, ..., C k ; S700, re-determine the second central pixel of each cluster, and use the k second central pixels as the first central pixel, and re-execute step S600 until the second central pixel remains unchanged. The second central pixel at this time is called the third central pixel; S900, comparing k third central pixels, selecting the third central pixel with the largest first data, the cluster where the third central pixel belongs is called the standard cluster, and the other clusters are called comparison clusters; S1000 , calculating the JM distance between each comparison cluster and the standard cluster, and selecting comparison clusters with a JM distance greater than 1.

8.

2. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 1, characterized in that: In step S200, the preprocessing includes filtering noise points, generating DEM data and DSM data based on point cloud data, and then obtaining CHM data by subtracting DEM data from DSM data.

3. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 2 is characterized in that: In step S300, homogeneity is expressed as: The contrast ratio is expressed as: Among them, P(i,j) represents the gray level co-occurrence matrix, i and j represent the gray level values, which are taken from the gray level range of the image, and their values ​​range from 1 to N, where N is the number of all possible gray levels in the image.

4. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 1, characterized in that: In step S400, all pixels constitute a data set X, X = {x1, x2, ..., x n }, each pixel x i is a three-dimensional vector (x i1 、x i2 、x i3 ), where x i1 is the height of the i-th pixel, x i2 is the homogeneity of the i-th pixel, x i3 is the contrast of the i-th pixel; In step S600, the Euclidean model is: Among them, d(x i ,μ j ) represents the classified pixel x i With the first central pixel μ j The Euclidean distance between il and μ jl Classify pixels x separately i and the first central pixel μ j The value in the first dimension.

5. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 4 is characterized in that: In step S700, the second central pixel of each cluster is determined according to the central model. j The second central pixel is μ' j , Among them, μ' j Cluster C j The second center pixel of |C j | is cluster C j The number of pixels in .

6. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claims 1-5, characterized in that: Also includes the steps: S800, with k=2, 3, ..., n, respectively executing steps S500-S700, thereby obtaining multiple cluster groups, obtaining the sum of square errors of each cluster group through a sum of square error model, and determining the best cluster group; The sum of squared error model is: Where SSE represents the within-cluster sum of squares error, k is the number of clusters in the cluster group, and x i Belongs to cluster C j Pixels, μ” j Cluster C j The third center pixel of ||x i -μ" j || 2 For x i With μ" j The square of the Euclidean distance between them; S900 , comparing k third central pixels in the best cluster group.

7. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 6, characterized in that: In step S800 , a function curve of the sum of squared errors and k is constructed, and the k value corresponding to the elbow of the function curve is the k value of the optimal cluster group.

8. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claims 1-5, characterized in that: In step S1000, the JM distance is JMD ij =2(1-e -α ) Among them, α is the Bhattacharyya distance, μ” i is the comparison cluster mean vector, μ” j is the mean vector of the standard cluster, σ i is the corresponding matrix covariance of the contrast cluster, σ j is the covariance of the corresponding matrix of the standard cluster; Filter out 1.8<JMD ij These clusters are called suspected stress clusters.

9. The method for detecting vegetation stress caused by natural gas micro-leakage based on LiDAR data according to claim 8, characterized in that: Also includes the steps: S1100: among all suspected stress clusters, interconnected suspected stress clusters are merged to form connected candidate regions, circular connected candidate regions are screened out, and the circular connected candidate regions are target regions. The stress center and stress range are determined according to the target regions and the objective function model. The objective function model is: Among them, (x i ,y i ) is the coordinate of the point at the edge of the target area, (x c ,y c ) are the coordinates of the stress center, and r is the radius of the stress area.

Images