Lidar cloud detection method based on bidirectional reconstruction of backscattered signals

Abstract

The present invention discloses a laser radar cloud detection algorithm based on bidirectional reconstruction backscattering signal, including: preprocessing the laser radar backscattering signal; selecting corresponding forward interpolation points according to the backscattering signal characteristics under ideal clean atmospheric conditions, obtaining forward reconstructed signals by Akima interpolation method, combining adaptive threshold and empirical threshold to perform cloud stratification, eliminate misjudgment and determine cloud interval; on this basis, further selecting backward interpolation points and obtaining backward reconstructed signals by Akima interpolation, and obtaining reconstructed cloud-free signals by averaging with forward reconstructed signals; dividing the interval into two sections with the cloud peak as the boundary, selecting several elements in the preprocessed signal with an intensity less than the simulated signal in each section, and determining the cloud base or cloud top height according to the number of elements. The present invention divides the cloud layer more accurately, and achieves a higher cloud layer detection rate and inversion accuracy.

Description

Technical Field

[0001] The present invention relates to a laser radar cloud detection method, and in particular to a laser radar cloud detection method based on bidirectional reconstruction of backscattered signals. Background Art

[0002] Clouds are polymers floating in the atmosphere that are generated by physical reactions of water vapor in the atmosphere. They cover about two-thirds of the Earth's surface and are the tangible result of the Earth's huge water cycle. Due to their physical properties, clouds have become the most important regulator for maintaining energy balance in the Earth-atmosphere system, and are also one of the biggest uncertainties affecting climate change. Therefore, cloud observation plays an important role in local and global climate research and human production and life, and has strategic value.

[0003] As a ground-based active remote sensing device with high temporal and spatial resolution, lidar is mainly used to determine information such as cloud height. In recent years, the algorithm of automatically inverting cloud parameters using lidar has been closely following the development of atmospheric science to meet the needs of scientific application research. At present, cloud detection algorithms based on single-channel backscatter signals can be roughly divided into differential methods, statistical algorithms, and sliding window methods. However, when such methods extract clouds from signals, they require empirical thresholds to directly judge clouds, which are greatly affected by system parameters and subjectivity. In addition, the inversion of cloud top height is prone to the situation where effective cloud tops rather than real cloud tops appear.

[0004] The existence of the above problems seriously restricts the accuracy of lidar measurement of the vertical structure of clouds. Summary of the invention

[0005] Purpose of the invention: The purpose of the present invention is to provide a lidar cloud detection method based on bidirectionally reconstructed backscattered signals, which can achieve accurate inversion of the vertical structure of clouds.

[0006] Technical solution: The laser radar cloud detection method of the present invention comprises the following steps:

[0007] S1, preprocess the backscatter signal of the laser radar to obtain a new backscatter signal B LN ;

[0008] S2, starting from the lowest detection distance, screens out the pre-processed backscatter signal B LN The monotonically decreasing elements in form the sequence B 1 ; According to sequence B 1 The heights corresponding to the elements in are divided into several feature segments T(k) according to the interval continuity, where k = 1, 2, 3..., and the left endpoint of the feature segment T(k) is defined as the left feature point The right endpoint is defined as the right feature point Sequence B 1 As the forward interpolation point, Akima interpolation is performed to obtain the forward reconstructed signal B A1 ;

[0009] S3, the new backscattered signal B LN With the forward reconstructed signal B A1 Subtract and get the residual signal R 1 ; Using adaptive threshold T h Determine the number of cloud layers and use the left and right feature points in step S2 to determine the cloud peaks of each cloud layer interval, where j = 1, 2, 3, ...;

[0010] S4, using empirical thresholds to eliminate misjudgments and merge clouds;

[0011] S5, at Yunfeng The cloud interval is determined based on the interval, and the left and right endpoints of the cloud interval are recorded as and

[0012] S6, starting from the highest detection distance, filter out the pre-processed backscatter signal B LN The monotonically increasing elements in the sequence B 2 ; Sequence B 2 As the backward interpolation point, Akima interpolation is performed to obtain the backward reconstructed signal B A2 ;

[0013] S7, forward reconstructed signal B A1 and the backward reconstructed signal B A2 Average, get the reconstructed simulated cloud-free signal B A ; Compare with cloud-free signal B A With the new backscattered signal B LN The signal strength at each location within the cloud layer is used to determine the cloud base and cloud top heights.

[0014] Furthermore, in step S1, the expression for preprocessing the backscattered signal is as follows:

[0015]

[0016] Where P(r) represents the intensity of the backscattered signal of the lidar at a height of r, O l (r) is the geometric overlap factor correction caused by field of view conflict, C is the system constant of the lidar, E is the emission energy, β(r) is the backscattering coefficient formed by all types of atmospheric scattering, σ(r') is the atmospheric extinction coefficient, It represents the atmospheric transmittance, P b is the background noise, P ap(r) represents the post-pulse noise, D tc (P(r)) is the dead time correction for the detector due to the temporary saturation of photons; Defined as the new backscatter signal B LN .

[0017] Further, in step S2, the new backscattered signal B is first LN The signal strengths are reorganized in ascending order, and the new backscattered signals B within the first k resolution distances at the highest detection distance are LN The minimum intensity value of is used as the standard, and elements smaller than the standard are removed; then the new backscattered signal B is screened out from the lowest detection distance. LN The elements that satisfy the condition that "with increasing height, the backscatter signal strength gradually decreases" form a monotonically decreasing sequence B 1 As the forward interpolation point, and with the spatial resolution as the interval, Akima interpolation is used to obtain the same value as B LN The forward reconstructed signal B corresponding to each position in A1 ;

[0018] The heights corresponding to the elements in each feature segment differ by a distance of spatial resolution; when the length of the feature segment is 1, the element serves as both a left feature point and a right feature point.

[0019] Further, in step S3, when the number of cloud peaks is determined by using the adaptive threshold, the backscatter signal B LN With the forward reconstructed signal B A1 The residual signal R 1 , the expression is as follows:

[0020] R 1 =B LN -B A1

[0021] The residual signal R 1 The mean of the elements greater than zero in is denoted as m 1 , take m 1 A times of is used as the adaptive threshold T for cloud peak detection h , then:

[0022] T h =a*m 1 .

[0023] Further, in step S4, using the distance threshold H 1 To judge a single-layer multi-peak cloud, assume that the xth candidate cloud peak is at and The x+1th candidate cloud peak is between and Between (0≤x≤j-1, 0≤m,n≤k-1), and They are and The corresponding height (km), then the distance threshold H 1 The expression is:

[0024]

[0025] Using the height threshold H 2 Eliminate misjudgment of each cloud peak interval, and the height threshold H 2 The expression is as follows:

[0026]

[0027] Further, in step S5, each cloud peak The cloud interval is the previous characteristic segment With the next feature segment between.

[0028] Furthermore, in step S6, when filtering the backscattered interpolation points, the elements within the cloud interval need to be excluded first; then the new backscattered signal B is filtered out from the highest detection distance. LN The elements that satisfy the condition that "with decreasing height, the backscatter signal strength gradually increases" form a monotonically increasing sequence B 2 As the forward interpolation point, and with the spatial resolution as the interval, Akima interpolation is used to obtain the new backscatter signal B LN The backward reconstructed signal B corresponding to each position in A2 .

[0029] Furthermore, in step S7, the reconstructed cloud-free signal B A The expression is as follows:

[0030]

[0031] Further, in step S7, the candidate cloud intervals are divided into C P Divided into two parts, according to the strong scattering characteristics of clouds to light, in C L -C P Invert the cloud base height within the segment and convert B LN Medium less than B A The element is denoted as S LP (i), the number of elements is denoted by N LP , the corresponding height is recorded as Cloud base height The expression is as follows:

[0032]

[0033] in, Indicates that in C L The cloud base height at

[0034] Similarly, C P -C R Segment B LN Medium less than B A The element is denoted as S PR (i), the number of elements is denoted by N PR , the corresponding height is recorded as Cloud top height The expression is as follows:

[0035]

[0036] Among them, S LP (N LP ) is C L -C P The position closest to the cloud peak in the segment, S LP (1) is C P -C R The location closest to Yunfeng in the section.

[0037] Compared with the prior art, the present invention has the following significant effects:

[0038] 1. The cloud-free signal is simulated by the interpolation and reconstruction method. Based on the empirical threshold, an adaptive threshold is added to more accurately identify and divide the cloud layer.

[0039] 2. Within the candidate cloud interval, the interval is divided into two sections with the cloud peak as the boundary. The intensity of the reconstructed simulated cloud-free signal is used as the reference value. The cloud base height and cloud top height are determined by comparing the preprocessed signal intensity with the corresponding reference value. This solves the problem that when the traditional algorithm uses the cloud base backscatter signal intensity to determine the cloud top height, the inversion result is too low than the true value, and can achieve a higher cloud detection rate and inversion accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 It is the overall flow chart of the present invention;

[0041] Figure 2(a) shows the distance-corrected backscattered signal B used by the traditional algorithm. NR Schematic diagram;

[0042] FIG. 2( b ) is a backscattered signal B after preprocessing in the present invention. LN Schematic diagram;

[0043] FIG. 3( a ) is a sequence B composed of forward interpolation points selected in the present invention. 1and the excluded element D 1 Schematic diagram;

[0044] FIG3( b ) shows the forward reconstructed signal B obtained by Akima interpolation using the selected forward interpolation points in the present invention. A1 Schematic diagram;

[0045] FIG. 4( a ) shows the residual signal R used in the present invention. 1 and adaptive threshold T h Schematic diagram of screening cloud peaks;

[0046] FIG4( b ) is a schematic diagram of determining the cloud peak range by using left and right feature points and combining the empirical threshold to merge the cloud peak intervals and eliminate misjudgments in the present invention;

[0047] FIG. 5( a ) A sequence B composed of backward interpolation points selected in the present invention 2 and the excluded element D 2 Schematic diagram;

[0048] FIG5(b) shows the backward reconstructed signal B obtained by Akima interpolation using the selected backward interpolation points in the present invention. A2 Schematic diagram;

[0049] Figure 6(a) is the cloud-free signal B reconstructed in the present invention A , and compare the pre-processed signal strength to determine the cloud base height and cloud top height of each layer of cloud;

[0050] FIG6( b ) is a schematic diagram of the cloud base height and cloud top height of the first layer of clouds inverted using laser radar measured data in the present invention;

[0051] FIG6( c ) is a schematic diagram of the cloud base height and cloud top height of the second layer of clouds inverted using laser radar measured data in the present invention;

[0052] FIG6( d ) is a schematic diagram of the cloud base height and cloud top height of the third layer of clouds inverted using laser radar measured data in the present invention;

[0053] Figure 7(a) is a schematic diagram of the inversion results of one day's data using the improved differential zero-crossing algorithm;

[0054] Figure 7(b) is a schematic diagram of the inversion results of the value distribution equalization algorithm for one day of data;

[0055] Figure 7(c) is a schematic diagram of the inversion results of one day's data using the improved simple multi-scale algorithm;

[0056] FIG7( d ) is a schematic diagram of the inversion result of one day's data using the algorithm proposed in the present invention;

[0057] Figure 8(a) is a schematic diagram of the correlation coefficient and root mean square error between the cloud base height inverted by the differential zero crossing algorithm and the station results, with the station results as the standard;

[0058] Figure 8(b) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud base height inverted by the value distribution equalization algorithm and the station results, with the station results as the standard;

[0059] Figure 8(c) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud base height inverted by the improved simple multi-scale algorithm and the station results, with the station results as the standard;

[0060] FIG8( d ) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud base height inverted by the algorithm proposed in the present invention and the station results, taking the station results as the standard;

[0061] Figure 9(a) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud top height inverted by the differential zero crossing algorithm and the station results, with the station results as the standard;

[0062] Figure 9(b) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud base height inverted by the value distribution equalization algorithm and the station results, with the station results as the standard;

[0063] Figure 9(c) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud top height inverted by the improved simple multi-scale algorithm and the station results, with the station results as the standard;

[0064] FIG9( d ) is a schematic diagram showing the correlation coefficient and root mean square error between the cloud top height inverted by the algorithm proposed in the present invention and the station results, taking the station results as the standard;

[0065] Fig.10 This is a schematic diagram showing the error comparison between the inversion results of the algorithm proposed in the present invention and other traditional algorithms within one year. DETAILED DESCRIPTION

[0066] The present invention is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0067] In combination with the characteristics of the measured laser radar backscatter signal, this paper proposes a cloud detection algorithm based on bidirectional reconstruction of the backscatter signal, which selects B from the two directions of the lowest and highest detection distances respectively. LN Perform Akima interpolation on the relevant interpolation points in B LN The forward and backward reconstructed signals corresponding to the detection positions are averaged to obtain the simulated signal in the cloudless condition, and its intensity is used as B LN The algorithm first reconstructs the forward signal and B LNThe difference between clouds and aerosols is distinguished by combining adaptive thresholds and empirical thresholds to determine the candidate cloud interval. The candidate cloud interval is then divided into two segments based on the cloud peak to screen B. LN Elements that are smaller than the corresponding reference values ​​are selected, and finally the vertical structure of the cloud is inverted according to the number of screened elements.

[0068] like Figure 1 As shown in FIG. 1 , it is a flow chart of the cloud layer detection method of the present invention, and the specific implementation process includes the following steps:

[0069] Step 1: preprocess the backscatter signal of the laser radar using a logarithmic function to obtain a new backscatter signal B LN .

[0070] The micro-pulse lidar equation can be expressed as:

[0071]

[0072] The present invention moves the geometric overlap correction factor, the dead time correction factor, and the post-pulse noise to the left term of formula (1), and uses a logarithmic function for processing without performing distance correction. The preprocessing of the backscattered signal includes geometric overlap factor correction, dead time factor correction, post-pulse noise correction, and logarithmic function processing, and its expression is as follows:

[0073]

[0074] Where P(r) represents the intensity of the backscattered signal of the lidar at a height of r, O l (r) is the geometric overlap factor correction caused by field of view conflict, C is the system constant of the lidar, E is the emission energy, β(r) is the backscattering coefficient formed by all types of atmospheric scattering, σ(r') is the atmospheric extinction coefficient, It represents the atmospheric transmittance, P b is the background noise, P ap (r) represents the post-pulse noise, D tc (P(r)) is the dead time correction for the detector due to the temporary photon saturation phenomenon. The right side of equation (2) is defined as the new backscatter signal B LN .

[0075] Figure 2(a) and Figure 2(b) show the comparison results of the distance-corrected backscatter signal and the pre-processed backscatter signal. Figure 2(a) is the normalized relative backscatter signal B after distance correction. NR (r), Figure 2(b) shows the preprocessed signal B LN From the figure, we can see that the signals before and after preprocessing have obvious cloud information at 3-4km, 7-8km and 10-11km, but the preprocessed signal BLN The trend of the backscattered signal decreasing with increasing altitude is maintained, effectively overcoming the influence of far-end noise on the detection results.

[0076] Step 2: Starting from the lowest detection distance, filter out the monotonically decreasing element composition sequence B in the pre-processed backscattered signal 1 , B 1 The selection of needs to consider the influence of noise, ambient light and other factors in actual measurement, that is, B LN The minimum value of B often does not appear at the maximum height. LN The signals are reassembled in ascending order to construct a new sequence P S (i), B LN The length of the signal is recorded as n, and the i-th element in the new sequence is in B LN The position in is denoted as I(i), and its expression is as follows:

[0077] P S (i) = B LN (I(i)), i=1,2,3...n, (3)

[0078] The k resolution ranges before the highest detection distance are B LN The minimum strength is denoted as M s , remove B LN Medium s Small elements (excluding the highest detection distance), the excluded element position is recorded as D 1 (i). From formula (3), we can know that D 1 (i) should contain the values ​​of the first several consecutive positions in I(i), namely:

[0079] D 1 (i)=I(i)(ifP s (i)<M s ,i=1,2,3...n) (4)

[0080] After removing the new backscatter signal B LN Middle D 1 (i) After the elements at position B, start from the lowest detection distance and select LN The elements that satisfy the condition that "with increasing height, the backscatter signal intensity gradually decreases" form a monotonically decreasing sequence B 1 As the forward interpolation point, and use Akima interpolation to get the same value as B at the interval of spatial resolution LN The forward reconstructed signal B corresponding to each position in A1 As shown in FIG3(a), in the present invention, k is set to 10. 1 Perform Akima interpolation as the forward interpolation point to obtain the forward reconstructed signal BA1 , as shown in Figure 3(b).

[0081] Sequence B 1 The heights corresponding to the elements in are divided into several feature segments T(k) (k=1,2,3...) according to their continuity. That is, the heights corresponding to the elements in each feature segment differ by a distance of spatial resolution (0.03km). The left endpoint of the feature segment is defined as the left feature point The right endpoint is defined as the right feature point If the length of the feature segment is 1, the element is used as both the left feature point and the right feature point.

[0082] Step 3: Convert the backscattered signal B LN With the forward reconstructed signal B A1 The residual signal R is obtained by subtracting 1 . Using the adaptive threshold T h Determine the number of cloud layers and use the left and right feature points in step 2 to determine the cloud peaks of each cloud layer interval.

[0083] When using adaptive threshold to determine the number of cloud peaks, the preprocessed signal B is used. LN With the forward reconstructed signal B A1 The residual R 1 , which is expressed as follows:

[0084] R 1 =B LN -B A1 (5)

[0085] R 1 The elements greater than zero in are the ones obtained by Akima interpolation and B LN The residual at the corresponding position in the 1 , take m 1 A times of is used as the adaptive threshold T for cloud peak detection h In the present invention, a is 3, and the expression is as follows:

[0086] T h =a*m 1 (6)

[0087] Therefore, R 1 Greater than the threshold T h Each maximum point of is a potential cloud peak, defined as As shown in Figure 4(a), each cloud peak appears between two characteristic segments, that is, each cloud peak In a previous feature segment With the next feature segment between.

[0088] Step 4: Use the empirical threshold to eliminate misjudgments and merge clouds.

[0089] It is generally believed that the distance between two layers of clouds is at least greater than 0.05 km. In order to avoid the situation where a single-layer multi-peak cloud is mistakenly classified as a multi-layer cloud, it is also necessary to judge the distance of adjacent candidate cloud peaks based on the spatial resolution of the lidar and classify the single-layer multi-peak cloud as a single layer of cloud. and The x+1th candidate cloud peak is between and Between (0≤x≤j-1, 0≤m,n≤k-1). and They are and If a single-layer multi-peak cloud occurs, the cloud peak intervals are merged and the cloud peak positions are updated. The data spatial resolution used in this invention is 0.03 km, and the distance threshold H is used. 1 The expression for merging a single-layer multi-peak cloud is as follows:

[0090]

[0091] Although the residual signal and adaptive threshold can detect most of the cloud peaks, the residual signal R within about 1 to 2 km near the detection end is 1 It is easily affected by thick aerosols, and strong backscattering signals similar to those of clouds appear, resulting in misjudgment. The residual signal R at the far end of the detection is about greater than 17 km. 1 , it may be mistaken for a cloud peak due to the low signal-to-noise ratio. The difference from the real cloud layer is that in these two cases, the left feature point T L With the right feature point T R The distance between them is generally very small. Therefore, the present invention uses the thickness threshold H 2 Eliminate misjudgments of each candidate cloud peak interval. For the xth candidate cloud peak in and Between 2 The expression is as follows:

[0092]

[0093] The final determined cloud peak and cloud peak interval are shown in Figure 4(b).

[0094] Step 5: Determine the cloud layer interval based on the cloud peak interval. The left and right endpoints of the cloud layer interval are recorded as and

[0095] Since the elements in each feature segment also meet the monotonically decreasing condition, they meet the characteristics of the cloud peak-cloud top segment. Therefore, based on the cloud peak interval determined in step 2 and step 3, the cloud layer interval is expanded to the right feature point of the previous feature segment and the right feature point of the next feature segment, that is, each cloud peak The cloud interval is the previous characteristic segment With the next feature segment between, recorded as the left endpoint of the cloud interval With right endpoint

[0096] Step 6: Starting from the highest detection distance, filter out the monotonically increasing element composition sequence B in the pre-processed backscatter signal 2 , use it as the backward interpolation point to perform Akima interpolation to obtain the backward reconstructed signal B A2 .

[0097] The set of elements in each cloud interval is denoted as D 2 (i) After removing B LN Middle D 2 (i) After the elements at position B, start from the highest detection distance and select LN The elements that satisfy the condition that "with decreasing height, the backscatter signal strength gradually increases" form a monotonically increasing sequence B 2 As the backward interpolation point, as shown in Figure 5(a). 2 The Akima interpolation is used to obtain the B LN The backward reconstructed signal B corresponding to each position in A2 , as shown in Figure 5(b).

[0098] Step 7: B A1 With B A2 After averaging, we get the reconstructed cloud-free signal B A , whose intensity is B LN Then the candidate cloud layer interval is divided into two sections based on the cloud peak, and the B LN Elements that are smaller than the corresponding reference values ​​are selected, and finally the cloud base and cloud top heights are determined based on the number of screened elements.

[0099] The forward reconstructed signal B A1 and the backward reconstructed signal B A2 The reconstructed signal B is obtained by averaging A , as shown in Figure 6(a), its expression is as follows:

[0100]

[0101] Reconstructed signal B A Simulated BLN The signal under cloud-free conditions will reconstruct signal B A The strength of B LN For the inversion of cloud base and cloud top height, the present invention divides the candidate cloud layer interval into C P Divide into two, respectively in C L -C P The cloud base height is inverted within the segment, in C P -C R Invert the cloud top height within the segment. LN The determination of the mid-cloud base and cloud top positions is shown in Figure 7(b), (c), and (d), respectively.

[0102] Since the extinction coefficient of clouds is larger than that of aerosols, C L -C P Segment B LN The overall trend is to rise rapidly to the cloud peak, while C P -C R Segment B LN The overall trend is that the cloud peak position drops rapidly until it reaches the background noise level. L -C P B LN Medium less than B A The element is denoted as S LP (i), the number of elements is denoted by N LP , the corresponding height is recorded as Cloud base height The expression is as follows:

[0103]

[0104] in, Indicates that in C L The cloud base height at .

[0105] Similarly, C P -C R Segment B LN Medium less than B A The element is denoted as S PR (i), the number of elements is denoted by N PR , the corresponding height is recorded as Cloud top height The expression is as follows:

[0106]

[0107] In formula (10), S LP (N LP ) is C L -C PThe closest position to Yunfeng in the segment is B LN Medium less than B A And the distance from C P The nearest position is defined as the cloud base C B In formula (11), S LP (1) is C P -C R The closest position to Yunfeng in the segment is B LN Medium less than B A And the distance C P The nearest location is defined as Genting C T If C L -C P or C P -C R Segment does not meet B LN Less than B A This condition is the element of C L As C B or C R As C T .

[0108] To verify the applicability of the method of the present invention, the following verification is performed:

[0109] Verification Example 1

[0110] The 200th profile in the VAP of the Southern Great Plains (SGP) station of the U.S. Atmospheric Radiation Measurement (ARM) on April 22, 2020 was selected for inversion. Figure 6(a) shows the reconstructed simulated cloud-free signal and the preprocessed signal; Figure 6(b) shows the inversion result of the first layer of clouds, with the cloud base and cloud top heights of 3.35km and 4.04km respectively; Figure 6(c) shows the inversion result of the second layer of clouds, with the cloud base and cloud top heights of 7.10km and 8.00km respectively; Figure 6(d) shows the inversion result of the third layer of clouds, with the cloud base and cloud top heights of 9.89km and 10.94km respectively. When inverting the cloud base height of this profile, the error was controlled within 3 resolution distances, and the error of the cloud top height was controlled within 7 resolution distances.

[0111] Verification Example 2

[0112] A total of 24 hours and 2880 profiles were selected from the SGP station from 0:00 on January 3, 2020 to 24:00 on January 3, and the improved differential zero crossing method (IDZC), the value distribution equilibrium method (VDE), the improved simple multi-scale method (SMS-SF) and the algorithm of the present invention (BRBS) were used to invert the vertical structure of the cloud layer, as shown in the following figure. Figure 7(a) , 7(b) , 7(c) and 7(d), it can be seen that BRBS intuitively has the most satisfactory inversion results.

[0113] In order to qualitatively analyze the data quality of cloud information inverted by different algorithms, the cloud information provided by the SGP station is used as the standard, and the inversion results of each algorithm are matched with the station results. The Pearson correlation coefficient (R) and root mean square error (ERMS) are used to characterize the performance of the four algorithms. For the same detected cloud layer, the inversion performance of the cloud base height is as follows: Figure 8(a) , 8(b) , 8(c), and 8(d), the inversion performance of the cloud top height is as follows: Figure 9(a) , 9(b) , 9(c), 9(d). In Figure 8, the cloud base heights inverted by IDZC and VDE are generally higher than the actual situation, while the cloud base heights inverted by SMS-SF and BRBS are generally consistent. In comparison, the correlation coefficient and root mean square error of the cloud base of BRBS indicate that the algorithm proposed in this invention can achieve the requirement of being closer to the actual situation. Figure 9(a) , 9(b) The four algorithms in Figures 9(c) and 9(d) all show a phenomenon that they are lower than the actual situation overall. This is because the laser emitted by the lidar attenuates to the background noise level before it passes through the cloud layer, which is also an unavoidable situation in the inversion. However, compared with other algorithms, the inversion results of BRBS are the most concentrated compared with the actual situation, with a correlation coefficient of more than 0.9 and a root mean square error of less than 300m.

[0114] Verification Example 3

[0115] The cloud VAPs that passed the quality assurance within one year at the SGP station in 2020 were selected for error analysis. For each location indicated by the time resolution and height resolution, if both the algorithm inversion result and the station result indicate that there are clouds at that location, it is considered a detection; if the algorithm inversion result indicates that there are clouds at that location and the station result shows that there are no clouds, it is considered a false detection; if the algorithm inversion result indicates that there are no clouds at that location and the station result indicates that there are clouds, it is considered a missed detection. The detection results of the IDZC, VDE, SMS-SF, and BRBS algorithms were statistically analyzed to obtain the detection rate (Detection Rate, R) of each algorithm for clouds in all profiles at the SGP station in 2020. D ), Error detection Rate, R E ) and Miss Rate, R M ),like Fig.10As shown in the figure, IDZC has a detection rate of nearly 60%, but it is affected by too many zero crossings of the actual first-order differential signal, and has a high misjudgment rate; VDE has the lowest detection rate and misjudgment rate, and the highest missed detection rate, which is consistent with the situation in Figure 7(b) where the detected cloud thickness is smaller than the actual one. SMS-SF has a lower misjudgment rate than BRBS, because BRBS can detect more thick aerosols within 0-3km that are mistakenly classified as clouds, and because BRBS has gotten rid of the problems of traditional algorithms in defining cloud tops, it has the lowest missed detection rate in the one-year experiment. Comprehensive analysis shows that BRBS can achieve the highest detection rate, the lowest missed detection rate, and the misdetection rate is also controlled at around 3%, with the best inversion effect.

[0116] The embodiments of the present invention are described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A lidar cloud detection method based on bidirectional reconstruction of backscattered signals, characterized in that: The following steps are involved: S1, preprocess the backscatter signal of the laser radar to obtain a new backscatter signal B LN ; S2, starting from the lowest detection distance, screens out the pre-processed backscatter signal B LN The monotonically decreasing elements in the sequence B1 form a sequence B1; according to the heights corresponding to the elements in the sequence B1, they are divided into several feature segments T(k) according to the interval continuity, where k = 1, 2, 3..., and the left endpoint of the feature segment T(k) is defined as the left feature point The right endpoint is defined as the right feature point Take sequence B1 as the forward interpolation point to perform Akima interpolation and obtain the forward reconstructed signal B A1 ; First, the new backscattered signal B LN The signal strengths are reorganized in ascending order, and the new backscattered signals B within the first k resolution distances at the highest detection distance are LN The minimum intensity value of is used as the standard, and the elements smaller than this standard are removed; Then, the new backscattered signal B is screened out from the lowest detection distance. LN The elements that satisfy the "backscatter signal intensity gradually decreases with increasing height" form a monotonically decreasing sequence B1 as the forward interpolation point, and use Akima interpolation to obtain the same value as B with the spatial resolution as the interval. LN The forward reconstructed signal B corresponding to each position in A1 ; The heights of the elements in each feature segment differ by a distance of spatial resolution. When the length of the feature segment is 1, the element is used as both a left feature point and a right feature point. S3, the new backscattered signal B LN With the forward reconstructed signal B A1 The residual signal R1 is obtained by performing the difference; using the adaptive threshold T h Determine the number of cloud layers and use the left and right feature points in step S2 to determine the cloud peaks of each cloud layer interval, where j = 1, 2, 3, ...; When using adaptive threshold to determine the number of cloud peaks, the backscatter signal B LN With the forward reconstructed signal B A1 The residual signal R1 is expressed as follows: R1=B LN -B A1 The mean of the elements greater than zero in the residual signal R1 is recorded as m1, and a times of m1 is taken as the adaptive threshold T for cloud peak detection. h , then: T h =a*m1; S4, using empirical thresholds to eliminate misjudgments and merge clouds; S5, at Yunfeng The cloud interval is determined based on the interval, and the left and right endpoints of the cloud interval are recorded as and S6, starting from the highest detection distance, filter out the pre-processed backscatter signal B LN The monotonically increasing elements in the sequence B2 form the sequence B2; the sequence B2 is used as the backward interpolation point for Akima interpolation to obtain the backward reconstructed signal B A2 ; S7, forward reconstructed signal B A1 and the backward reconstructed signal B A2 The reconstructed simulated cloud-free signal B is obtained by averaging A ; Compare with cloud-free signal B A With the new backscattered signal B LN The signal strength at each location within the cloud layer is used to determine the cloud base and cloud top heights.

2. The method for detecting cloud layer by laser radar based on bidirectional reconstruction of backscattered signals according to claim 1, characterized in that: In step S1, the expression for preprocessing the backscattered signal is as follows: Where P(r) represents the intensity of the backscattered signal of the lidar at a height of r, O l (r) is the geometric overlap factor correction caused by the field of view conflict, C is the system constant of the lidar, E is the emission energy, β(r) is the backscattering coefficient formed by all types of atmospheric scattering, σ(r') is the atmospheric extinction coefficient, It represents the atmospheric transmittance, P b is the background noise, P ap (r) represents the post-pulse noise, D tc (P(r)) is the dead time correction for the detector due to the temporary saturation of photons; Defined as the new backscatter signal B LN .

3. The method for detecting cloud layer by laser radar based on bidirectional reconstruction of backscattered signals according to claim 1, characterized in that: In step S5, each cloud peak The cloud interval is the previous characteristic segment With the next feature segment between.

4. The method for detecting cloud layer by laser radar based on bidirectional reconstruction of backscattered signals according to claim 1, characterized in that: In step S6, when screening the backscattered interpolation points, it is first necessary to exclude the elements within the cloud interval; then, the new backscattered signal B is screened from the highest detection distance. LN The elements that satisfy the "backscatter signal intensity gradually increases with decreasing height" form a monotonically increasing sequence B2 as the forward interpolation point, and use Akima interpolation to obtain the new backscatter signal B with the spatial resolution as the interval. LN The backward reconstructed signal B corresponding to each position in A2 .

5. The method for detecting cloud layer by laser radar based on bidirectional reconstruction of backscattered signals according to claim 4 is characterized in that: In step S7, the reconstructed cloud-free signal B A The expression is as follows:

6. The method for detecting cloud layer by laser radar based on bidirectional reconstruction of backscattered signals according to claim 1, characterized in that: In step S7, the candidate cloud intervals are divided into P Divided into two parts, according to the strong scattering characteristics of clouds to light, in C L -C P Invert the cloud base height within the segment and convert B LN Medium less than B A The element is denoted as S LP (i), the number of elements is denoted by N LP , the corresponding height is recorded as Cloud base height The expression is as follows: in, Indicates that in C L The cloud base height at Similarly, C P -C R Segment B LN Medium less than B A The element is denoted as S PR (i), the number of elements is denoted by N PR , the corresponding height is recorded as Cloud top height The expression is as follows: Among them, S LP (N LP ) is C L -C P The position closest to the cloud peak in the segment, S LP (1) is C P -C R The location closest to Yunfeng in the section.

Images