CN113406594A - Single photon laser fog penetration method based on double-quantity estimation method - Google Patents

Single photon laser fog penetration method based on double-quantity estimation method Download PDF

Info

Publication number
CN113406594A
CN113406594A CN202110607936.0A CN202110607936A CN113406594A CN 113406594 A CN113406594 A CN 113406594A CN 202110607936 A CN202110607936 A CN 202110607936A CN 113406594 A CN113406594 A CN 113406594A
Authority
CN
China
Prior art keywords
target
signal
echo
smoke
photon
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.)
Granted
Application number
CN202110607936.0A
Other languages
Chinese (zh)
Other versions
CN113406594B (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.)
Harbin Institute Of Technology Beijing Industrial Technology Innovation Research Institute Co ltd
Harbin Institute of Technology
Original Assignee
Harbin Institute Of Technology Beijing Industrial Technology Innovation Research Institute Co ltd
Harbin Institute of Technology
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 Harbin Institute Of Technology Beijing Industrial Technology Innovation Research Institute Co ltd, Harbin Institute of Technology filed Critical Harbin Institute Of Technology Beijing Industrial Technology Innovation Research Institute Co ltd
Priority to CN202110607936.0A priority Critical patent/CN113406594B/en
Publication of CN113406594A publication Critical patent/CN113406594A/en
Application granted granted Critical
Publication of CN113406594B publication Critical patent/CN113406594B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/48Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S17/00
    • G01S7/4802Details 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
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/006Theoretical aspects
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G01S17/88Lidar systems specially adapted for specific applications
    • G01S17/89Lidar systems specially adapted for specific applications for mapping or imaging
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information 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)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

The invention discloses a single photon laser fog penetration method based on a double quantity estimation method. Step 1: establishing a fog model based on Gamma distribution; step 2: establishing a target echo signal model; and step 3: detecting total echo photons; and 4, step 4: and (3) extracting a target signal in the total echo photons detected in the step (3) through the fog model in the step (1) and the target echo signal model in the step (2). The method is used for solving the problem that the target signal extraction capability of a full parameter estimation method, a single parameter estimation method and a traditional peak value method is weak.

Description

Single photon laser fog penetration method based on double-quantity estimation method
Technical Field
The invention belongs to the field of Gm-APD laser radar detection; in particular to a single photon laser fog penetration method based on a double quantity estimation method.
Background
With the development of new applications such as automobile laser radars and unmanned aerial vehicles, people are more and more interested in high-resolution imaging of targets in a visual degradation environment, particularly in a dense fog environment. When the laser propagation environment is a highly scattering medium, the main limitation of active optical imaging is the absorption and scattering of the laser light, which results in significant signal attenuation over a relatively short propagation distance. When depth imaging is carried out in dense fog, the number of photons of a target echo is too small, so that a target signal cannot be extracted from the dense fog signal. Therefore, it is necessary to develop an algorithm for separating a target signal from a smoke background, meet the requirements of depth image detection and short-time imaging of a target behind outdoor dense smoke, and improve the weather adaptability of the Gm-APD laser radar.
Conventional Peak method (Peak selection algorithm, PSA): and selecting the peak point of the signal in each pixel point as the distance position of the target.
Single Parameter Estimation Algorithm (SPEA): and reconstructing a depth image by adopting a single parameter estimation construction algorithm (SPEA) according to the actually measured fog attenuation coefficient and the Gamma distribution model of the smog.
Full parameter estimation (APEA): and approximately recognizing the target signal as noise, and directly estimating the photon distribution histogram by adopting maximum likelihood estimation so as to reconstruct the depth image.
Conventional Peak method (Peak selection algorithm, PSA): the peak method can only extract signals of the targets behind the low-concentration smoke. When the smoke concentration is increased, the smoke echo intensity is higher, the target echo intensity is lower, a target peak is almost submerged in the falling edge of the smoke echo signal, and the target signal cannot be found basically.
Single Parameter Estimation Algorithm (SPEA): due to the density, the dynamic property and the non-uniformity of the fog, the true value of the fog is difficult to obtain accurately when the attenuation coefficient of the fog is measured, so that parameter estimation is carried out by using the fog, the error is large, and the scheme is not suitable for outdoor experiments.
Full parameter estimation (APEA): the echo signal comprises a smoke signal and a target signal, and the parameter estimation error is larger by directly performing the parameter estimation method on the echo photon histogram by neglecting the target signal.
Disclosure of Invention
The invention provides a single photon laser fog penetration method based on a double quantity estimation method, which can solve the problem of weak target signal extraction capability in the three methods, improve the estimation precision of the smoke signal of the latter two algorithms, and can better extract weak target signals in an extreme environment.
The invention is realized by the following technical scheme:
a single photon laser fog penetration method based on a double quantity estimation method comprises the following steps:
step 1: establishing a fog model based on Gamma distribution;
step 2: establishing a target echo signal model;
and step 3: detecting total echo photons;
and 4, step 4: and (3) extracting the target signal in the total echo photons detected in the step (3) through the fog model in the step (1) and the target echo signal model in the step (2).
Further, in the step 1, specifically, in the triggering of the Gm-APD detector, when the laser pulse energy is extremely weak, the initial numbers of electrons generated by background noise, dark counting noise and echo photons all obey Poisson distribution; the echo photons include target echo, backscatter, and ambient light 0; at t according to Poisson distribution statistics1To t2The probability of generating q photoelectric events within a detection interval is:
Figure BDA0003094754920000021
wherein Nr (t)1,t2) Is at t1To t2The average number of photons in the detection interval, expressed as:
Figure BDA0003094754920000022
wherein λ (t) is a rate function of the initial electrons in the GM-APD detector; at a detection time t1To t2The probability of not generating photoelectrons is P (0) ═ exp (-Nr (t)1,t2) So that the probability of generating photoelectrons is 1-exp (-Nr (t)1,t2) ); since the Gm-APD detector only produces one detection within the range gate, the probability of producing a detection at the jth slot is:
Figure BDA0003094754920000023
according to a detection probability distribution curve of the photon number output by the Gm-APD detector on a time axis, the average photon number Nr (j) of the j-th time slot reaching the focal plane of the detector is obtained by reverse estimation, and the average photon number Nr (j) is expressed as:
Figure BDA0003094754920000024
because smoke particles have strong scattering and absorption effects on laser light, continuous scattering events exist in photons in the forward transmission process of the laser light; if the initial photon number per scattering event is N0The change of the photon number N (t) after the primary scattering meets the relation with the time;
N(t)=N0exp(-αct) (5)
wherein alpha is the attenuation coefficient of the laser in the fog, and c is the light speed;
the photon probability density function in a primary scattering event is:
Figure BDA0003094754920000031
the detection time tau of each scattered photon is due to the independence of the scattering eventsiAll satisfy the exponential distribution; the photon detection time after multiple scattering is
Figure BDA0003094754920000032
Wherein k is the number of scattering events; gamma distribution is used for solving the problem of time required from the beginning to the occurrence of k random events, and the sum of k independent exponential random variables obeys the Gamma distribution, T-GAMMA (k, beta); the density function of the number of echo photons as a function of time t is thus obtained when the laser is transmitted in the smoke:
Figure BDA0003094754920000033
where r (k) is a Gamma function, k and β are shape and inverse scale parameters; when k is 1, the Gamma distribution is exponential:
GT(t|k=1,β)=βexp(-βt) (8)
let β ═ α c, calculated using equations (4) and (6), and the number of photons reaching the detector focal plane after continuous scattering by smoke is:
Figure BDA0003094754920000034
where r is related to the backscatter coefficient of smoke.
Further, the step 2 specifically includes taking a gaussian function as a fitting model of the target echo signal, where the fitting model is expressed as:
Figure BDA0003094754920000035
wherein, the central position ttargetRelated to the target distance, ttarget2d/c, wherein d is the target distance; introducing a FWHM into the model to characterize the characteristics of the echo signal; FWHM of laser emission pulse after interaction with target surface is taup,τpProportional to σ, obtained by the following equation:
Figure BDA0003094754920000036
Figure BDA0003094754920000037
Figure BDA0003094754920000038
the photon distribution of the target signal photons at time t is obtained from equations (11) to (13):
Figure BDA0003094754920000041
where s is related to the backscatter coefficient of the target.
Further, step 3 specifically includes that the echo signal intensity of the smoke is greater than the target signal intensity due to the fact that the smoke is close to the radar system and the concentration of the smoke is large, detected echo signals are distributed in a double-peak mode, and the echo peak of the target signal exists on the falling edge of the smoke signal; the number of echo photons detected at time t is as follows:
Figure BDA0003094754920000042
further, the step 4 specifically includes the following steps:
step 4.1: estimating the time profile of the echo signal;
step 4.2: determining a fog signal location based on the time profile estimate of step 4.1;
step 4.3: based on the fog signal position of step 4.2, estimating a fog signal;
step 4.4: based on the smoke signal estimation of step 4.3, a target signal estimation is performed.
Further, the step 4.1 is specifically that each pixel point of the photon counting radar records the flight time of the received photon, so that a flight time histogram of the echo photon is obtained, the detection probability in each time bin is obtained according to the number of imaging frames, further, the photon distribution nf (t) reaching the focal plane of the detector is obtained through calculation by a formula (4), and then, filtering is performed by using a gaussian function, so as to obtain the time profile estimation.
Further, the step 4.2 is specifically to estimate peak positions of the smoke and the target echo signal; carrying out target and smoke peak detection on the estimated time contour signal by using continuous wavelet change;
the continuous wavelet transform formula is as follows:
Figure BDA0003094754920000043
wherein CoCsIs the CWT coefficient, s (t) is the signal obtained by Gaussian fitting,
Figure BDA0003094754920000044
is a mother wavelet, p and q are scaling and shifting factors,
Figure BDA0003094754920000045
are scaled and translated wavelets.
Further, the step 4.3 is specifically to use the calculated Nf(t) carrying out GTParameter estimation of (t | k, β); the smoke signal estimation mainly comprises the following steps:
step 4.3.1: obtaining the initial position tau of the smoke echo signal according to CWT transformationaAnd end position τbWill fit the range τaTo taubThe number of inner echo photons is extracted to obtain Nfab);
Step 4.3.2: n from step 4.3.1fab) Normalization processing is carried out, and the mean value E (tau) is obtained by calculationab) Sum variance D (τ)ab). As known from GAMMA distribution function, if t-GAMMA (k, beta), the mean and variance are respectively E (t)=k/β,D(t)=k/β2. Thus, at time τaTo taubbeta-E (T)ab)/D(τab);
Step 4.3.3: according to the beta obtained in the step 4.3.2, utilizing maximum likelihood estimation to obtain an estimated value of a distribution parameter k, and finally obtaining an estimated smoke echo signal distribution GT(t|k,β);
Step 4.3.4: smoke echo signal N from CWTfab) The peak intensity and position of step 4.3.3 is GT(t | k, β) according to Nfab) The maximum values of (A) are unified in scale, the variation scale is r in formula (9), and the smoke echo signal G is obtained after peak value translationTr(t|k,β)。
Further, the step 4.4 is specifically to calculate the total number of echo photons Nf(t) and estimated Smoke Signal GTrSubtracting the two signals to obtain the initial photon number distribution S of the target signalT(t); fitting to obtain FWHM tau by using reflected echo signals after interaction of laser emission pulses and target surfacepFurther, obtaining a system response function R with a Gaussian shape according to the formula (14); extracting depth information of a target from the histogram of each pixel by adopting a cross-correlation method; for each pixel, the cross-correlation coefficient Cr(t) is a histogram S of passing target echoesT(t) is calculated from the system response function R of each pixel, and the expression is as follows:
Figure BDA0003094754920000051
calculating the distance position of the maximum point in the result Cr (t) as a target;
and finally obtaining the depth image of the whole target by estimating the distance of the target in each pixel point.
The invention has the beneficial effects that: (excellence, much good)
Drawings
Fig. 1 is a standard image of an object without the presence of smoke.
Fig. 2 is a histogram distribution of a standard image of an object in which no smoke is present.
FIG. 3 is a diagram of a position distribution of an experimental scene and a target according to the present invention.
FIG. 4 is a diagram of the distribution and signal extraction of the echo signal of a single target pixel under different attenuation coefficients.
Fig. 5 is a comparison graph of target depth image reconstruction results under different attenuation lengths, where (a) is a graph for extracting a target signal by using a peak method for time profile estimation, (b) is a graph for estimating a smoke echo signal according to an attenuation length value obtained by measurement, (c) is a graph for estimating a smoke echo signal according to histogram raw data of photon flight time obtained by measurement, and (d) is a graph for extracting a target signal extracted by the present invention.
FIG. 6 is a comparison of the reconstruction results of the target depth images at different acquisition times, wherein F represents the acquisition frame number (a) in the comparison graphs at 3000, 7500, 12500 and 17500 for the acquisition frame number when estimating the time profile by the peaking method, (b), (c), (d)
Fig. 7 is a graph showing the variation of TR of the reconstructed image with the number of acquired frames when AL is 1.2, (b) a graph showing the variation of RARE of the reconstructed image with the number of acquired frames when AL is 1.2, (c) a graph showing the variation of RARE of the reconstructed image with the number of acquired frames when AL is 3.6, and (d) a graph showing the variation of TR of the reconstructed image with the number of acquired frames when AL is 3.6.
Fig. 8 is a histogram distribution diagram of different algorithm-reconstructed images with a small number of acquisition frames, where (a) the attenuation length AL is 1.2, F is 3000, A, B and Background correspond to the position maps of the target A, B and the Background in the standard image, respectively, (b) the attenuation length AL is 2.7, F is 3000, A, B and Background correspond to the position maps of the target A, B and the Background in the standard image, respectively.
FIG. 9 is a flow chart of the method of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the accompanying drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, are within the scope of the present invention.
A single photon laser fog penetration method based on a double quantity estimation method comprises the following steps:
step 1: establishing a fog model based on Gamma distribution;
step 2: establishing a target echo signal model;
and step 3: detecting total echo photons;
and 4, step 4: and (3) extracting the target signal in the total echo photons detected in the step (3) through the fog model in the step (1) and the target echo signal model in the step (2).
Further, the step 1 is specifically that a pulse laser radar based on a Gm-APD detector measures the distance to a target by a photon time-of-flight method. As an important component of the photon counting radar, a laser emits periodic light beams to a target, a time-to-digital converter (TDC) is started to start timing, and when the TDC stops timing and is triggered when an echo pulse is received, the target distance is obtained through the flight time recorded by the TDC; in the triggering of the Gm-APD detector, when the laser pulse energy is extremely weak, the initial electron numbers generated by background noise, dark counting noise and echo photons are subjected to Poisson distribution; the echo photons include target echo, backscatter, and ambient light 0; at t according to Poisson distribution statistics1To t2The probability of generating q photoelectric events within a detection interval is:
Figure BDA0003094754920000061
wherein Nr (t)1,t2) Is at t1To t2The average number of photons in the detection interval, expressed as:
Figure BDA0003094754920000071
wherein λ (t) is a rate function of the initial electrons in the GM-APD detector; at a detection time t1To t2The probability of not generating photoelectrons is P (0) ═ exp (-Nr (t)1,t2) So that the probability of generating photoelectrons is 1-exp (-Nr (t)1,t2) ); since the Gm-APD detector only produces one detection within the range gate, the probability of producing a detection at the jth slot is:
Figure BDA0003094754920000072
according to a detection probability distribution curve of the photon number output by the Gm-APD detector on a time axis, the average photon number Nr (j) of the j-th time slot reaching the focal plane of the detector is obtained by reverse estimation, and the average photon number Nr (j) is expressed as:
Figure BDA0003094754920000073
because smoke particles have strong scattering and absorption effects on laser light, continuous scattering events exist in photons in the forward transmission process of the laser light; if the initial photon number per scattering event is N0The change of the photon number N (t) after the primary scattering meets the relation with the time;
N(t)=N0exp(-αct) (5)
wherein alpha is the attenuation coefficient of the laser in the fog, and c is the light speed;
the photon probability density function in a primary scattering event is:
Figure BDA0003094754920000074
the detection time tau of each scattered photon is due to the independence of the scattering eventsiAll satisfy the exponential distribution; the photon detection time after multiple scattering is
Figure BDA0003094754920000075
Wherein k is the number of scattering events; gamma distribution is used for solving the problem of time required from the beginning to the occurrence of k random events, and the sum of k independent exponential random variables obeys the Gamma distribution, T-GAMMA (k, beta); the density function of the number of echo photons as a function of time t is thus obtained when the laser is transmitted in the smoke:
Figure BDA0003094754920000076
where r (k) is a Gamma function, k and β are shape and inverse scale parameters; when k is 1, the Gamma distribution is exponential:
GT(t|k=1,β)=βexp(-βt) (8)
let β ═ α c, calculated using equations (4) and (6), and the number of photons reaching the detector focal plane after continuous scattering by smoke is:
Figure BDA0003094754920000077
where r is related to the backscatter coefficient of smoke.
Further, the step 2 is specifically that the received echo signal is a convolution result of the laser emission pulse and the intensities of the target surface reflection signals at different distances, and the emission pulse and the target surface reflection signals of most radar systems approximately satisfy a gaussian distribution; taking a Gaussian function as a fitting model of the target echo signal, wherein the fitting model is expressed as:
Figure BDA0003094754920000081
wherein, the central position ttargetRelated to the target distance, ttarget2d/c, wherein d is the target distance; since the FWHM parameter of the echo signal reflects the surface roughness, some other surfaces areWhat characteristics and the width of the transmit pulse; introducing a FWHM into the model to characterize the characteristics of the echo signal; FWHM of laser emission pulse after interaction with target surface is taup,τpProportional to σ, obtained by the following equation:
Figure BDA0003094754920000082
Figure BDA0003094754920000083
Figure BDA0003094754920000084
the photon distribution of the target signal photons at time t is obtained from equations (11) to (13):
Figure BDA0003094754920000085
where s is related to the backscatter coefficient of the target.
Further, in step 3, specifically, the application scenario of the lidar system is considered that smoke exists between the target and the system, and the smoke concentration is relatively high and is in a nonlinear and unstable state. The echo photons mainly comprise smoke echo photons, target echo photons and noise photons; the noise intensity is extremely weak relative to the smoke intensity, so the model of the invention does not consider noise photons; the smoke is close to the radar system and has high concentration, so that the intensity of an echo signal of the smoke is larger than that of a target signal, the detected echo signals are distributed in a double peak mode, and the echo peak of the target signal exists on the falling edge of the smoke signal; the number of echo photons detected at time t is as follows:
Figure BDA0003094754920000086
further, since each pixel of the area array photon counting radar has time resolution capability, the echo signal in each pixel comprises smoke backscattered light and signal light, and the aim of the radar is to separate the target signal from the smoke background signal. The step 4 specifically comprises the following steps:
step 4.1: estimating the time profile of the echo signal;
step 4.2: determining a fog signal location based on the time profile estimate of step 4.1;
step 4.3: based on the fog signal position of step 4.2, estimating a fog signal;
step 4.4: based on the smoke signal estimation of step 4.3, a target signal estimation is performed.
Further, the step 4.1 is specifically that each pixel point of the photon counting radar records the flight time of the received photon, so that a flight time histogram of the echo photon is obtained, the detection probability in each time bin is obtained according to the number of imaging frames, further, the photon distribution nf (t) reaching the focal plane of the detector is obtained through calculation by a formula (4), and then, filtering is performed by using a gaussian function, so as to obtain the time profile estimation.
Further, the step 4.2 is specifically to cause a phenomenon that wave peaks of the echo signals overlap due to the small distance between the target and the smoke. Estimating peak positions of smoke and a target echo signal in order to separate the target signal from a background smoke signal; wavelet transform has been used as a signal analysis tool with great success in peak detection; continuous Wavelet Transform (CWT) can realize the detection of the quantity and the position of the components of the full-waveform echo signal of the laser radar. Therefore, the estimated time profile signal is subjected to peak detection of the target and the smoke by using continuous wavelet change;
the continuous wavelet transform formula is as follows:
Figure BDA0003094754920000091
wherein CoCsIs the CWT coefficient, s (t) is the signal obtained by Gaussian fitting,
Figure BDA0003094754920000092
is a mother wavelet, p and q are scaling and shifting factors,
Figure BDA0003094754920000093
are scaled and translated wavelets.
Further, in the step 4.3, specifically, the number of echo photons of the laser transmitted through the smoke satisfies the Gamma distribution with time t, so that N obtained by calculation is usedf(t) carrying out GTParameter estimation of (t | k, β); due to Nf(t) contains smoke echo and target echo information, ignores target signal, and directly compares Nf(t) the method of performing parameter estimation results in large parameter estimation errors. To increase GTThe accuracy of the (t | k, β) parameter estimation is calculated only for smoke echo signals at the time of parameter estimation. The smoke signal estimation mainly comprises the following steps:
step 4.3.1: obtaining the initial position tau of the smoke echo signal according to CWT transformationaAnd end position τbWill fit the range τaTo taubThe number of inner echo photons is extracted to obtain Nfab);
Step 4.3.2: n from step 4.3.1fab) Normalization processing is carried out, and the mean value E (tau) is obtained by calculationab) Sum variance D (τ)ab). As can be seen from the GAMMA distribution function, if t GAMMA (k, β), the mean and variance are respectively E (t) k/β, D (t) k/β2. Thus, at time τaTo taubbeta-E (T)ab)/D(τab);
Step 4.3.3: according to the beta obtained in the step 4.3.2, utilizing maximum likelihood estimation to obtain an estimated value of a distribution parameter k, and finally obtaining an estimated smoke echo signal distribution GT(t|k,β);
Step 4.3.4: smoke echo signal N from CWTfab) The peak intensity and position of step 4.3.3 is GT(t | k, β) according to Nfab) The maximum values of (A) are unified in scale, the variation scale is r in formula (9), and the smoke echo signal G is obtained after peak value translationTr(t|k,β)。
Further, the step 4.4 is specifically to calculate the total number of echo photons Nf(t) and estimated Smoke Signal GTrSubtracting the two signals to obtain the initial photon number distribution S of the target signalT(t); due to factors such as unsteadiness of smoke distribution and instability of the system, the target echo photons S of each pixel point are causedT(t) there is a non-uniform distribution. FWHM tau is obtained by fitting reflected echo signals after interaction of laser emission pulses and target surfacepFurther, obtaining a system response function R with a Gaussian shape according to the formula (14); extracting depth information of a target from the histogram of each pixel by adopting a cross-correlation method; for each pixel, the cross-correlation coefficient Cr(t) is a histogram S of passing target echoesT(t) is calculated from the system response function R of each pixel, and the expression is as follows:
Figure BDA0003094754920000101
calculating the distance position of the maximum point in the result Cr (t) as a target;
and finally obtaining the depth image of the whole target by estimating the distance of the target in each pixel point.
In order to test the stability of the radar system, firstly, a target without smoke is imaged before an experiment is started, and whether a histogram of an echo signal meets the experiment requirement or not is analyzed. Imaging results as shown in fig. 1, the contour information of objects a and B is well matched to the actual object situation in fig. 3 and is completely separated from the background. The histogram results of fig. 1 are shown in fig. 2, with photon flight time on the abscissa, 1ns for a time Bin, and the number of photons for that time Bin on the ordinate. The curve mainly comprises three peaks, wherein the difference between the first peak and the second peak is 3 times Bin, namely 3 ns. According to the speed of the laser in the air is about3×108m/s, the distance Δ L between two peaks can be calculated1=3×10-9×3×1080.45 m, 38 times Bin between the first and third peaks, and a distance Δ L between the peaks2About 5.7 m. Calculated result Δ L1、ΔL2The distance from the actual target A, B is equal to the distance between the target a and the background plate. Due to the time jitter of instruments and equipment such as a laser, a detector and the like, the difference between the measured target position and the actual position is 3 times Bin, but the relative position of the target is in accordance with the actual position, so the measurement result of the system meets the experimental requirement. The invention takes the figure 1 as a standard image for judging the subsequent smoke image processing result. Fig. 2 shows histogram distribution of a standard image, three peaks are sequentially represented as a target a, a target B and a background plate, and fig. 3 shows an experimental scene and position distribution of the target.
To test our algorithm's ability to extract weak target signals through smoke, the system imaged target B for 20000 frames. Selecting a pixel point in a target range, respectively obtaining a target echo signal when the attenuation coefficient AL is increased from 1.2 to 3.6, and extracting smoke and a target signal according to the target signal extraction algorithm in section 3, wherein the result is shown in fig. 4, a red dotted line is a fitted smoke signal, and a green solid line is a target signal. Along with the increase of the attenuation coefficient, the number of photons penetrating through smoke is reduced due to the increase of the backscattering capacity of the smoke on laser, the intensity of smoke echo signals is gradually increased, target echo signals are gradually reduced, and the superposed echo signals FT(t) changes from bimodal to unimodal distribution with gradually decreasing echo pulse width. When the attenuation coefficient AL is 3.6, the target signal is completely covered by the smoke signal, the SBR value is-23.2 dB, which is reduced by 16.9dB compared to AL 1.2, but our algorithm can still extract a weak target signal.
FIG. 4 shows the echo signal distribution and signal extraction of a single target pixel under different attenuation coefficients. Histogram of the initial reflected photons calculated. Gaussian fit the curve fit of the signal obtained by Gaussian fitting. Gamma fit background signal estimation. Signal, the difference of the histogram and the Gamma fit, namely the extracted target Signal photon. Target is the Target signal estimation. Model fit-Gamma fitting and sum of target signals.
The system carries out continuous imaging 20000 frames on the target A, B under different attenuation lengths, the corresponding acquisition time is 1s, different algorithms are utilized to carry out target signal extraction, and the reconstructed target depth image is shown in fig. 5. With the attenuation length AL increased from 1.2 to 3.6, the interference capability of smoke on the target signal is gradually enhanced, so that the TR values of the target depth image reconstructed by the four algorithms are smaller and larger, and the TR values of the target depth image reconstructed by the four algorithms are reduced by 99.3%, 86.9%, 48.1% and 65.4%, respectively. Although the percent reduction in TR value for the method of the present invention is 17.3% greater relative to the APEA algorithm, it is reduced by 33.9% and 21.5% relative to PSA and SPEA, respectively. When AL is 3.6, TR obtained by the method of the present invention is respectively improved by 0.2300 and 0.0408, and RARE is respectively reduced by 0.3793 and 0.0231, compared with SPEA and APEA, the restored target contour is more complete, and the extracted target position is more accurate. Therefore, the method has better stability and target signal extraction capability when the attenuation length is changed.
Fig. 5 shows comparison of reconstruction results of target depth images under different attenuation lengths, where different columns show reconstruction results of different algorithms under the same attenuation length, and different rows show reconstruction results of the same algorithm under different attenuation lengths. TR is the target recovery and RARE is the relative average range error. (a) Target signal extraction is performed for time profile estimation using the peak method. (b) The smoke echo signal is estimated from the measured attenuation length value. (c) The smoke echo signals are estimated according to the measured histogram raw data of the photon flight time. (d) Is the target signal extracted by the method of the invention.
When the attenuation length AL is 1.2, continuous multi-frame imaging is performed on the target, the number of acquisition frames F is 3000, 7500, 12500, and 17500, and the acquisition times corresponding to the acquisition frames F are 0.150s, 0.375s, 0.625s, and 0.875s, respectively, target signal extraction is performed by using different algorithms, and a reconstructed target depth image is shown in fig. 6. Although the depth image reconstructed by the method of the present invention has more noise than SPEA and APEA at F-3000, the resulting TR is raised by 0.2361 and 0.2583, respectively, and RARE is reduced by 0.1299 and 0.1137, respectively. The reconstructed target profile of the method of the invention gradually tends to be complete with increasing acquisition time, whereas the SPEA and APEA algorithms are opposite. Therefore, the method of the invention has the best target recovery capability under shorter acquisition time.
In order to study the influence of the acquisition time on the target depth image reconstruction result, the target signal extraction is carried out on the images with the attenuation lengths AL of 1.2 and 3.6 respectively, the imaging frame number is increased from 3000 to 20000, and the variation curves of TR and RARE of the reconstructed image along with the acquisition frame number are shown in FIG. 5. When the smoke concentration is low (AL ═ 1.2), TR of the four algorithm reconstructed images gradually increases and RARE gradually decreases as the number of acquisition frames increases. Compared with the SPEA and APES algorithms, TR of the reconstructed image is averagely improved by 0.3271 and 0.3327, and RARE is averagely reduced by 0.1356 and 0.1146. When the smoke concentration is high (AL ═ 3.6), the TR of the reconstructed image of APEA and the method of the invention increases gradually as the number of acquisition frames increases, whereas the TR of SPEA increases and then decreases and its value approaches zero. Although the RARE values obtained by the method of the present invention are almost equal to the result of the APEA algorithm, the TR is improved by 0.0613 on average. In summary, the method has the maximum TR value and the minimum RARE value, that is, the method has better stability and target signal extraction capability when reconstructing a high scattering medium or a target depth image with very short acquisition time.
And 7, comparing the evaluation indexes of the depth image reconstruction results under different acquisition frame numbers. (a) And (b) curves of TR and RARE of the reconstructed image with the number of frames acquired when AL is 1.2, respectively, (c) and curves of TR and RARE of the reconstructed image with the number of frames acquired when AL is 3.6, respectively.
The function of reconstructing the target depth image is to distinguish the targets in the image according to the distance positions. Therefore, the histogram statistics of the reconstructed target depth image with a small number of acquisition frames is performed, and the result is shown in fig. 8. When the smoke attenuation length AL is 1.2 or AL is 2.7, the target distance value obtained by the APEA and SPEA algorithms is greatly different from the target distance value of the standard image, but the target depth image histogram reconstructed by the method of the present invention is consistent with the standard image and includes three echo peaks. Target search is performed on the histograms of the three algorithms by using a peak method to obtain a target A, B in a reconstructed image and a distance error between a background and a real target, which is specifically shown in table 1. When the attenuation length AL is 1.2 or 2.7, the average error of the method of the present invention is 0.33Bin, which is at least 4Bins higher than the conventional algorithm. The result shows that the method can distinguish multiple targets according to the distance value under the condition of dense smoke or short acquisition time, and has the highest sensitivity and the highest distance measurement precision.
FIG. 8 shows histogram distribution of different algorithms reconstructed images for a small number of acquisition frames, where A, B and Background correspond to the position of the object A, B and the Background in the standard image, respectively. (a) Attenuation length AL 1.2, F3000, (b) attenuation length AL 2.7, F3000.
TABLE 1 comparison of target distance and true distance error for different algorithm reconstructions
Figure BDA0003094754920000121

Claims (9)

1. A single photon laser fog penetration method based on a double quantity estimation method is characterized by comprising the following steps:
step 1: establishing a fog model based on Gamma distribution;
step 2: establishing a target echo signal model;
and step 3: detecting total echo photons;
and 4, step 4: and (3) extracting a target signal in the total echo photons detected in the step (3) through the fog model in the step (1) and the target echo signal model in the step (2).
2. The single photon laser fog penetration method based on the dual-quantity estimation method as claimed in claim 1, characterized in thatSpecifically, in the step 1, in the triggering of the Gm-APD detector, when the laser pulse energy is extremely weak, the initial numbers of electrons generated by background noise, dark count noise and echo photons are all subjected to Poisson distribution; the echo photons include target echo, backscatter, and ambient light 0; at t according to Poisson distribution statistics1To t2The probability of generating q photoelectric events within a detection interval is:
Figure FDA0003094754910000011
wherein Nr (t)1,t2) Is at t1To t2The average number of photons in the detection interval, expressed as:
Figure FDA0003094754910000012
wherein λ (t) is a rate function of the initial electrons in the GM-APD detector; at a detection time t1To t2The probability of not generating photoelectrons is P (0) ═ exp (-Nr (t)1,t2) So that the probability of generating photoelectrons is 1-exp (-Nr (t)1,t2) ); since the Gm-APD detector only produces one detection within the range gate, the probability of producing a detection in the jth slot is:
Figure FDA0003094754910000013
according to a detection probability distribution curve of the photon number output by the Gm-APD detector on a time axis, the average photon number Nr (j) of the j-th time slot reaching the focal plane of the detector is obtained by reverse estimation, and the average photon number Nr (j) is expressed as:
Figure FDA0003094754910000014
due to strong scattering and absorption of laser light by smoke particlesSo that there is a continuous scattering event of photons during forward transmission of the laser; if the initial photon number per scattering event is N0The change of the photon number N (t) after the primary scattering meets the relation with the time;
N(t)=N0exp(-αct) (5)
wherein alpha is the attenuation coefficient of the laser in the fog, and c is the light speed;
the photon probability density function in a primary scattering event is:
Figure FDA0003094754910000021
the detection time tau of each scattered photon is due to the independence of the scattering eventsiAll satisfy the exponential distribution; the photon detection time after multiple scattering is
Figure FDA0003094754910000022
Wherein k is the number of scattering events; gamma distribution is used for solving the problem of time required from the beginning to the occurrence of k random events, and the sum of k independent exponential random variables obeys the Gamma distribution, T-GAMMA (k, beta); the density function of the number of echo photons as a function of time t is thus obtained when the laser is transmitted in the smoke:
Figure FDA0003094754910000023
where r (k) is a Gamma function, k and β are shape and inverse scale parameters; when k is 1, the Gamma distribution is an exponential distribution:
GT(t|k=1,β)=βexp(-βt) (8)
let β ═ α c, calculated using equations (4) and (6), and the number of photons reaching the detector focal plane after continuous scattering by smoke is:
Figure FDA0003094754910000024
where r is related to the backscatter coefficient of smoke.
3. The single photon laser defogging method based on the dual quantity estimation method according to claim 2, wherein the step 2 is specifically to use a gaussian function as a fitting model of the target echo signal, and the fitting model is represented as:
Figure FDA0003094754910000025
wherein, the central position ttargetRelated to the target distance, ttarget2d/c, wherein d is the target distance; introducing a FWHM into the model to characterize the characteristics of the echo signal; FWHM of laser emission pulse after interaction with target surface is taup,τpProportional to σ, obtained by the following equation:
Figure FDA0003094754910000026
Figure FDA0003094754910000027
Figure FDA0003094754910000028
the photon distribution of the target signal photons at time t is obtained from equations (11) to (13):
Figure FDA0003094754910000029
where s is related to the backscatter coefficient of the target.
4. The single photon laser fog penetration method based on the dual quantity estimation method as claimed in claim 3, wherein the step 3 is specifically that the echo signal intensity is greater than the target signal intensity due to the fact that the smoke is close to the radar system and the concentration is greater, the detected echo signals are distributed in a double peak mode, and the target signal echo peak exists at the falling edge of the smoke signal; the number of echo photons detected at time t is as follows:
Figure FDA0003094754910000031
5. the single photon laser fog penetration method based on the dual quantity estimation method as claimed in claim 4, wherein the step 4 specifically comprises the following steps:
step 4.1: estimating the time profile of the echo signal;
step 4.2: determining a fog signal location based on the time profile estimate of step 4.1;
step 4.3: based on the fog signal position of step 4.2, estimating a fog signal;
step 4.4: based on the smoke signal estimation of step 4.3, a target signal estimation is performed.
6. The single photon laser fog penetration method based on the double quantity estimation method as claimed in claim 5, wherein the step 4.1 is specifically that each pixel point of the photon counting radar records the flight time of the received photon, so that a flight time histogram of the echo photon is obtained, the detection probability in each time bin is obtained according to the number of imaging frames, further the photon distribution Nf (t) reaching the focal plane of the detector is obtained through calculation of a formula (4), and then filtering is performed by using a Gaussian function to obtain the time profile estimation.
7. The single photon laser fog penetration method based on the dual quantity estimation method as claimed in claim 5, wherein the step 4.2 is to estimate the peak positions of the smoke and the target echo signal; carrying out target and smoke peak detection on the estimated time contour signal by using continuous wavelet change;
the continuous wavelet transform formula is as follows:
Figure FDA0003094754910000032
wherein CoCsIs the CWT coefficient, s (t) is the signal obtained by Gaussian fitting,
Figure FDA0003094754910000033
is a mother wavelet, p and q are scaling and shifting factors,
Figure FDA0003094754910000034
are scaled and translated wavelets.
8. The single photon laser fog penetration method based on the dual-quantity estimation method as claimed in claim 5, wherein the step 4.3 is to use the calculated Nf(t) carrying out GTParameter estimation of (t | k, β); the smoke signal estimation mainly comprises the following steps:
step 4.3.1: obtaining the initial position tau of the smoke echo signal according to CWT transformationaAnd end position τbWill fit the range τaTo taubThe number of inner echo photons is extracted to obtain Nfab);
Step 4.3.2: n from step 4.3.1fab) Normalization processing is carried out, and the mean value E (tau) is obtained by calculationab) Sum variance D (τ)ab). As can be seen from the GAMMA distribution function, if t GAMMA (k, β), the mean and variance are respectively E (t) k/β, D (t) k/β2. Thus, at time τaTo taubbeta-E (T)ab)/D(τab);
Step 4.3.3: beta obtained according to step 4.3.2Obtaining an estimated value of a distribution parameter k by utilizing maximum likelihood estimation, and finally obtaining the estimated smoke echo signal distribution GT(t|k,β);
Step 4.3.4: smoke echo signal N from CWTfab) The peak intensity and position of step 4.3.3 is GT(t | k, β) according to Nfab) The maximum values of (3) are unified in scale, the variation scale is r in the formula (9), and the smoke echo signal G is obtained after peak value translationTr(t|k,β)。
9. The single photon laser fog penetration method based on the dual-quantity estimation method as claimed in claim 8, wherein the step 4.4 is to calculate the total number of echo photons Nf(t) and the estimated smoke signal GTrSubtracting the two signals to obtain the initial photon number distribution S of the target signalT(t); FWHM tau is obtained by fitting reflected echo signals after interaction of laser emission pulses and target surfacepFurther, obtaining a system response function R with a Gaussian shape according to the formula (14); extracting depth information of a target from the histogram of each pixel by adopting a cross-correlation method; for each pixel, the cross-correlation coefficient Cr(t) is a histogram S of passing target echoesT(t) is calculated from the system response function R of each pixel, and the expression is as follows:
Figure FDA0003094754910000041
calculating the distance position of the maximum point in the result Cr (t) as a target;
and finally obtaining the depth image of the whole target by estimating the distance of the target in each pixel point.
CN202110607936.0A 2021-06-01 2021-06-01 Single photon laser fog penetrating method based on double-quantity estimation method Active CN113406594B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110607936.0A CN113406594B (en) 2021-06-01 2021-06-01 Single photon laser fog penetrating method based on double-quantity estimation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110607936.0A CN113406594B (en) 2021-06-01 2021-06-01 Single photon laser fog penetrating method based on double-quantity estimation method

Publications (2)

Publication Number Publication Date
CN113406594A true CN113406594A (en) 2021-09-17
CN113406594B CN113406594B (en) 2023-06-27

Family

ID=77675670

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110607936.0A Active CN113406594B (en) 2021-06-01 2021-06-01 Single photon laser fog penetrating method based on double-quantity estimation method

Country Status (1)

Country Link
CN (1) CN113406594B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113820726A (en) * 2021-09-30 2021-12-21 中国科学院光电技术研究所 Noise suppression method based on multi-dimensional filtering in non-vision field target detection
CN115097484A (en) * 2022-06-23 2022-09-23 哈尔滨工业大学 double-Gamma estimation-based single photon laser radar fog-penetration imaging method

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08122436A (en) * 1994-10-27 1996-05-17 Kansei Corp Method and equipment for measuring distance
US20120075615A1 (en) * 2009-06-22 2012-03-29 Toyota Motor Europe Nv/Sa Pulsed light optical rangefinder
CN105096272A (en) * 2015-08-19 2015-11-25 常州工学院 De-hazing method based on dual-tree complex wavelet
CN105303532A (en) * 2015-10-21 2016-02-03 北京工业大学 Wavelet domain Retinex image defogging method
CN109791202A (en) * 2016-09-22 2019-05-21 苹果公司 Laser radar with irregular pulse train
CN111079304A (en) * 2019-12-26 2020-04-28 哈尔滨工业大学 Calculation method for farthest detection distance of Gm-APD laser radar
CN111999742A (en) * 2020-07-14 2020-11-27 哈尔滨工业大学 Single-quantity estimation-based Gm-APD laser radar fog-penetrating imaging reconstruction method
US20200380710A1 (en) * 2019-05-31 2020-12-03 University Of Connecticut System and method for optical sensing, visualization, and detection in turbid water using multi-dimensional integral imaging

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08122436A (en) * 1994-10-27 1996-05-17 Kansei Corp Method and equipment for measuring distance
US20120075615A1 (en) * 2009-06-22 2012-03-29 Toyota Motor Europe Nv/Sa Pulsed light optical rangefinder
CN105096272A (en) * 2015-08-19 2015-11-25 常州工学院 De-hazing method based on dual-tree complex wavelet
CN105303532A (en) * 2015-10-21 2016-02-03 北京工业大学 Wavelet domain Retinex image defogging method
CN109791202A (en) * 2016-09-22 2019-05-21 苹果公司 Laser radar with irregular pulse train
US20200380710A1 (en) * 2019-05-31 2020-12-03 University Of Connecticut System and method for optical sensing, visualization, and detection in turbid water using multi-dimensional integral imaging
CN111079304A (en) * 2019-12-26 2020-04-28 哈尔滨工业大学 Calculation method for farthest detection distance of Gm-APD laser radar
CN111999742A (en) * 2020-07-14 2020-11-27 哈尔滨工业大学 Single-quantity estimation-based Gm-APD laser radar fog-penetrating imaging reconstruction method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GUY SATAT等: "Towards photography through realistic fog", 《IEEE》 *
WU QINQIN等: "Continuous wavelet transform and iterative decrement algorithm for the Lidar full-waveform echo decomposition", 《APPLIED OPTICS》 *
高尚: "Gm-APD激光雷达雾天成像重构算法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113820726A (en) * 2021-09-30 2021-12-21 中国科学院光电技术研究所 Noise suppression method based on multi-dimensional filtering in non-vision field target detection
CN113820726B (en) * 2021-09-30 2023-06-13 中国科学院光电技术研究所 Noise suppression method based on multidimensional filtering in non-visual field target detection
CN115097484A (en) * 2022-06-23 2022-09-23 哈尔滨工业大学 double-Gamma estimation-based single photon laser radar fog-penetration imaging method

Also Published As

Publication number Publication date
CN113406594B (en) 2023-06-27

Similar Documents

Publication Publication Date Title
CN101509972B (en) Wideband radar detecting method for correcting correlation matrix based on high resolution target distance image
CN110554404B (en) Gm-APD array laser radar imaging method and system under strong background noise
CN108304781B (en) Area array Geiger APD laser imaging radar image preprocessing method
CN113406594A (en) Single photon laser fog penetration method based on double-quantity estimation method
CN104914446A (en) Three-dimensional distance image time domain real-time denoising method based on photon counting
CN110031821B (en) Vehicle-mounted obstacle avoidance laser radar waveform extraction method, laser radar and medium
CN114089366A (en) Water body optical parameter inversion method of satellite-borne single photon laser radar
CN105954733A (en) Time-domain filtering method based on photon flight time correlation
CN111999742B (en) Single-quantity estimation-based Gm-APD laser radar fog-penetrating imaging reconstruction method
CN111323765A (en) Satellite-borne photon counting laser radar echo signal processing and target extraction method
Wu et al. Improvement of detection performance on single photon lidar by EMD-based denoising method
CN114636991A (en) Time of flight calculation using bin delta estimation
Zhang et al. Dual-parameter estimation algorithm for Gm-APD Lidar depth imaging through smoke
CN113253240B (en) Space target identification method based on photon detection, storage medium and system
Zhang et al. Three-dimensional imaging of ships in the foggy environment using a single-photon detector array
US20210302554A1 (en) System and method for lidar defogging
Jutzi et al. Measuring and processing the waveform of laser pulses
US20230194666A1 (en) Object Reflectivity Estimation in a LIDAR System
CN115616608B (en) Single photon three-dimensional imaging distance super-resolution method and system
CN108008374B (en) Sea surface large target detection method based on energy median
CN115097484A (en) double-Gamma estimation-based single photon laser radar fog-penetration imaging method
DE10001015C2 (en) Method for measuring the distance of objects, atmospheric particles and the like by means of lidar or laser radar signals
Castorena et al. Compressive sampling of lidar: Full-waveforms as signals of finite rate of innovation
CN112767284B (en) Laser three-dimensional imaging cloud and fog back scattering filtering method and system based on photon counting entropy
CN116930125B (en) Method for measuring attenuation coefficient of backward scattering full-gating imaging water body

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