CN113325440B - Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition - Google Patents

Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition Download PDF

Info

Publication number
CN113325440B
CN113325440B CN202110491517.5A CN202110491517A CN113325440B CN 113325440 B CN113325440 B CN 113325440B CN 202110491517 A CN202110491517 A CN 202110491517A CN 113325440 B CN113325440 B CN 113325440B
Authority
CN
China
Prior art keywords
signal
data
height
aerosol
cloud
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110491517.5A
Other languages
Chinese (zh)
Other versions
CN113325440A (en
Inventor
殷振平
易帆
张天澈
周军
何芸
柳付超
张云飞
熊梓稀
谢文超
王磊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202110491517.5A priority Critical patent/CN113325440B/en
Publication of CN113325440A publication Critical patent/CN113325440A/en
Application granted granted Critical
Publication of CN113325440B publication Critical patent/CN113325440B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/95—Lidar systems specially adapted for specific applications for meteorological use
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88—Lidar systems specially adapted for specific applications
    • G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Electromagnetism (AREA)
  • Optical Radar Systems And Details Thereof (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

The invention relates to a method and a system for inverting polarized laser radar data based on image recognition and signal feature decomposition. The method comprises the steps of meteorological data downloading, radar data preprocessing, cloud layer identification, reference height selection and data inversion, wherein the meteorological data downloading automatically downloads radiosonde data or model reanalysis data related to observation time for calculating the scattering intensity of atmospheric molecules, the atmospheric molecule scattering intensity is compared with laser radar data processed by a laser radar data preprocessing module, after a section and an aerosol clean area are respectively removed by cloud layer identification and reference height selection, the residual data enter the radar data inversion, and information such as aerosol extinction coefficient and height distribution of particle depolarization ratio is obtained. The invention solves the difficulty of manually identifying cloud layers and aerosol clean areas in the data processing of the atmospheric sounding laser radar, fully automates the traditional data analysis process, and can meet the requirements of various practical applications on the calculation efficiency and the product quality.

Description

Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition
Technical Field
The invention belongs to the field of laser atmospheric remote sensing inversion, and particularly relates to a polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition.
Background
And (3) inverting the data of the atmospheric sounding laser radar, and inputting a height area with negligible aerosol concentration and a section for filtering cloud influence. For the identification of aerosol clean zones, there are two main types of processes: (1) directly adopting signals at the height of the stratosphere; (2) the Ra yleigh fit determination method (bairs et al, 2016) utilizes the comparison of the lidar signal and the atmospheric molecular scattering signal. Both methods have their own disadvantages, and method 1 is difficult to use for inversion of daytime lidar data. Because the signal during the day is subject to the signal-to-noise ratio, the signal at stratospheric height is dominated by noise. Method 2 can be used for daytime aerosol clean zone screening, but because method 2 requires a clean zone search window assuming a certain width, the method fails when the clean zone is narrow. Moreover, since method 2 is a height-by-height comparison method, it is computationally inefficient to traverse all height gates.
Typical methods for cloud layer identification are: signal gradient methods (Wang and Sassen,2001), wavelet transform methods (Baars et al, 2016), Douglas-Peucker methods (Gong et al, 2011), machine learning methods (Binietoglou et al, 2018), and signal feature enhanced image recognition methods (Zhao et al, 2014). Wherein, the signal gradient method and the wavelet change method both need the calibrated radar signal profile. In practice, the system efficiency of radar varies with the laser energy, the system transmittance, and the quantum efficiency of the photomultiplier tube, and thus cannot be used in an automated data processing scheme. The Douglas-Peucker method can decompose the aerosol layer and cloud layer in the signal profile, but the decomposition is very susceptible to signal noise, so the detection effect on thin cloud and daytime cloud layer is not good (Mao et al, 2013). The machine learning method can overcome the above difficulties, but the machine learning method needs a pre-classified data set for model training, and thus needs extra workload.
In addition, the laser radar data has important effects on understanding physical processes occurring in the atmosphere, real-time environmental monitoring, later-stage evaluation of aerosol content change and the like. Conventional lidar data processing is typically based on manual analysis, which is time consuming. This solution is not desirable when the data volume is large or when multiple sites are observed jointly, so an automatic inversion method is needed to solve these problems. Foreign atmospheric lidar automatic inversion methods are generally based on raman lidar systems and therefore do not address the inversion of elastic lidar data (Baars et al, 2016; Thorsen and Fu, 2015). The satellite-borne CALIPO data processing algorithm is successfully applied to satellite-borne elastic radar data processing, but the ground-based radar and the satellite-borne radar have different observation angles, so that the inversion algorithm of the satellite-borne radar cannot be directly applied to ground-based data processing (Winker et al, 2009). The processing of ground-based elastic lidar data still needs to be addressed by other methods.
Disclosure of Invention
In order to solve the data inversion problem, the invention provides the automatic laser radar inversion method combining the image recognition algorithm and the signal characteristic decomposition, can solve the problems of low efficiency in manual inversion, difficulty in selecting reference heights in low signal-to-noise ratio in daytime, dependence of the traditional cloud layer recognition algorithm on the efficiency of a laser radar system, fixed length of a search window of the traditional reference height recognition algorithm, low calculation efficiency and the like, can realize the characteristics of automatic recognition of cloud layers and reference heights, noise filtering, intelligent recognition of signal characteristics, high calculation efficiency and the like, and has the advantages of no need of any manual intervention, simplicity in setting, high universality and the like.
In order to achieve the above purpose, the invention adopts the following technical scheme to realize: a structure diagram of the whole method is shown in figure 1, and the method comprises polarized laser radar data preprocessing, cloud layer recognition, non-cloud section segmentation, reference height recognition, data inversion and inversion result output, wherein meteorological data input is required in the reference height recognition and data inversion processes.
Step 1, preprocessing polarization laser radar data to correct a pulse accumulation effect in the laser radar data;
step 2, identifying cloud layers detected in the laser radar single-section data;
step 3, screening out all signal profiles containing cloud layers to obtain profiles without cloud layers, and dividing the non-cloud profiles into continuous subsets with fixed time lengths;
step 4, inputting meteorological data;
step 5, calculating a molecular scattering signal by combining the meteorological data input in the step 4, identifying the reference height of the cloudless profile segmentation radar signal in the step 3, and finding out a characteristic height section closest to the molecular scattering signal as the reference height;
step 6, combining the meteorological data input in the step 4, and performing inversion from the continuous cloud-free signal profile in the step 5 to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio;
and 7, finally outputting the inversion result into a scientific data format file with a standard format.
Further, in step 1, the pulse pile-up effect in the laser radar data is corrected, and the correction formula is as follows:
Figure BDA0003052433420000031
wherein C' is the corrected photon counting rate; c is the photon count rate affected by the pulse pile-up effect; τ is the dead time; after pulse accumulation effect correction, the two-channel data of photon counting and analog sampling are spliced, and the splicing formula is as follows:
Figure BDA0003052433420000032
where C "represents the spliced signal, m represents the weight of the analog sampled signal, s is the slope, o is the offset, A is the intensity of the analog sampled signal, C min And C max Respectively the range of the photon counting rate of the splicing;
and finally, removing the background signal from the spliced signal, wherein the calculation formula of the background signal is as follows:
Figure BDA0003052433420000033
wherein P is bg For background signal, k is the number of range gates, n is the initial range gate of the background signal, P (z) i ) Is the signal strength at the ith height gate.
Further, the cloud layer identification method in step 2 is to identify k 1 The multiple signal uncertainty is used for distinguishing jitter generated by natural change of signals from echo signal jitter generated by a real characteristic layer, and adjacent signal deviation larger than k is reserved from the near end to the far end of the signals in sequence 1 A point of multiple signal uncertainty, wherein the signal uncertainty is calculated by the formula:
Figure BDA0003052433420000034
the same process is carried out from the far end to the near end again according to the reverse sequence, the signals obtained twice are averaged to obtain a semi-discretized signal profile, and the signals obtained twice are averaged to obtain a semi-discretized signal profile; then, a histogram equalization algorithm is adopted to detect a height interval higher than a reference signal as a characteristic layer, and finally, a threshold value is utilized to judge whether the maximum value and the minimum value of the signal in the layer exceed k or not 2 To distinguish whether the feature layer is a cloud layer.
Further, the meteorological data input in step 4 includes data downloading and data interpolation, wherein the data downloading is used for downloading real-time or historical radiosonde data, 1 ° × 1 ° meteorological field data output by a global data assimilation system GDAS, or a meteorological field output by a re-analysis data set ERA5 of a european mesoscale meteorological forecast model of the fifth generation, and the data interpolation is used for interpolating a global-grid meteorological field onto observation grid points of the laser radar by a bilinear interpolation method.
Further, the reference height identification in the step 5 comprises signal feature decomposition and Rayleigh fit judgment, wherein the signal feature decomposition divides a radar signal vertical section into height sections with the same discrete signal features by using a Douglas-Peucker algorithm, then initial detection is carried out on the divided signal feature height sections, the height sections comprise the vertical thickness of the feature height sections and the positions of the feature height sections, then Rayleigh fit judgment is carried out on the feature height sections meeting the height range and minimum thickness limitation, and the feature height section closest to the molecular scattering signal is found out to be used as the reference height; wherein the molecular scattering signal is calculated according to meteorological data, and the calculation formula is as follows:
Figure BDA0003052433420000041
wherein P is m Represents the molecular scattering signal, beta m And alpha m Representing the molecular backscattering coefficient and extinction coefficient; the molecular backscattering coefficient and extinction coefficient can be calculated by the following formulas:
Figure BDA0003052433420000042
Figure BDA0003052433420000043
wherein p (z) is atmospheric pressure in units of hectopascal, T (z) is atmospheric temperature in units of Kelvin, λ 0 The wavelength of the excitation light is in nanometers.
Further, the specific implementation manner of obtaining the aerosol extinction coefficient, the backscattering coefficient and the particle depolarization ratio by inversion in the step 6 is as follows;
for an emission wavelength of λ 0 The radar equation describing the received signal of the elastic scattering lidar of (1) is shown in equation (8):
Figure BDA0003052433420000044
where P (z) is the strength of the radar echo signal at height z, C is the receiving efficiency of the entire system, O (z) is the overlap factor of the system, and β a Is the backscattering coefficient, alpha, of the aerosol a Z' represents an integral variable for the extinction coefficient of the aerosol; the backscattering coefficient and extinction coefficient of aerosol cannot be solved by an equation, and the traditional solution method is to assume that the backscattering coefficient and the extinction coefficient satisfy a linear relationship, as shown in equation (9):
Figure BDA0003052433420000051
wherein S is a Is the aerosol extinction-backscattering ratio, also known as the lidar ratio; based on this relationship, a numerical solution for the backscattering coefficient of the aerosol can be obtained by solving the bernoulli equation as follows:
Figure BDA0003052433420000052
wherein z is ref For the reference height, obtained by step 5, z "represents an integral variable;
after obtaining the backscattering coefficient of the aerosol, the particle depolarization ratio can be obtained by inverting the following formula:
Figure BDA0003052433420000053
Figure BDA0003052433420000054
where K is the ratio of the system efficiencies of the parallel and perpendicular channels, calibrated by the Δ 90 method, P || And P ⊥ Signal strength, delta, for parallel and vertical channels, respectively m Is the molecular depolarization ratio, a quantity related only to the receive channel filter bandwidth and atmospheric temperature.
The invention also provides a polarization laser radar data inversion system based on image recognition and signal characteristic decomposition, which comprises the following modules:
the preprocessing module is used for preprocessing the polarization laser radar data so as to correct the pulse accumulation effect in the laser radar data;
the cloud layer identification module is used for identifying cloud layers detected in the laser radar single-section data;
the cloud-free section segmentation module is used for screening out all signal sections containing cloud layers to obtain sections without the cloud layers, and segmenting the cloud-free sections into continuous subsets with fixed time lengths;
the meteorological data input module is used for inputting meteorological data;
the reference height identification module is used for calculating a molecular scattering signal by combining the input meteorological data, carrying out reference height identification on the cloudless profile segmentation radar signal and finding out a characteristic height section closest to the molecular scattering signal as a reference height;
and the data inversion module is used for inverting from the continuous non-cloud signal profile by combining the input meteorological data to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio.
And the inversion result output module is used for outputting the inversion result as a scientific data format file with a standard format.
The beneficial effects of the invention are: compared with the prior art, the method can solve the problems that the efficiency of a manual inversion process is low, the reference height in the daytime cannot be identified, the cloud layer identification algorithm depends on a laser radar calibration result, and the calculation efficiency in the reference height identification process is low. The method has the advantages of obtaining inversion products from the polarization laser radar data in a high-efficiency and automatic mode.
Drawings
The invention is further illustrated with reference to the following figures and examples.
FIG. 1 is an algorithmic work flow diagram of the present invention;
FIG. 2 is a laser radar signal pre-processing presentation diagram;
fig. 3 is a diagram showing an identification process of actual data in a cloud layer identification process, (a) the diagram shows typical non-cloud and cloud laser radar signal profiles, wherein a solid line shows an actually measured non-cloud radar signal profile, a dotted line shows a signal profile with a low-layer liquid water cloud and a high-layer rolling cloud, and a signal marked by a thick line segment is a cloud layer identified by the cloud layer identification algorithm; (b) represents the signal profile after each half-discretization, wherein the dotted line represents three times the corresponding signal uncertainty; (c) the solid and dashed implementations represent the signal profile equalized by the value diagram, respectively, and the dashed gray lines represent the corresponding standard signals;
fig. 4 is a diagram showing the actual data processing results of the lidar of the present invention throughout the day, and (a) is a time height diagram of the range correction signal. Wherein the black solid line is the aerosol backscattering coefficient obtained by inversion, and the bold solid line is the reference height range obtained automatically; (b) the cloud layer identification result is obtained; the signal profile containing the cloud layer consists of a middle strip.
Fig. 5 is a diagram showing the identification process of actual data in the reference height identification process, and (a) represents the measured distance correction signal. The black dot-dash line represents signal uncertainty; (b) the middle dots distinguish new number segments which are identified by a Douglas-Peucker algorithm and have the same characteristics; (c) the signal marked by the medium black is the signal segment closest to the molecular scattering signal;
Detailed Description
In order to make the technical means, features and functions of the invention easy to understand, the invention is further explained below with reference to the specific drawings.
Step 1, preprocessing polarization laser radar data to correct a pulse accumulation effect in the laser radar data, wherein a correction formula is as follows:
Figure BDA0003052433420000071
wherein C' is the corrected photon counting rate; c is the photon count rate affected by the pulse pile-up effect; τ is the dead time. After pulse accumulation effect correction, the two-channel data of photon counting and analog sampling are spliced, and the splicing formula is as follows:
Figure BDA0003052433420000072
where C "represents the spliced signal, m represents the weight of the analog sampled signal, s is the slope, o is the offset, A is the intensity of the analog sampled signal, C min And C max Respectively the range of the photon counting rate of the splicing. The weights, slopes, and offsets may be calculated from the signal strengths, typically as known quantities. The splicing range is related to the signal of the used acquisition instrument, and is generally selected to be 1-50 MHz. And finally, removing the background signal from the spliced signal, wherein the calculation formula of the background signal is as follows:
Figure BDA0003052433420000073
wherein P is bg For background signal, k is the number of range gates, n is the initial range gate of the background signal, P (z) i ) Is the signal strength at the ith height gate. Generally, n is a height gate at a very far end, and the effective echo signal at such a height is extremely weak, and almost all received signals are background signals.
And 2, identifying the cloud layer for identifying the cloud layer detected in the laser radar single-section data. The identification method is characterized in that triple signal uncertainty is used for distinguishing jitter generated by natural change of a signal from echo signal jitter generated by a real characteristic layer, the signal is sequentially kept from a near end to a far end at a point of signal uncertainty of which the adjacent signal deviation is more than three times (other values can be taken, generally more than or equal to 3, and if the value is too small, a noise signal is included), wherein the signal uncertainty can be calculated by the following formula:
Figure BDA0003052433420000074
and the same process is performed from the far end to the near end again according to the reverse sequence, the signals obtained twice are averaged to obtain a semi-discretized signal profile, and the signals obtained twice are averaged to obtain the semi-discretized signal profile. And then detecting a height interval higher than the reference signal by adopting a histogram equalization algorithm commonly used in image processing as a characteristic layer. And finally, judging whether the maximum value and the minimum value of the signals in the layer exceed four (other values can be also taken, and 4 is the most reasonable threshold) by utilizing the threshold to distinguish whether the characteristic layer is a cloud layer.
And 3, screening all the signal profiles containing the cloud layer to obtain the profiles without the cloud layer, and dividing the non-cloud profiles into continuous subsets with fixed time lengths.
And 4, meteorological data input, including data downloading and data interpolation, wherein the data downloading is used for downloading real-time or historical radiosonde data, meteorological field data of 1 degree multiplied by 1 degree output by a global data assimilation system GDAS or a meteorological field output by a re-analysis data set ERA5 of a fifth generation European mesoscale meteorological forecasting model (ECMWF), and the data interpolation is to interpolate the global gridded meteorological field to observation grid points of the laser radar by a bilinear interpolation method.
Step 5, calculating a molecular scattering signal by combining the meteorological data input in the step 4, identifying the reference height of the cloudless profile segmentation radar signal in the step 3, and finding out a characteristic height section closest to the molecular scattering signal as the reference height; the reference altitude identification includes signal feature decomposition and Rayleigh fit determination. The signal characteristic decomposition divides a radar signal vertical section into height sections with the same signal characteristics by using a Douglas-Peucker algorithm, then performs initial detection on the height sections of the divided signal characteristics, including the vertical thickness of the height sections of the characteristics and the positions of the height sections of the characteristics, then performs Rayleigh fit judgment on the height sections of the characteristics meeting the limitation of height range and minimum thickness, and finds out the height section of the characteristics closest to the molecular scattering signal as a reference height. Wherein the molecular scattering signal is calculated according to meteorological data, and the calculation formula is as follows:
Figure BDA0003052433420000081
wherein P is m Represents the molecular scattering signal, beta m And alpha m Representing the molecular backscattering coefficient and extinction coefficient. The molecular backscattering coefficient and extinction coefficient can be calculated by the following formulas:
Figure BDA0003052433420000082
Figure BDA0003052433420000083
where p (z) is atmospheric pressure in units of hectopascal. T (z) is the atmospheric temperature in Kelvin. Lambda 0 The wavelength of the excitation light is in nanometers.
And 6, performing data inversion to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio by performing inversion on the continuous cloud-free signal profile obtained in the step 5.
For an emission wavelength of λ 0 The radar equation describing the received signal of the elastic scattering lidar of (1) is shown in equation (8):
Figure BDA0003052433420000091
wherein P (z) is the strength of radar echo signal at height z, C is the receiving efficiency of the whole system, and O (z) is that of the systemOverlap factor, beta a Is the backscattering coefficient, alpha, of the aerosol a Is the aerosol extinction coefficient. The backscattering coefficient and extinction coefficient of aerosol cannot be solved by an equation, and the traditional solution method is to assume that the backscattering coefficient and the extinction coefficient both satisfy a linear relationship, as shown in equation (9):
Figure BDA0003052433420000092
wherein S is a Is the extinction-backscattering ratio of the particles, also known as the laser radar ratio. Based on this relationship, a numerical solution for the backscattering coefficient of the aerosol can be obtained by solving the bernoulli equation as follows:
Figure BDA0003052433420000093
wherein z is ref For reference height, it can be obtained by step 5. After obtaining the aerosol backscattering coefficient, the aerosol extinction coefficient can be calculated by equation (9). Because there is a large difference between the scattering characteristics of clouds in the atmosphere and aerosols (Ansmann et al, 1992), and the scattering signals of cloud layers are affected by multiple scattering of cloud layers (Hogan,2006), in order to satisfy the assumption that the ratio of lidar is constant in the inversion process, all the sections containing cloud layers need to be filtered, and the sections without clouds can be obtained through step 2 and step 3. After obtaining the backscattering coefficient of the aerosol, the particle depolarization ratio delta a (z) can be obtained by inversion of the following equation:
Figure BDA0003052433420000094
Figure BDA0003052433420000095
where K is the ratio of the system efficiencies of the parallel and perpendicular channels, can be calibrated by the Δ 90 method (Freu)dent haler,et al.,2009)，P || And P ⊥ The signal strength of the parallel and vertical channels, respectively. Delta m Is the molecular depolarization ratio, a quantity related only to the receive channel filter bandwidth and atmospheric temperature.
And 7, finally outputting the inversion result into a scientific data format file with a standard format.
The overall working principle of the algorithm of the invention is as follows: the laser radar collected signal is input into the laser radar data processing module for preliminary signal correction and signal splicing, the result is output into the cloud layer identification module, the sections with cloud layer interference are identified and removed from the original sections, then inputting the cloudless profile into a cloudless profile segmentation module, segmenting the separated signal profile into a plurality of continuous signal profiles with a fixed time length, the signals of the profiles are accumulated in time and then input into a reference height identification module, the areas closest to the molecular scattering signals are identified as reference heights through comparison with the molecular scattering signals, and then inputting the reference height and the signal into a data inversion module to calculate an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio, and finally outputting the aerosol extinction coefficient, the backscattering coefficient and the particle depolarization ratio as a scientific data format file through an inversion result output module.
The embodiment of the invention also provides a polarization laser radar data inversion system based on image recognition and signal characteristic decomposition, which comprises the following modules:
the preprocessing module is used for preprocessing the polarization laser radar data so as to correct the pulse accumulation effect in the laser radar data;
the cloud layer identification module is used for identifying cloud layers detected in the laser radar single-section data;
the cloud-free section segmentation module is used for screening out all signal sections containing cloud layers to obtain sections without the cloud layers, and segmenting the cloud-free sections into continuous subsets with fixed time lengths;
the meteorological data input module is used for inputting meteorological data;
the reference height identification module is used for calculating a molecular scattering signal by combining the input meteorological data, carrying out reference height identification on the cloudless profile segmentation radar signal and finding out a characteristic height section closest to the molecular scattering signal as a reference height;
and the data inversion module is used for inverting from the continuous non-cloud signal profile by combining the input meteorological data to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio.
And the inversion result output module is used for outputting the inversion result as a scientific data format file with a standard format.
The specific implementation manner and the steps of each module correspond, and the invention is not described.
The first embodiment is as follows: the collected raw data of the 532nm polarization lidar includes two types of photon counting data and analog sampling data, and the two types of data have different dynamic ranges as shown in fig. 2. Photon count data tends to saturate with low-level echo data because there is a pulse pile-up response when the signal is too strong. The analog sampling data has higher electric background noise, so the display effect on the echo signal which is weak at high altitude is poor. In order to facilitate the subsequent radar signal processing, the original laser radar signal passes through a polarization laser radar data preprocessing module, the photon counting pulse accumulation effect is corrected, and the two paths of output signals are integrated into a path of splicing signal which is saturated and has a wider dynamic range. The integrated signal is subjected to cloud interference detection through a cloud layer identification algorithm module, and typical identification results of sections with clouds and sections without clouds are shown in fig. 3. Fig. 3(a) shows the results of the signal profile affected by the high-level cirrus and the low-level liquid cloud (liquid cloud) and the cloudless signal profile, and after noise filtering and semi-discretization, the noise in the two profiles is filtered, and the result is shown in fig. 3 (b). Finally, the contrast of the signal is improved by a threshold distribution histogram method, and the region where the signal is increased is highlighted in the signal profile, so that a cloud layer and an aerosol layer are highlighted, as shown in fig. 3(c), and finally, the cloud layer is identified by a simple cloud layer and aerosol layer identification algorithm, and the identification result is shown in fig. 4 (b). The results are then divided into successive observation intervals of 1 hour duration, the signal passage times of the individual intervals being accumulated into individual signal profiles of high signal-to-noise ratio, which are fed into a reference height recognition module. In the reference altitude identification module, the altitude segments with signal intensity less than three times the signal uncertainty are first filtered, as shown in fig. 5 (a). The entire signal is then decomposed into different feature height segments using the Dou glas-Peucker algorithm, as shown in fig. 5 (b). The height segment closest to the molecular scattering signal is then identified as the reference height using the Rayleigh fit algorithm, as shown in FIG. 5 (c). The reference height segment and the signal profile are input into an inversion algorithm module, and finally a height profile of the aerosol backscattering coefficient, the extinction coefficient and the particle depolarization ratio is obtained, wherein the backscattering coefficient profile is shown in fig. 4 (a).
It should be understood that portions not set forth in detail in this specification pertain to existing algorithms or techniques.
It should be understood that the above description of the preferred embodiments is given in some detail, and not as a limitation on the scope of the invention, and that various alternatives and modifications can be devised by those skilled in the art without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (5)

1. A method for inverting polarized laser radar data based on image recognition and signal feature decomposition is characterized by comprising the following steps:
step 1, preprocessing polarization laser radar data to correct a pulse accumulation effect in the laser radar data;
step 2, identifying cloud layers detected in the laser radar single-section data;
step 3, screening out all signal profiles containing cloud layers to obtain profiles without the cloud layers, and dividing the non-cloud profiles into continuous subsets with fixed time lengths;
step 4, inputting meteorological data;
step 5, calculating a molecular scattering signal by combining the meteorological data input in the step 4, identifying the reference height of the cloudless profile segmentation radar signal in the step 3, and finding out a characteristic height section closest to the molecular scattering signal as the reference height;
the reference height identification in the step 5 comprises signal feature decomposition and Rayleigh fit judgment, wherein the signal feature decomposition divides a radar signal vertical section into height sections with the same discrete signal features by using a Douglas-Peucker algorithm, then the initial detection is carried out on the height sections of the divided signal features, the height sections comprise the vertical thickness of the feature height sections and the positions of the feature height sections, then the Rayleigh fit judgment is carried out on the feature height sections meeting the limitation of the height range and the minimum thickness, and the feature height section closest to the molecular scattering signal is found out to be used as the reference height; wherein the molecular scattering signal is calculated according to meteorological data, and the calculation formula is as follows:
Figure FDA0003672160490000011
wherein P is m Represents the molecular scattering signal, beta m And alpha m Representing the molecular backscattering coefficient and extinction coefficient; the molecular backscattering coefficient and extinction coefficient can be calculated by the following formulas:
Figure FDA0003672160490000012
Figure FDA0003672160490000013
wherein p (z) is atmospheric pressure in hectopascal, T (z) is atmospheric temperature in Kelvin, λ 0 Is the wavelength of the excitation light in nanometers;
step 6, combining the meteorological data input in the step 4, and performing inversion from the continuous cloud-free signal profile in the step 5 to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio;
the specific implementation mode of obtaining the aerosol extinction coefficient, the backscattering coefficient and the particle depolarization ratio through inversion in the step 6 is as follows;
for an emission wavelength of λ 0 The radar equation describing the received signal of the elastic scattering lidar of (1) is shown in equation (8):
Figure FDA0003672160490000021
where P (z) is the strength of the radar echo signal at height z, C is the receiving efficiency of the entire system, O (z) is the overlap factor of the system, and β a Is the backscattering coefficient of the aerosol, alpha a Is the extinction coefficient of the aerosol, z' represents the integral variable; the backscattering coefficient and extinction coefficient of aerosol cannot be solved by an equation, and the traditional solution method is to assume that the backscattering coefficient and the extinction coefficient satisfy a linear relationship, as shown in equation (9):
Figure FDA0003672160490000022
wherein S is a Is the aerosol extinction-backscattering ratio, also known as the lidar ratio; based on this relationship, a numerical solution for the backscattering coefficient of the aerosol can be obtained by solving the bernoulli equation as follows:
Figure FDA0003672160490000023
wherein z is ref For the reference height, obtained by step 5, z "represents an integral variable;
after obtaining the backscattering coefficient of the aerosol, the particle depolarization ratio can be obtained by inverting the following formula:
Figure FDA0003672160490000024
Figure FDA0003672160490000025
where K is the ratio of the system efficiencies of the parallel and perpendicular channels, calibrated by the Δ 90 method, P || And P ⊥ Signal strength, delta, for parallel and vertical channels, respectively m Is the molecular depolarization ratio, which is a quantity related only to the bandwidth of the optical filter of the receiving channel and the atmospheric temperature;
and 7, outputting the inversion result into a scientific data format file with a standard format.
2. The method of claim 1, wherein the method comprises the following steps: in step 1, the pulse pile-up effect in the laser radar data is corrected, and the correction formula is as follows:
Figure FDA0003672160490000031
wherein C' is the corrected photon counting rate; c is the photon count rate affected by the pulse pile-up effect; τ is the dead time; after pulse accumulation effect correction, the two-channel data of photon counting and analog sampling are spliced, and the splicing formula is as follows:
Figure FDA0003672160490000032
where C "represents the spliced signal, m represents the weight of the analog sampled signal, s is the slope, o is the offset, A is the intensity of the analog sampled signal, C min And C max Respectively the range of the photon counting rate of the splicing;
and finally, removing the background signal from the spliced signal, wherein the calculation formula of the background signal is as follows:
Figure FDA0003672160490000033
wherein P is bg For background signal, k is the number of range gates, n is the initial range gate of the background signal, P (z) i ) Is the signal strength at the ith height gate.
3. The method of claim 1, wherein the method comprises the following steps: the cloud layer identification method in the step 2 is to identify k 1 The multiple signal uncertainty is used for distinguishing jitter generated by natural change of signals from echo signal jitter generated by a real characteristic layer, and adjacent signal deviation larger than k is reserved from the near end to the far end of the signals in sequence 1 A point of multiple signal uncertainty, wherein the signal uncertainty is calculated by the formula:
Figure FDA0003672160490000034
the same process is carried out from the far end to the near end again according to the reverse sequence, the signals obtained twice are averaged to obtain a semi-discretized signal profile, and the signals obtained twice are averaged to obtain a semi-discretized signal profile; then, a histogram equalization algorithm is adopted to detect a height interval higher than a reference signal as a characteristic layer, and finally, a threshold value is utilized to judge whether the maximum value and the minimum value of the signal in the layer exceed k or not 2 To distinguish whether the feature layer is a cloud layer.
4. The method of claim 1, wherein the method comprises the following steps: the meteorological data input in the step 4 comprises data downloading and data interpolation, wherein the data downloading is used for downloading real-time or historical radiosonde data, 1-degree multiplied by 1-degree meteorological field data output by a global data assimilation system GDAS or a meteorological field output by a re-analysis data set ERA5 of a fifth generation European mesoscale meteorological forecasting model, and the data interpolation is to interpolate the global-grid meteorological field to an observation grid point of the laser radar by a bilinear interpolation method.
5. A polarized laser radar data inversion system based on image recognition and signal feature decomposition is characterized by comprising the following modules:
the preprocessing module is used for preprocessing the polarization laser radar data so as to correct the pulse accumulation effect in the laser radar data;
the cloud layer identification module is used for identifying cloud layers detected in the laser radar single-section data;
the cloud-free section segmentation module is used for screening out all signal sections containing cloud layers to obtain sections without the cloud layers, and segmenting the cloud-free sections into continuous subsets with fixed time lengths;
the meteorological data input module is used for inputting meteorological data;
the reference height identification module is used for calculating a molecular scattering signal by combining the input meteorological data, carrying out reference height identification on the cloudless profile segmentation radar signal and finding out a characteristic height section closest to the molecular scattering signal as a reference height;
the method comprises the steps of identifying reference height, wherein the identification of the reference height comprises signal characteristic decomposition and Rayleigh fit judgment, the signal characteristic decomposition divides a radar signal vertical section into height sections with the same signal characteristics by using a Douglas-Peucker algorithm, then initial detection is carried out on the height sections of the divided signal characteristics, the height sections of the signal characteristics comprise the vertical thickness of the height sections of the characteristics and the positions of the height sections of the characteristics, then Rayleigh fit judgment is carried out on the height sections of the characteristics meeting the limitation of the height range and the minimum thickness, and the height section of the characteristics closest to a molecular scattering signal is found out and used as the reference height; wherein the molecular scattering signal is calculated according to meteorological data, and the calculation formula is as follows:
Figure FDA0003672160490000041
wherein P is m Represents the molecular scattering signal, beta m And alpha m Representing the molecular backscattering coefficient and extinction coefficient; the molecular backscattering coefficient and extinction coefficient can be calculated by the following formulas:
Figure FDA0003672160490000042
Figure FDA0003672160490000043
wherein p (z) is atmospheric pressure in units of hectopascal, T (z) is atmospheric temperature in units of Kelvin, λ 0 Is the wavelength of the excitation light in nanometers; the data inversion module is used for inverting from the continuous cloud-free signal profile by combining the input meteorological data to obtain an aerosol extinction coefficient, a backscattering coefficient and a particle depolarization ratio;
the specific implementation mode is as follows;
for an emission wavelength of λ 0 The radar equation describing the received signal of the elastic scattering lidar of (1) is shown in equation (8):
Figure FDA0003672160490000051
where P (z) is the strength of the radar echo signal at height z, C is the receiving efficiency of the entire system, O (z) is the overlap factor of the system, and β a Is the backscattering coefficient, alpha, of the aerosol a Z' represents an integral variable for the extinction coefficient of the aerosol; the backscattering coefficient and extinction coefficient of aerosol cannot be solved by an equation, and the traditional solution method is to assume that the backscattering coefficient and the extinction coefficient satisfy a linear relationship, as shown in equation (9):
Figure FDA0003672160490000052
wherein S is a As extinction-backscattering ratio of aerosols, also known as laser lasersUp to; based on this relationship, a numerical solution for the backscattering coefficient of the aerosol can be obtained by solving the bernoulli equation as follows:
Figure FDA0003672160490000053
wherein z is ref Is a reference height, obtained by a reference height identification module, and z' represents an integral variable;
after obtaining the backscattering coefficient of the aerosol, the particle depolarization ratio can be obtained by inverting the following formula:
Figure FDA0003672160490000054
Figure FDA0003672160490000055
where K is the ratio of the system efficiencies of the parallel and perpendicular channels, calibrated by the Δ 90 method, P || And P ⊥ Signal strength, delta, for parallel and vertical channels, respectively m Is the molecular depolarization ratio, which is a quantity related only to the bandwidth of the optical filter of the receiving channel and the atmospheric temperature;
and the inversion result output module is used for outputting the inversion result as a scientific data format file with a standard format.
CN202110491517.5A 2021-05-06 2021-05-06 Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition Active CN113325440B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110491517.5A CN113325440B (en) 2021-05-06 2021-05-06 Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110491517.5A CN113325440B (en) 2021-05-06 2021-05-06 Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition

Publications (2)

Publication Number Publication Date
CN113325440A CN113325440A (en) 2021-08-31
CN113325440B true CN113325440B (en) 2022-08-16

Family

ID=77414239

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110491517.5A Active CN113325440B (en) 2021-05-06 2021-05-06 Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition

Country Status (1)

Country Link
CN (1) CN113325440B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117590360B (en) * 2024-01-15 2024-04-16 合肥中科光博量子科技有限公司 Color ratio calibration method for multi-wavelength laser radar
CN118112545A (en) * 2024-04-30 2024-05-31 武汉大学 Atmospheric laser radar reference height extraction method and system based on Rayleigh fitting

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5241315A (en) * 1992-08-13 1993-08-31 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Micro pulse laser radar
US6404494B1 (en) * 1998-12-22 2002-06-11 University Of Washington Measurement of the lidar ratio for atmospheric aerosols using a 180 degree-backscatter nephelometer
CN102628946A (en) * 2012-04-11 2012-08-08 南京信息工程大学 Atmospheric sulfur dioxide and ozone profile Raman-Rayleigh/Lamy multifunctional laser radar measuring device and detection method
CN103630908A (en) * 2013-12-08 2014-03-12 中国科学技术大学 Laser frequency spectrum retrieval method and measurement calibration method for molecular scattering anemometry laser radar
CN107092020A (en) * 2017-04-19 2017-08-25 北京大学 Merge the surface evenness monitoring method of unmanned plane LiDAR and high score image
CN112505651A (en) * 2020-12-23 2021-03-16 北京遥测技术研究所 Automatic processing method for atmospheric detection laser radar

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5241315A (en) * 1992-08-13 1993-08-31 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Micro pulse laser radar
US6404494B1 (en) * 1998-12-22 2002-06-11 University Of Washington Measurement of the lidar ratio for atmospheric aerosols using a 180 degree-backscatter nephelometer
CN102628946A (en) * 2012-04-11 2012-08-08 南京信息工程大学 Atmospheric sulfur dioxide and ozone profile Raman-Rayleigh/Lamy multifunctional laser radar measuring device and detection method
CN103630908A (en) * 2013-12-08 2014-03-12 中国科学技术大学 Laser frequency spectrum retrieval method and measurement calibration method for molecular scattering anemometry laser radar
CN107092020A (en) * 2017-04-19 2017-08-25 北京大学 Merge the surface evenness monitoring method of unmanned plane LiDAR and high score image
CN112505651A (en) * 2020-12-23 2021-03-16 北京遥测技术研究所 Automatic processing method for atmospheric detection laser radar

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Application of the shipborne remote sensing supersite OCEANET for profiling of Arctic aerosols and clouds during Polarstern cruise PS106;HannesJ.Griesche 等;《Atmospheric Measurement Techniques》;20191211;第1-37页 *
基于偏振激光雷达和CALIPSO对武汉上空沙尘气溶胶的观测研究;何芸;《中国优秀博硕士学位论文全文数据库(博士)工程科技Ⅰ辑》;20160715(第7期);第27-33页 *
激光雷达反演大气边界层高度的优化方法;于思琪 等;《光学学报》;20210430;第41卷(第7期);第1-9页 *
钠多普勒激光雷达大气参数反演方法;徐丽 等;《红外与激光工程》;20090228;第38卷(第1期);第140-143页 *

Also Published As

Publication number Publication date
CN113325440A (en) 2021-08-31

Similar Documents

Publication Publication Date Title
CN113325440B (en) Polarization laser radar data inversion method and system based on image recognition and signal characteristic decomposition
CN110082842B (en) Precipitation estimation method and device
CN109581546B (en) Rainfall type identification method based on microwave link attenuation and polarization information
CN111832518B (en) Space-time fusion-based TSA remote sensing image land utilization method
Lee et al. Variability of drop size distributions: Noise and noise filtering in disdrometric data
CN116108008A (en) Decorative material formaldehyde detection data processing method
CN107703554B (en) Multichannel millimeter wave radiometer temperature and humidity profile inversion system and inversion method thereof
CN109657363B (en) Space-time continuous PM2.5 inversion method
CN111191673B (en) Ground surface temperature downscaling method and system
Virtanen et al. Collocation mismatch uncertainties in satellite aerosol retrieval validation
Teinilä et al. A study of the sea-salt chemistry using size-segregated aerosol measurements at coastal Antarctic station Neumayer
CN109255100A (en) A kind of Urban Rain inversion algorithm based on microwave attenuation characteristic response fingerprint recognition
CN108732313A (en) Urban air pollution object concentration intelligence observation system
CN113189014A (en) Ozone concentration estimation method fusing satellite remote sensing and ground monitoring data
CN117390378B (en) Intelligent management method and system for dual-carbon platform data
CN113887324A (en) Fire point detection method based on satellite remote sensing data
Liu et al. Graphics algorithm for deriving atmospheric boundary layer heights from CALIPSO data
He et al. Use of the C‐Band Microwave Link to Distinguish between Rainy and Dry Periods
IL297153A (en) Method for calculating a stream of at least one gas emitted by a source into the atmosphere, measurement method, and associated system and kit
CN113378476A (en) Global 250-meter resolution space-time continuous leaf area index satellite product generation method
CN109374484B (en) Data quality control method of surface laser raindrop spectrometer under strong wind condition
CN116466368A (en) Dust extinction coefficient profile estimation method based on laser radar and satellite data
CN114296103B (en) Airborne high-spectral-resolution laser radar extinction coefficient inversion method
CN112698354B (en) Atmospheric aerosol and cloud identification method and system
CN113625369A (en) Miniaturized intelligent measuring system and method for atmospheric visibility

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant