Iterative inversion method for measuring visibility by laser radar with abnormal value self-correction function

Abstract

The invention discloses an abnormal value self-correction laser radar measurement visibility iterative inversion method, and relates to the technical field of meteorological detection equipment, and the method comprises the steps: S1, carrying out the measurement of atmospheric optical characteristics through a laser radar, and obtaining the measurement data of the laser radar; s2, performing dark count correction on the group of data; s3, partitioning the group of data by using an innovative method, automatically eliminating abnormal values, and adaptively filling proper echo signal values; s4, utilizing a slope method to obtain an extinction coefficient; s5, substituting the extinction coefficient obtained by the slope method in S4 into a Klett formula for iteration, and stopping iteration until a difference value between a certain result and a corresponding iteration initial value is smaller than a set threshold value; and S6, calculating the atmospheric visibility according to the relation between the extinction coefficient and the visibility. The method has the advantages that the phenomenon that the visibility is abnormally reduced due to abnormal data caused by smoke dust, agglomerate fog or hard objects on a detection path can be avoided, and the like.

Description

Technical Field

[0001] The present invention relates to the technical field of meteorological detection equipment, and particularly relates to an iterative inversion method for lidar-measured visibility with outlier self-correction. Background Art

[0002] Visibility is the maximum horizontal distance at which a normal human eye can identify a black target from the sky as the background under specific weather conditions at a certain moment. Atmospheric visibility is an important indicator characterizing atmospheric transparency and has important significance in fields such as human health, transportation, aviation, and navigation.

[0003] When using the slope method to invert lidar data, due to the inhomogeneity of the atmosphere, it cannot comprehensively reflect atmospheric information; the Klett iterative method can solve this problem. However, when there are dust, hard targets, cluster fog, etc. on the detection path, the echo signal will increase abnormally. At this time, due to the characteristic of the Klett iterative method that it can comprehensively reflect the atmospheric situation, this outlier will cause the visibility to decrease abnormally, resulting in an incorrect measurement result.

[0004] The outlier processing methods known to those skilled in the art have a risk of failure in lidar power signal processing. The risk stems from the data characteristics of lidar signals: non-linear attenuation of echo signals, interference from multiple scattering, and noise superposition caused by atmospheric aerosol turbulence, etc. In addition, the existing visibility inversion technologies related to outlier identification generally have two major defects: First, the optimal linear region screening mechanism is adopted, resulting in the abandonment of valid data in other regions, seriously reducing the utilization rate of detection data, leading to inaccurate inversion results and being unable to fully reflect the true atmospheric situation outside the region; Second, the visibility inversion based on the slope method can only obtain an approximate solution in the local linear interval, while the visibility inversion by the Klett method will cause the visibility to decrease abnormally due to outliers in the region, and the visibility inversion result for the selected linear region itself is also not accurate enough and has low robustness. Summary of the Invention

[0005] To solve the above technical problems, the present invention provides an iterative inversion method for lidar-measured visibility with outlier self-correction. On the basis of retaining the comprehensive reflection of the atmosphere by the original Klett algorithm, an outlier dynamic discrimination mechanism is creatively introduced: by calculating the median absolute deviation of the extinction coefficient iteration residual in real time and using this value to perform zoning processing on the detection region. For each different region, abnormal nodes beyond the threshold range are removed and filled, eliminating the pollution of non-atmospheric intrinsic noise to the inversion link, thereby significantly improving the robustness of the visibility inversion result, avoiding the abnormal decrease in visibility caused by the accumulation of outliers in the traditional method, and at the same time being able to more accurately and comprehensively reflect the visibility data before and after the abnormal points.

[0006] The principle of outlier rejection in the present invention lies in: constructing a transformation function (where \(i = 0, 1, \cdots, n - 1\), and \(n\) is the number of lidar sampling points) to achieve the conversion from the power signal domain to the logarithm - square domain (where is the power of the echo signal detected by the lidar, is the corresponding distance), by calculating the difference sequence of adjacent transformation values , store (where \(i = 0, 1, \cdots, n - 1\), and \(n\) is the number of lidar sampling points) as one - bit data, denoted as the array , where the superscript \(T\) represents vector transpose, represents the \(N\) - dimensional real - number space. Using the stationary characteristics to perform outlier judgment on : Denote as (\( is the median operation), denote as , if (\(n\) takes values from 3 to 5, and preferably takes 3 according to experiments), it indicates that is abnormal, otherwise it indicates that is normal. Based on the robust statistics constructed by the median \(med\) and \(mad\), hierarchically partition the detection - distance space to obtain the regions before and after mutation. After separately re - performing the above outlier judgment on the regions before and after mutation, the visibility data before and after mutation are obtained through the following operations: When and only when the consecutive difference values and are both abnormal, determine that the original power value is an outlier, the corresponding distance is , remove this outlier and fill in the echo - signal power value suitable for the current data at the corresponding position, and the filled - in corrected echo - signal power data is . Take the updated echo - signal power data as the input data, use the Klett iterative method to calculate the extinction coefficients and visibility data of the two regions before and after mutation respectively, and output the mutation - position interval at the same time.

[0007] To achieve the above - mentioned purpose, the present invention adopts the following technical solution:

[0008] An iterative inversion method for lidar - measured visibility with outlier self - correction, comprising the following steps:

[0009] Step S1, use a lidar to measure the horizontal atmospheric optical properties and obtain the measurement data of the lidar;

[0010] Step S2, correct this set of data;

[0011] Step S3: Partition the data and perform outlier judgment and outlier removal on this set of data, including: identifying outliers by calculating the median absolute deviation of the differences between adjacent data; according to the outlier discrimination result, removing outliers and updating the echo signal power and distance array;

[0012] Step S4: Obtain the extinction coefficient using the slope method;

[0013] Step S5: Take the extinction coefficient calculated by the slope method as the initial value of the first iteration, substitute it into the Klett formula, and perform iteration until the difference between a certain result and the corresponding iteration initial value is less than the set threshold, then stop the iteration;

[0014] Step S6: Calculate the atmospheric visibility according to the relationship between the extinction coefficient and the visibility.

[0015] Beneficial effects:

[0016] On the premise of inheriting the advantage of the Klett algorithm in comprehensively reflecting the real atmosphere, this invention adaptively removes outliers and performs partition processing for the entire data space, breaking through the inherent defects of the traditional local screening method, and accurately reflecting the visibility with obvious differences before and after mutation. By introducing a dynamic threshold generator based on the median absolute deviation of the residual distribution, real-time outlier discrimination and filling are performed on the discrete sampling points of the lidar inversion link, and interference nodes such as instrument transient noise and local foreign object occlusion are accurately removed. Based on the characteristics of the lidar original echo signal data, this invention transforms the problem of abnormal detection of the power signal under multiple interference couplings into the problem of outlier identification of stable distribution parameters, thereby effectively suppressing the situation where the inversion value of visibility is abnormal caused by the transmission of abnormal data in the traditional algorithm, as well as the possible waste of measurement data and incomplete reflection of the atmosphere caused by the local screening method, and systematically ensuring the stability of the atmospheric detection results. Description of the Drawings

[0017] Figure 1 It is a flowchart of an iterative inversion method for lidar measurement of visibility with self-correction of outliers according to the present invention;

[0018] Figure 2 It is a schematic diagram of real echo signal data provided by an embodiment of the present invention. Detailed Embodiment

[0019] In order to make the purpose, technical solution and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention and are not used to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other. Such as Figure 1As shown in the figure, an iterative inversion method for lidar measurement visibility with outlier self-correction according to the present invention includes the following steps:

[0020] Step S1: Receive the echo signal power data generated by the backscattering of atmospheric molecules, and strictly establish a one-to-one mapping relationship between the power value (unit: watt, W) and the detection path distance (unit: kilometer, km); the inversion data selection threshold is set as the calibration distance point when the spatial overlap factor is equal to 1. The present invention only performs extinction coefficient inversion on the discrete sampling sequence of the echo signal corresponding to the effective detection interval where the overlap factor is stably 1 on the distance axis.

[0021] Step S2: Perform dark count correction on the received echo signal power data. The dark count is the data value received by the detector when there is no lidar receiving backscattering. Based on the effective photon signal data set pre-calibrated by the dark count, systematically integrate it and construct it into an array structure The detection distance data is constructed into an array structure where represents the value of the echo signal power detected for the i-th time after dark count correction, represents The corresponding detection distance. i = 0, 1…n, where n is the number of lidar sampling points.

[0022] Step S3: Perform partitioning, and perform outlier judgment and outlier rejection on this set of data, including: identifying outliers by calculating the median absolute deviation of the differences between adjacent data; according to the outlier discrimination result, rejecting outliers and updating the echo signal power and distance arrays, including:

[0023] Step S3.1: Perform distance correction processing on the received echo signal power data and distance data. Let , and store (i = 0, 1…n - 1, where n is the number of lidar sampling points) as one-bit data, denoted as array , where the superscript T represents vector transpose, represents the N-dimensional real number space, calculate the median of the array , median( , denoted as med; calculate the absolute deviation between each data point in the array and med, ; calculate the median of these absolute deviations, , denoted as mad; use mad as the threshold to identify outliers; for each value in the array, if > mad n (where n ranges from 3 to 5, preferably 3 according to experiments; i ), indicating abnormality; otherwise it is a normal value.

[0024] Step S3.2: Based on the robust statistic constructed from the median med and mad obtained in S3.1, hierarchically partition the detection distance space, including:

[0025] 1. Division of the main stable region: Traverse forward the array, and identify ten consecutive more than half of which are normal values. Denote the distance corresponding to the first normal value as . Traverse backward, and also identify ten consecutive more than half of which are normal values. Mark the distance corresponding to the last normal value as , thereby determining the main stable region = ;

[0026] 2. Division of the secondary region and the mutation region: Width of the forward remaining interval: , width of the backward remaining interval: . If , recalculate the med and mad values within the region, and perform forward / backward traversal according to the above process to determine the secondary region , and the mutation position interval = , ; conversely, if , then perform the same operation on to determine the secondary region , and the mutation position interval = , . Obtain the main and secondary stable regions and , as well as the mutation position interval .

[0027] Step S3.3: Perform the operations of Step S3.1 on the main and secondary stable regions and obtained in Step S3.2 respectively. After performing abnormality judgment on each item of the array within the two regions, perform the operation of filling abnormal values within the two regions. The following method will be described in detail for the main stable region, and the implementation process for the secondary region is the same. For the main stable region = ​​Echo signal power data and distance data, distance correction ( 、 ) processing, if and only if the continuous difference value and are both abnormal, determine that the original power value is an abnormal point, and the corresponding distance is . Remove this abnormal value and fill in the echo signal value suitable for the current data at the corresponding position. The corrected echo signal power data filled in is . Region Based on the effective photon signal data set after outlier self-correction processing, systematically integrate it and construct it into an array structure At the same time, integrate and construct the processed detection distance data set into an array structure In. Similarly, region The processed effective photon signal data set is systematically integrated and constructed into an array structure At the same time, integrate and construct the processed detection distance data set into an array structure In.

[0028] Step S4, perform least squares fitting on the 、 and data output by S3.3 respectively, and the slopes obtained are denoted as and , and The initial extinction coefficients of the two regions are respectively: 、 .

[0029] Step S5, calculate the final extinction coefficients and of the two regions and respectively, and the formula is as follows:

[0030] ;

[0031] where is the upper limit of the range distance calculated by the previous slope method; is the currently calculated distance, ; is the distance correction signal at ; is the distance correction signal at ; k represents the extinction backscattering ratio, and its value range is 0.67-1, usually taking the value of 1; the initial is the initial extinction coefficient calculated in step S5 ; For the extinction coefficient at a detection distance of .

[0032] Calculate respectively the extinction coefficients of each sampling point on the path, and average them to obtain the average extinction coefficient value; when the value of the average extinction coefficient is greater than a predetermined threshold value compared with the initial value in the iterative process, take as the initial value for the next iteration , and repeat the iteration. When the value of the average extinction coefficient is less than the predetermined threshold value compared with the initial value in the iterative process, that is , then stop the iteration and output the current average extinction coefficient , which is the final extinction coefficient, The final extinction coefficient of the area is denoted as , The final extinction coefficient of the area is denoted as .

[0033] Step S6, substitute the final extinction coefficient obtained in the previous step and into the extinction coefficient visibility equation to obtain and , and the equations are as follows:

[0034] ;

[0035] where V is the visibility; is the atmospheric extinction coefficient, which is the and obtained in S5, ; is the laser wavelength; is the correction coefficient, V>50km, q = 1.6; 6km<V<50km, q = 1.3; V<6km, q = 0.585 .

[0036] As Figure 2 shown is the schematic diagram of the real echo signal data provided by the embodiment of the present invention. The ordinate S is the value of the echo signal power after distance correction, and the abscissa r is the detection distance.

[0037] The output result of the present invention for calculating the real echo signal data is: the extinction coefficient in the main stable area (1.1km - 1.7025km) is 0.0789714, and the visibility is 51.7188km; the secondary stable area The extinction coefficient for (2.5725 km - 5 km) is 0.501226, and the visibility is 8.14861 km. The mutation region is 1.7025 km - 2.5725 km.

[0038] It should be noted that in the present invention, relational terms such as "first" and "second" are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprising", "including" or any other variant thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article or device. Without further limitation, an element defined by the phrase "comprising a..." does not exclude the presence of additional identical elements in the process, method, article or device comprising the said element.

[0039] The above are only specific embodiments of the present invention, enabling those skilled in the art to understand or implement the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined in the present invention can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown in the present invention, but will conform to the widest scope consistent with the principles and novel features of the present invention.

Claims

1. An iterative inversion method for lidar measurement visibility with outlier self-correction, characterized in that, It includes the following steps: Step S1: Use lidar to measure the horizontal or slant-path atmospheric optical properties and obtain the measurement data of the lidar; Step S2: Correct this set of data; Step S3: Perform zoning and judge and update the data for this set of data, including: identifying outliers by calculating the median absolute deviation of the differences between adjacent data; according to the outlier discrimination result, perform zoning, and then remove outliers from each region and update the echo signal power and distance array; Step S4: Obtain the extinction coefficient using the slope method; Step S5: Take the extinction coefficient calculated by the slope method as the initial value of the first iteration, substitute it into the Klett formula, and perform iteration until the difference between a certain result and the corresponding iteration initial value is less than the set threshold, and then stop the iteration; Step S6: Calculate the atmospheric visibility according to the relationship between the extinction coefficient and the visibility; 2. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 1, characterized in that The said Step S1 includes: receiving the laser signal power backscattered by atmospheric molecules, and this data should correspond one by one to the detection point distance.

3. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 1, characterized in that The said Step S2 includes: performing dark count correction on the received echo signal power data.

4. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 1, wherein The said Step S3 includes: Step S3.1 performs distance correction on the received echo signal power data and distance data processing; Step S3.2: Based on the robust statistic constructed by the median (med) and mad obtained in Step S3.1, perform hierarchical zoning on the detection distance space; Step S3.3 performs the operation of Step S3.1 on the primary and secondary stable regions obtained in Step S3.2 and respectively. After performing anomaly judgment on each item of the array within the two regions an anomaly value filling operation is performed on each of the two regions respectively.

5. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 4, characterized in that The said Step S3.1 includes: Let , store as a single-bit data, denoted as array , where \(i = 0, 1, \ldots, n - 1\), \(n\) is the number of lidar sampling points, the superscript \(T\) represents vector transpose, represents the \(N\)-dimensional real space. Calculate the median of the array , that is, \(median( ), denoted as \(med\); calculate the absolute deviation between each data point in the array and \(med\), ; calculate the median of these absolute deviations, that is , denoted as \(mad\); use \(mad\) as a threshold to identify outliers; for each value in the array, if > \(mad n\), \(n\) is taken as 3 according to the experiment; \(i ), it indicates that it is an outlier; otherwise it is a normal value.

6. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 5, characterized in that The said Step S3.2 includes: Perform the division of the main stability region and traverse forward the array. If ten consecutive more than half of which are normal values, record the distance corresponding to the first normal value as ; Traverse backward. Similarly, if ten consecutive more than half of which are normal values, mark the distance corresponding to the last normal value as , thereby determining the main stability region = ; Perform secondary region and mutation region division, the width of the front remaining interval is , the width of the rear remaining interval is ; If , recalculate the med and mad values within the region, perform forward / backward traversal according to the above process to determine the secondary region , the mutation position interval = , ; Otherwise, if , then perform the same operation on to determine the secondary region , the mutation position interval = , , obtain the primary and secondary stable regions and , and the mutation position interval ;​​ For the above region and , perform the operation of S3.1 again respectively; calculate the median of ΔS, that is, median(ΔS); calculate the absolute deviation between each data point in ΔS and the median median(ΔS), |ΔS j -median(ΔS)|; calculate the median of these absolute deviations, that is, obtain MAD, MAD = median(|ΔS j -median(ΔS)|); use MAD as the threshold to identify outliers. If >mad n, where n is taken as 3 according to the experiment; i ), it indicates abnormal; otherwise it is a normal value.

7. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 6, wherein In the step S3.3, for the echo signal power data and distance data of the main stable region = , perform distance correction ( 、 )processing. When and only when the continuous difference values and are both abnormal, determine that the original power value is an abnormal point, the corresponding distance is . Remove this abnormal value and fill in the echo signal value suitable for the current data at the corresponding position. The corrected echo signal power data after filling is ; Region Based on the effective photon signal data set after outlier self-correction processing, systematically integrate it and construct it into an array structure , and at the same time integrate the processed detection distance data set into an array structure ; Similarly, region The processed effective photon signal data set is systematically integrated and constructed into an array structure , and at the same time integrate the processed detection distance data set into an array structure .

8. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 7, characterized in that The said Step S4 includes: For the outputs of step S3.3 and and data are respectively subjected to least square fitting, and the obtained slopes are respectively denoted as and , and The initial extinction coefficients of the two regions are respectively: , .

9. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 8, characterized in that The said Step S5 includes: Calculated by the Klett iterative method and the final extinction coefficients of the two regions and , the formula is as follows: ； Among them, is the upper limit of the range distance calculated by the slope method in the previous step; is the currently calculated distance, ; is the distance correction signal at; is the distance correction signal at; k represents the extinction backscattering ratio, and its value range is 0.67 - 1, usually taking the value of 1; the initial is the initial extinction coefficient calculated in step S5 ; is the extinction coefficient at the detection distance of ; Calculate separately the extinction coefficients of each sampling point on the path, and average them to obtain the average extinction coefficient value; when the average extinction coefficient value and the initial value in the iteration process are greater than a predetermined threshold, take as the initial value for the next iteration , and repeat the iteration; when the average extinction coefficient value and the initial value in the iteration process are less than a predetermined threshold, that is , then stop the iteration and output the current average extinction coefficient , which is the final extinction coefficient, the final extinction coefficient of the region is denoted as , the final extinction coefficient of the region is denoted as .

10. The iterative inversion method for lidar measurement visibility with outlier self-correction according to claim 1, wherein The said Step S6 includes: The final extinction coefficient obtained in this step S5 and are substituted into the following equation to obtain and which is ； Among them, V is the visibility; is the atmospheric extinction coefficient, which is obtained from S5 and , ; is the laser wavelength; is the correction coefficient, when V > 50km, q = 1.6; when 6km < V < 50km, q = 1.3; when V < 6km, q = 0.585 .

Images