Satellite-borne single-photon laser radar sea surface ranging error correction method

Abstract

The invention discloses a satellite-borne single-photon laser radar sea surface ranging error correction method. The method comprises the following steps: (1) correcting a single-photon laser radar internal clock error, a transmission pulse time delay and a pulse detection error; (2) accumulating multiple photon events through a static target to generate a photon histogram, and extracting time gravity center compensation system drift; (3) converting laser elevation data and potential height, performing interpolation in combination with NCEP related data, and resolving surface atmospheric pressure at a laser foot point to complete NCEP atmospheric pressure data correction; (4) adopting a CFA2.2 model to construct an atmospheric delay mapping function suitable for different laser pointing nadir angles, calculating atmospheric dry term retardation and wet term retardation, and correcting a laser measurement distance; (5) completing sea surface point cloud data coarse denoising according to sea surface point cloud density distribution characteristics, and extracting a sea surface fluctuation function; and (6) constructing a dynamic parameter optimized DBSCAN denoising algorithm model, completing fine denoising of the point cloud data, and realizing extraction of the effective distance measurement value of the sea surface point cloud.

Description

Technical Field

[0001] The invention relates to a method for correcting sea surface ranging errors of a satellite-borne single-photon laser radar, and belongs to the field of laser optics. Background Art

[0002] The photon counting detection system has unique advantages in sea surface detection due to its micro-pulse, multi-beam, high sensitivity and high repetition rate. The sea surface ranging error of satellite-borne single-photon lidar includes the error of the transceiver link, the photon detection probability error and the atmospheric refraction error. At present, the construction of the atmospheric refraction error mapping function of the satellite-borne laser altimetry system does not take into account the influence of atmospheric refraction. As the laser pointing angle increases, the delay error of this method increases exponentially, and the influence of atmospheric refraction is difficult to ignore. The photon detection probability error is usually calculated using laboratory static targets. In actual on-orbit dynamic ranging, the introduction of various noises will affect the solution of the laser ranging value, and error correction needs to be performed in combination with the characteristics of the sea surface point cloud.

[0003] At present, the main means to improve the ranging accuracy of satellite-borne single-photon lidar focuses on the refined processing of point cloud data, denoising single photons and clustering of signal and noise photons. Commonly used unsupervised classification methods include DBSCAN algorithm, LOF algorithm, OPTICS algorithm, PSO algorithm, etc., and supervised classification methods include random forest method, support vector machine, conditional random field, etc. The unsupervised classification method clusters according to the point cloud density, and the key lies in determining the clustering parameters related to the target characteristics. Affected by factors such as sea breeze and atmosphere, there are a large number of noise photons in the single-photon data of the sea surface, and the point cloud signal of the sea surface shows large fluctuations, and the point cloud density of signal photons varies along the orbital direction. Extracting effective point cloud information on the sea surface and correcting errors are the basis for laser products to exert their application efficiency. Summary of the invention

[0004] The technical problem solved by the present application is: to overcome the deficiencies of the prior art, to provide a method for correcting the sea surface ranging error of a satellite-borne single-photon lidar, and to realize the ranging error correction for the original point cloud data of the sea surface; to overcome the problem in the prior art that "when performing dynamic ranging on orbit, the sea surface is affected by factors such as sea breeze and atmosphere, the sea surface point cloud signal presents large fluctuations, and the point cloud density varies to a certain extent along the orbit direction, and the introduction of various noises produces ranging errors."

[0005] This patent designs a method for correcting the sea surface ranging error of satellite-borne single-photon lidar, which comprehensively considers the error of the transceiver link, the error of the photon detection probability, the atmospheric refraction error and the sea surface point cloud ranging characteristics. An atmospheric refraction mapping function model is constructed to reduce the delay calculation error when the laser pointing angle is large; a DBSCAN denoising algorithm with dynamic parameter optimization is designed to adapt to the calculation of sea surface ranging values.

[0006] The technical solutions provided by this application are as follows:

[0007] A method for correcting sea surface ranging errors of a spaceborne single-photon laser radar comprises the following steps:

[0008] Step 1: The single photon flight time before correction is t. According to the internal clock error Δt1, the transmission pulse delay Δt2 and the pulse detection error Δt3 of the single photon lidar, the single photon flight time after the transceiver link error correction is It can be expressed as

[0009] Step 2: Use the average value t of the single-photon detection results accumulated from 200 consecutive emissions of the same target mean With the true value t tar Calculate the system drift error Δt4 = t mean -t tar , combined with the calculation results in step 1, the photon flight time of a single photon after instrument correction is

[0010] Step 3: Interpolate the 1°×1° latitude and longitude grid sampling data provided by NCEP. The conversion relationship between the potential height H and the laser footprint elevation h is: Where φ is the geographic latitude, R e =6371009m is the radius of the earth, g0=9.80665m / s 2 is the average gravitational acceleration in the atmosphere, and the constant term g eq =9.7803267715m / s 2 , k e =0.001931851353, e 2 =0.00669438002290. In the atmosphere NCEP, the upper and lower potential layers of the two standard atmospheric pressure layers adjacent to the potential height H are taken as T0, T1 and T2, respectively. H0, the temperature, relative humidity and potential height of the lower layer are recorded as T1, H1, temperature and relative humidity vary linearly with the potential height between different standard atmospheric pressure layers, that is Among them, T, R h are the temperature and relative humidity of the standard pressure layer respectively;

[0011] Step 4: Calculate the water vapor partial pressure P w (R h ,T)=R h ·P s (T), where T, R hare the temperature and relative humidity of the standard pressure layer after interpolation, P s is the saturated water vapor pressure, which is calculated as follows:

[0012]

[0013] In the formula, the constant P b =1000hPa, E s (δ) is the Chebyshev polynomial, coefficient a s ={2794.027,1430.604,-18.234,7.674,-0.022,0.263,0.146,0.055,0.033,0.015,0.013}, where s = 0,1,2,…,10, T is the temperature of the standard pressure layer after interpolation, and the temperature constant T max =648K,T min =273K.

[0014] Step 5: Write the atmospheric pressure hydrostatic equation in the form of potential height to construct the atmospheric pressure P and potential height H gt The relationship equation is as follows:

[0015]

[0016] in, They represent the inverse of the compressibility of dry air and water vapor, respectively. R = 8314.51 J / (kmol*K) is the gas constant, M d 、M w are the molecular weight of dry air and water vapor, respectively, d =28.9644 g / mol, M w =18.0152 g / mol, T is the temperature of the interpolated standard pressure layer, P w is the partial pressure of water vapor. The fourth-order Runge-Kutta formula is used to solve the above equation to obtain the surface atmospheric pressure P at the laser foot point. surf .

[0017] Step 6: Calculate the atmospheric delay mapping function and use the CFA2.2 model (the mapping function model for calculating atmospheric delay values) to determine the dry and wet term mapping functions m d (ε), m w (ε) is as follows:

[0018]

[0019] Wherein, the subscript j = d, w is used to distinguish between dry and wet terms, ε is the laser beam altitude angle, when the laser points to the nadir angle θ = 90-ε < 10°, a j 、b j 、cj The value of is much smaller than 1, and the above formula can be simplified to As the laser pointing nadir angle θ increases, the parameter a j =1.185×10 ；3 [1+6.071×10 ；5 (P surf -1000)-1.471×10 ；4 P w +3.072×10 ；3 (T c -20)+1.956×10 ；2 (ξ+6.5)-5.645×10 ；3 (A j -11.231)],b j =1.144×10 ；3 [1+1.164×10 ；5 (P surf -1000)-2.795×10 ；4 P w +3.109×10 ；3 (T c -20)+3.038×10 ；2 (ξ+6.5)-1.217×10 ；2 (A j -11.231)],c j =-0.009, where ξ is the temperature attenuation rate with altitude, here ξ = -6.5℃ / km, A j (j=d,w) is the height parameter of dry and wet items, A d =40.136+0.14872·T c , A w =11, T c is the temperature in °C, with T c =T-273.16, P w is the water vapor partial pressure, P surf is the surface atmospheric pressure.

[0020] The single-photon lidar of this patent emits multi-beam lasers, and there are beams with a nadir angle greater than 10°; using the CFA2.2 model, the accuracy of the atmospheric delay error ΔL is improved to the millimeter level, and the larger the nadir angle, the more obvious the improvement effect.

[0021] Step 7: Decompose the atmospheric delay into dry delay ΔL d and wet term delay ΔL w , the atmospheric delay error is Among them, ρ wa =1.0×103 kg / m 3 is the density of water, P w , P surf 、m d (ε), m w (ε) is obtained from the above calculation, M d 、M w , R, g0 are the known constants mentioned above, k1 and k2 are functions of the laser center wavelength λ, calculated according to the empirical formula given by Owens:

[0022]

[0023] When λ=1064nm, k1=0.80277K / Pa, k2=0.66388K / Pa; when λ=532nm, k1=0.77493K / Pa, k2=64873K / Pa;

[0024] According to the photon flight time corrected by the instrument Get the measured distance L after the laser is refracted by the atmosphere, c is the speed of light;

[0025] The actual laser straight line measurement distance is Measure distance based on actual laser straight line Obtain photon flight time after atmospheric delay error correction

[0026] Step 8: Convert the input sea surface laser point cloud data into single photon flight time after atmospheric delay error correction And the actual laser straight line measurement distance Perform photon flight time and laser measurement distance correction, and the processed point cloud data is recorded as point r (x r ,y r ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction, respectively, and r takes the value of 1, 2…T r , represents the point cloud data sequence, the total number of point clouds is T r . Press point r middle y r Take the value and mark the vertical layer as z p , where z is the vertical distance between the centers of each layer, p is the subscript of each layer, p∈[1,M], there are M height layers in total, and the number of photons in each height layer is recorded as num p ;

[0027] Step 9: Set the coarse denoising threshold Among them, e∈(1,M / 3) and e is a positive integer, i is the subscript of each altitude layer, traverse each altitude layer, if the corresponding altitude layer num p >th is established, corresponding to the altitude layer z p The data in the middle is the photon point cloud data of the sea surface after rough denoising, denoted as point rc (x rc ,y rc ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction respectively, rc takes the value of 1, 2…Tm, representing the point cloud data sequence, and the total number of point clouds after coarse denoising is Tm;

[0028] Step 10: Based on the photon point cloud data point after rough denoising in step 9 rc Medium rc The value is taken, and the segment along the track is marked as w l , where w is the distance along the track from the center of each segment, l is the subscript of each segment, l∈[1,V], there are V segments along the track, and the number of photons extracted from each segment is recorded as lum l , the median extraction result of the photon height of each segment is recorded as zh l If the segment w l lum l = 0, then the median photon height of the segment is Where Vn≤V, Vn is the number of segments with non-zero photon counts in the along-track segment, li is the segment subscript, and the average number of photons in the non-zero segments

[0029]

[0030] Step 11: Separate the segments l By subtracting adjacent values, we can get the median change function of photon height along the track. l∈[1,V-1], for the difference diff(w l ) is smoothed to obtain the sea surface fluctuation function slope(w l ).

[0031] Step 12: Set the neighborhood radius eps and the core point threshold minpts. For the continuous point cloud data of the sea surface, the neighborhood is set to an ellipse with the major and minor axes e respectively. a =5·eps,e b =0.4·eps, the direction of the major axis of the ellipse neighborhood is based on the slope(w l ) is determined, and the core point threshold minpts is determined according to the number of photons lum in each segment in step 10 l , the major axis and minor axis of the neighborhood are adjusted.

[0032] Step 13: Substitute the above parameter values. According to the sea surface fluctuation function slope(w l ), use the DBSCAN algorithm with dynamic parameter optimization to perform fine denoising on the photon point cloud data point rc (x rc ,y rc ). The denoised point cloud data is denoted as point rf (x rf ,y rf ), where x and y respectively represent the range information of the point cloud point along the track and in the vertical direction. The rf value is 1, 2... Tn, representing the point cloud data sequence. The total number of points in the fine-denoised point cloud is Tn. (x rf ,y rf ) is the screening result of the effective measurement distance of the point cloud.

[0033] For the above method for correcting the sea surface ranging error of a spaceborne single-photon lidar, in Step 1, the calculation method of the internal clock error Δt1 is as follows:

[0034] According to the average delay Δt f measured by each unit of the logic gate circuit delay carry chain in the photon counting electronic acquisition card, the standard deviation , and the number of units NE in the delay chain, calculate the internal clock error as Δt1 = n·Δt f , where n < NE and n is an integer.

[0035] For the above method for correcting the sea surface ranging error of a spaceborne single-photon lidar, in Step 1, the calculation method of the emission pulse delay Δt2 is as follows:

[0036] Δt2 = Δt o + Δt c

[0037] where Δt o is the delay from the internally generated emission signal to the laser receiving the emission signal, and Δt c is the delay from the laser receiving the emission signal to the pulse emission.

[0038] For the above method for correcting the sea surface ranging error of a spaceborne single-photon lidar, in Step 1, the calculation method of the pulse detection error Δt3 is as follows:

[0039] Set the emission pulse detection thresholds at two levels, and record the moments when the rising edge and falling edge of the emission pulse reach the two detection thresholds respectively Calculate the difference between the emission pulse peak values obtained by using the Gaussian fitting method and a single detection threshold at the four moment points:

[0040] Δt3 = t pos1 - t pos2=ma internal

[0041] Among them, t pos1 ,t pos2 They are the moments corresponding to the peak values ​​of the emission pulses obtained by Gaussian fitting method and single threshold, respectively.

[0042] In the above-mentioned method for correcting the sea surface ranging error of a satellite-borne single-photon laser radar, in step 5, the fourth-order Runge-Kutta formula is used to obtain the surface atmospheric pressure P at the laser foot point. surf Here’s how:

[0043]

[0044] Set the initial value P0 = P w , the number of iterations is N kr , Ls is the corresponding integration step, kr1, kr2, kr3, kr4 are Runge-Kutta formula parameters, H gti is the potential height at the i-th iteration, P i The potential height is H gti Atmospheric pressure at that time, take

[0045] In the above-mentioned method for correcting sea surface ranging errors of a satellite-borne single-photon laser radar, in step eight, the sea surface laser point cloud data correction method is as follows:

[0046] point or (x or ,y or ) represents the information of the original single photon event record, and defines the earliest received single photon in the extracted point cloud data as point1, point or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or and Tt1 are point or and the photon return time of point 1; vt is the flying speed of the satellite relative to the ground, y or By point or Satellite attitude and orbit at time, photon flight time t or and DEM are obtained by deriving and calculating the rigorous geometric model;

[0047] point or (x or ,y or ) Photon return time after instrument calibration Photon flight time after instrument correction Single photon point orThe flight time is t or ; After atmospheric delay error correction, single photon point or The actual laser straight line measurement distance is ΔL or For single photon point or Atmospheric delay error; Photon flight time after atmospheric delay error correction Based on satellite attitude and orbit, photon return time Photon flight time The DEM is updated by rigorous geometric model calculation or ,y or , get the corrected point cloud data and record point r (x r ,y r ).

[0048] In the above-mentioned method for correcting sea surface ranging errors of a satellite-borne single-photon laser radar, in step 12, the calculation method of the core point threshold minpts is as follows:

[0049]

[0050] Where Δw = w l:1 -w l , when the segmentation is uniform along the track direction, Δw does not change with the value of l, minpts(w l ) is the value for different track segments w l Set the core point threshold, The average number of photons in the non-zero segments.

[0051] In the above-mentioned method for correcting sea surface ranging errors of a space-borne single-photon laser radar, in step thirteen, the DBSCAN algorithm for dynamic parameter optimization is as follows:

[0052] Point cloud after rough denoising rc (x rc ,y rc ) is the center, rc takes the value of 1, 2…Tm, and the point cloud point is determined. rt (x rt ,y rt ),(rt≠rc,and rt∈[1,Tm]) whether the condition is satisfied:

[0053]

[0054] Among them, w l Point cloud point rc The track segment to which it belongs, x′ rt and y′ rtare the distance calculation process quantities along the major axis and minor axis of the ellipse neighborhood, respectively. If f(rc,rt)<1, it is considered that point rt At point rc Otherwise, point rt Not at point rc The statistics are based on point rc (x rc ,y rc ) as the center, the point that meets the above conditions rt The number of point clouds rum(rc), if rum(rc)>minpts(w l ), point rc is valid point cloud data; otherwise, point rc Invalid point cloud data will be removed during fine denoising.

[0055] Compared with the prior art, the present invention has the following beneficial effects:

[0056] (1) The present invention analyzes the single-photon lidar ranging transceiver link error, photon detection probability error and atmospheric refraction error, and performs error correction analysis based on the point cloud characteristics of the sea surface application scenario to achieve on-orbit dynamic ranging error correction.

[0057] (2) The present invention takes into account the influence of atmospheric refraction, uses the CFA2.2 model to construct dry and wet term mapping functions, introduces the calculation results of atmospheric parameters related to NCEP data to determine the mapping function parameters, and improves the calculation accuracy of the atmospheric refraction delay error when the laser pointing angle is large.

[0058] (3) The present invention extracts the density parameters in the along-track direction and combines them with the sea surface fluctuation function to design a DBSCAN denoising algorithm with dynamic parameter optimization, thereby realizing adaptive adjustment of the neighborhood direction and discrimination threshold. This effectively avoids the problem of uneven distribution of along-track point cloud density caused by sea breeze and fog and the loss of some signal point clouds in the extraction results of the point cloud density denoising algorithm due to sea surface fluctuations, and is more suitable for solving sea surface ranging values. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 A schematic diagram of ranging error analysis of a space-borne single-photon laser radar facing the sea surface provided in an embodiment of the present invention;

[0060] Figure 2 A schematic flow chart of a method for correcting sea surface ranging errors of a satellite-borne single-photon lidar provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0061] To make the objectives, technical solutions, and advantages of the present invention clearer, the following will further describe in detail the disclosed embodiments of the present invention with reference to the accompanying drawings.

[0062] The single-photon lidar emits a laser beam. After being reflected by the sea surface, the photons received by the single-photon lidar form photon events. Based on multiple photon events, the photon point cloud data of the sea surface is obtained.

[0063] Figure 1 The figure shows a schematic diagram of the ranging error analysis of the spaceborne single-photon lidar for the sea surface provided by the embodiment of the present invention. The effective ranging value of the sea surface is extracted in combination with the detection link and the application scenario. The detection link includes the transceiver link error, the photon detection probability model error, and the optical path difference caused by atmospheric refraction. Among them, the transceiver link error includes the internal clock error of the logic gate circuit delay, the delay from the internal clock to the emission pulse, and the detection error of the emission pulse moment. The photon event detection probability model analyzes the static target system drift error. The optical path difference caused by atmospheric refraction includes the analysis of the dry and wet term refractive indices and the mapping function. The application scenario for the sea surface includes rough denoising of the dynamic target point cloud and fine denoising of the dynamic target point cloud.

[0064] Figure 2 The figure shows a schematic flowchart of the ranging error correction method of the spaceborne single-photon lidar for the sea surface provided by the embodiment of the present invention, which is used to correct the photon point cloud data of the sea surface. The method includes the following steps:

[0065] Step 1: The single-photon flight time before correction is t. According to the internal clock error Δt1, the emission pulse delay Δt2, and the pulse detection error Δt3 of the single-photon lidar, the single-photon flight time after being corrected by the transceiver link error can be expressed as The specific error calculation process is as follows:

[0066] (1a). According to the average delay Δt measured by each unit of the carry chain of the logic gate circuit delay in the photon counting electronic acquisition card f , the standard deviation and the number of units NE in the delay chain, calculate the internal clock error as Δt1 = n·Δt f , where n < NE and is an integer;

[0067] (1b). The emission pulse delay Δt2 = Δt o +Δt c , where Δt o is the delay from the internally generated emission signal to the laser receiving the emission signal, and Δt c is the delay from the laser receiving the emission signal to the pulse emission;

[0068] (1c) Set two high and low detection thresholds for the transmitted pulse, and record the time when the rising and falling edges of the transmitted pulse reach the two detection thresholds respectively. Calculate the difference between the peak values ​​of the transmitted pulses obtained at four time points using the Gaussian fitting method and a single detection threshold:

[0069] Δt3=t pos1 -t pos2 =ma internal

[0070] Among them, t pos1 ,t pos2 They are the moments corresponding to the peak values ​​of the emission pulses obtained by Gaussian fitting method and single threshold, respectively.

[0071] Step 2: Use the average value t of the single-photon detection results accumulated from 200 consecutive emissions of the same target mean With the true value t tar Calculate the system drift error Δt4 = t mean -t tar , combined with the calculation results in step 1, the photon flight time of a single photon after instrument correction is

[0072] Step 3: Interpolate the 1°×1° latitude and longitude grid sampling data provided by NCEP. The conversion relationship between the potential height H and the laser footprint elevation h is: Where φ is the geographic latitude, R e =6371009m is the radius of the earth, g0=9.80665m / s 2 is the average gravitational acceleration in the atmosphere, and the constant term g eq =9.7803267715m / s 2 , k e =0.001931851353, e 2 =0.00669438002290. The upper and lower potential layers of the two standard atmospheric pressure layers adjacent to H in the atmosphere NCEP are taken as T0, H0, the temperature, relative humidity and potential height of the lower layer are recorded as T1, H1, temperature and relative humidity vary linearly with the potential height between different standard atmospheric pressure layers, that is

[0073] Step 4: Calculate the water vapor partial pressure P w (R h ,T)=R h ·P s (T), where T, R hare the temperature and relative humidity of the standard pressure layer after interpolation, P s is the saturated water vapor pressure, which is calculated as follows:

[0074]

[0075] In the formula, the constant P b =1000hPa, E s (δ) is the Chebyshev polynomial, coefficient a s ={2794.027,1430.604,-18.234,7.674,-0.022,0.263,0.146,0.055,0.033,0.015,0.013}, where s = 0,1,2,…,10, T is the temperature of the standard pressure layer after interpolation, and the temperature constant T max =648K,T min =273K.

[0076] Step 5: Write the atmospheric pressure hydrostatic equation in the form of potential height to construct the atmospheric pressure P and potential height H gt The relationship equation is as follows:

[0077]

[0078] in, They represent the inverse of the compressibility of dry air and water vapor, respectively. R = 8314.51 J / (kmol*K) is the gas constant, M d 、M w are the molecular weight of dry air and water vapor, respectively, d =28.9644 g / mol, M w =18.0152 g / mol, T is the temperature of the interpolated standard pressure layer, P w is the water vapor partial pressure. The fourth-order Runge-Kutta formula is used to calculate the surface atmospheric pressure P at the laser foot point surf Here’s how:

[0079]

[0080] Set the initial value P0 = P w , the number of iterations is N kr , here we take N kr = 10000 times, Ls is the corresponding integration step, kr1, kr2, kr3, kr4 are Runge-Kutta formula parameters, H gti is the potential height at the i-th iteration, P i The potential height is H gti Atmospheric pressure at that time, take

[0081] Step 6: Calculate the atmospheric delay mapping function and use the CFA2.2 model to determine the dry and wet term mapping functions m d (ε), m w (ε) is as follows:

[0082]

[0083] Wherein, the subscript j = d, w is used to distinguish between dry and wet terms, ε is the laser beam altitude angle, when the laser points to the nadir angle θ = 90-ε < 10°, a j 、b j 、c j The value of is much smaller than 1, and the above formula can be simplified to As the laser pointing nadir angle θ increases, the parameter a j =1.185×10 ；3 [1+6.071×10 ；5 (P surf -1000)-1.471×10 ；4 P w +3.072×10 ；3 (T c -20)+1.956×10 ；2 (ξ+6.5)-5.645×10 ；3 (A j -11.231)],b j =1.144×10 ；3 [1+1.164×10 ；5 (P surf -1000)-2.795×10 ；4 P w +3.109×10 ；3 (T c -20)+3.038×10 ；2 (ξ+6.5)-1.217×10 ；2 (A j -11.231)],c j =-0.009, where ξ is the temperature attenuation rate with altitude, here ξ = -6.5℃ / km, A j (j=d,w) is the height parameter of dry and wet items, A d =40.136+0.14872·T c , A w =11, T c is the temperature in °C, with T c =T-273.16, P w is the water vapor partial pressure, P surf is the surface atmospheric pressure.

[0084] Step 7: Decompose the atmospheric delay into dry delay ΔL d and wet term delay ΔL w , the atmospheric delay error is Among them, ρ wa =1.0×10 3 kg / m 3 is the density of water, P w , P surf 、m d (ε), m w (ε) is obtained from the above calculation, M d 、M w , R, g0 are the known constants mentioned above, k1 and k2 are functions of the laser center wavelength λ, calculated according to the empirical formula given by Owens:

[0085]

[0086] When λ=1064nm, k1=0.80277K / Pa, k2=0.66388K / Pa; when λ=532nm, k1=0.77493K / Pa, k2=64873K / Pa, according to the photon flight time after instrument correction Get the measured distance after the laser is refracted by the atmosphere c is the speed of light, and the actual laser straight line measurement distance is Measure distance based on actual laser straight line Obtain photon flight time after atmospheric delay error correction

[0087] Step 8: Convert the input sea surface laser point cloud data into single photon flight time after atmospheric delay error correction And the actual laser straight line measurement distance Perform photon flight time and laser distance measurement correction. point or (x or ,y or ) represents the information of the original single photon event record, and defines the earliest received single photon in the extracted point cloud data as point1, point or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or and Tt1 are point or and the photon return time of point 1; vt is the flying speed of the satellite relative to the ground, y or By point or Satellite attitude and orbit at time, photon flight time tor and DEM are obtained by deriving and calculating the rigorous geometric model;

[0088] point or (x or ,y or ) Photon return time after instrument calibration Photon flight time after instrument correction Single photon point or The flight time is t or ; After atmospheric delay error correction, single photon point or The actual laser straight line measurement distance is ΔL or For single photon point or Atmospheric delay error; Photon flight time after atmospheric delay error correction Based on satellite attitude and orbit, photon return time Photon flight time The DEM is updated by rigorous geometric model calculation or ,y or , get the corrected point cloud data and record point r (x r ,y r ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction, respectively, and r takes the value of 1, 2…T r , represents the point cloud data sequence, the total number of point clouds is T r . Press point r middle y r Take the value and mark the vertical layer as z p , where z is the vertical distance between the centers of each layer, p is the subscript of each layer, p∈[1,M], there are M height layers in total, and the number of photons in each height layer is recorded as num p ;

[0089] Step 9: Set the coarse denoising threshold Among them, e∈(1,M / 3) and e is a positive integer, i is the subscript of each altitude layer, traverse each altitude layer, if the corresponding altitude layer num p >th is established, corresponding to the altitude layer z p The data in the middle is the photon point cloud data of the sea surface after rough denoising, denoted as point rc (x rc ,y rc ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction respectively, rc takes the value of 1, 2…Tm, representing the point cloud data sequence, and the total number of point clouds after coarse denoising is Tm;

[0090] Step 10: Based on the photon point cloud data point after rough denoising in step 9 rc Medium rc The value is taken, and the segment along the track is marked as w l , where w is the distance along the track from the center of each segment, l is the subscript of each segment, l∈[1,V], there are V segments along the track, and the number of photons extracted from each segment is recorded as lum l , the median extraction result of the photon height of each segment is recorded as zh l If the segment w l lum l = 0, then the median photon height of the segment is Where Vn≤V, Vn is the number of segments with non-zero photon counts in the along-track segment, li is the segment subscript, and the average number of photons in the non-zero segments

[0091]

[0092] Step 11: Separate the segments l By subtracting adjacent values, we can get the median change function of photon height along the track. l∈[1,V-1], for the difference diff(w l ) is smoothed to obtain the sea surface fluctuation function slope(w l ).

[0093] Step 12: Set the neighborhood radius eps and the core point threshold minpts. For the continuous point cloud data of the sea surface, the neighborhood is set to an ellipse with the major and minor axes e respectively. a =5·eps,e b =0.4·eps, the direction of the major axis of the ellipse neighborhood is based on the slope(w l ) is determined according to the number of photons lum in each segment in step 10 l The calculation method of core point threshold minpts is as follows:

[0094]

[0095] Where Δw = w l:1 -w l , when the segmentation is uniform along the track direction, Δw does not change with the value of l, minpts(w l ) is the value for different track segments w l Set the core point threshold, The average number of photons in the non-zero segments.

[0096] Step 13: Substitute the above parameter values ​​and use the DBSCAN algorithm with dynamic parameter optimization to denoise the photon point cloud data point rc (xrc ,y rc ) for fine denoising, and the point cloud point after coarse denoising rc (x rc ,y rc ) is the center, rc takes the value of 1, 2…Tm, and the point cloud point is determined. rt (x rt ,y rt ),(rt≠rc,and rt∈[1,Tm]) whether the condition is satisfied:

[0097]

[0098] Among them, w l Point cloud point rc The track segment to which it belongs, x′ rt and y′ rt are the distance calculation process quantities along the major axis and minor axis of the ellipse neighborhood, respectively. If f(rc,rt)<1, it is considered that point rt At point rc Otherwise, point rt Not at point rc The statistics are based on point rc (x rc ,y rc ) as the center, the point that meets the above conditions rt The number of point clouds rum(rc), if rum(rc)>minpts(w l ), point rc is valid point cloud data; otherwise, point rc Invalid point cloud data will be removed during fine denoising.

[0099] The point cloud data after fine denoising is recorded as point rf (x rf ,y rf ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction respectively, rf takes the value of 1, 2…Tn, representing the point cloud data sequence, and the total number of point clouds after fine denoising is Tn, (x rf ,y rf ) is the point cloud effective measurement distance screening result.

[0100] The contents not described in detail in this application specification belong to the common knowledge of those skilled in the art.

[0101] The present application is described in detail above in conjunction with specific implementation methods and exemplary examples, but these descriptions cannot be understood as limiting the present application. Those skilled in the art understand that, without departing from the spirit and scope of the present application, a variety of equivalent replacements, modifications or improvements can be made to the technical solution of the present application and its implementation methods, all of which fall within the scope of the present application. The scope of protection of the present application shall be subject to the attached claims.

Claims

1. A method for correcting sea surface ranging errors of a spaceborne single-photon laser radar, characterized in that: include: S1: The single photon flight time is recorded as t. According to the internal clock error Δt1, the transmission pulse delay Δt2 and the pulse detection error Δt3 of the single photon lidar, the single photon flight time after the transceiver link error correction is obtained. S2: The average value t of the detection results obtained by accumulating multiple consecutive single photons emitted from the same target mean With the true value t tar The difference is the system drift error Δt4; the photon flight time of a single photon after instrument correction is S3: Interpolate the latitude and longitude grid sampling data to obtain the temperature T and relative humidity R of the standard pressure layer h ; S4: According to the temperature T and relative humidity R of the standard pressure layer h , calculate the water vapor partial pressure P w ; S5: Write the atmospheric pressure hydrostatic equation in the form of potential height, and construct the relationship between atmospheric pressure P and potential height H gt Relational equation; Solve the above equation to obtain the surface atmospheric pressure P at the laser foot point surf ; S6: According to the water vapor partial pressure P w , the surface atmospheric pressure P at the laser foot point surf , combined with the mapping function model for calculating the atmospheric delay value, determine the dry term mapping function m of the atmospheric delay d (ε) and wet term mapping function m w (ε); S7: Decomposition of atmospheric delay into dry delay ΔL d and wet term delay ΔL w , according to the stem mapping function m d (ε), wet term mapping function m w (ε), dry term delay ΔL d and wet term delay ΔL w , calculate the atmospheric delay error ΔL; According to the photon flight time corrected by the instrument Obtain the measured distance L of the laser after it is refracted by the atmosphere; c is the speed of light; The actual laser straight line measurement distance is Measure distance based on actual laser straight line Obtain photon flight time after atmospheric delay error correction S8: Obtain photon point cloud data of the sea surface based on single-photon laser radar, and calculate each data in the photon point cloud data according to the photon flight time after the single photon is corrected by the atmospheric delay error. And the actual laser straight line measurement distance Perform photon flight time and laser measurement distance correction to obtain the processed point cloud data point r (x r ,y r ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction, respectively, and r takes the value of 1, 2…T r , represents the point cloud data sequence, the total number of point clouds is T r ; S9: Processed point cloud data point r (x r ,y r ) is used for denoising to obtain effective photon point cloud data.

2. The method for correcting sea surface ranging errors of a satellite-borne single-photon laser radar according to claim 1, characterized in that: The internal clock error Δt1 = n·Δt f , where, according to the average time delay Δt measured by each unit of the logic gate circuit delay carry chain in the photon counting electronic acquisition card f and the number of units NE in the delay chain, n < NE and n is an integer.

3. The method for correcting sea surface ranging errors of a space-borne single-photon laser radar according to claim 1, characterized in that: The transmission pulse delay Δt2=Δt o +Δt c , where Δt o is the time delay from the internally generated transmit signal to the laser receiving the transmit signal, Δt c It is the time delay from when the laser receives the transmission signal to when the pulse is emitted.

4. The method for correcting sea surface ranging errors of a spaceborne single-photon laser radar according to claim 1, characterized in that: The pulse detection error Δt3 is calculated as follows: Set the high and low transmission pulse detection thresholds, and record the time when the rising and falling edges of the transmission pulse reach the two detection thresholds respectively. Calculate the difference between the peak values ​​of the transmitted pulses obtained at four time points using the Gaussian fitting method and a single detection threshold: Δt3=t pos1 -t pos2 =ma internal Among them, t pos1 ,t pos2 They are the moments corresponding to the peak values ​​of the transmitted pulse obtained by Gaussian fitting method and single threshold, respectively.

5. The method for correcting sea surface ranging errors of a spaceborne single-photon laser radar according to claim 1, characterized in that: In S6, according to the water vapor partial pressure P w , the surface atmospheric pressure P at the laser foot point surf , combined with the mapping function model for calculating the atmospheric delay value, determine the dry term mapping function m of the atmospheric delay d (ε) and wet term mapping function m w (ε), include: The mapping function model for calculating the atmospheric delay value is the CFA2.2 model; the subscript j = d, w is used to distinguish between dry and wet terms; ε is the laser beam altitude angle, when the laser points to the nadir angle θ = 90-ε < 10°, a j 、b j 、c j The value of is much smaller than 1, and the above formula can be simplified to As the laser pointing nadir angle θ increases, a j =1.185×10 ；3 [1+6.071×10 ；5 (P surf -1000)-1.471×10 ；4 P w +3.072×10 ；3 (T c -20)+1.956×10 ；2 (ξ+6.5)-5.645×10 ；3 (A j -11.231)],b j =1.144×10 ；3 [1+1.164×10 ；5 (P surf -1000)-2.795×10 ；4 P w + 3.109×10 ；3 (T c -20)+3.038×10 ；2 (ξ+6.5)-1.217×10 ；2 (A j -11.231)],c j =-0.009; where ξ is the temperature attenuation rate with altitude, ξ = -6.5℃ / km, A j (j=d,w) is the height parameter of dry and wet items, A d =40.136+0.14872·T c , A w =11, T c is the temperature in °C, with T c =T-273.16;P w is the water vapor partial pressure, P surf is the surface atmospheric pressure.

6. The method for correcting sea surface ranging errors of a spaceborne single-photon laser radar according to claim 1, characterized in that: In S8, the photon point cloud data of the sea surface is obtained according to the single-photon laser radar, and each data in the photon point cloud data is calculated according to the photon flight time after the single photon is corrected by the atmospheric delay error. And the actual laser straight line measurement distance Perform laser measurement distance correction to obtain the processed point cloud data point r (x r ,y r ),include: point or (x or ,y or ) represents the information of the original single photon event record, and defines the earliest received single photon in the extracted point cloud data as point1, point or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or and Tt1 are point or and the photon return time of point 1; vt is the flying speed of the satellite relative to the ground, y or By point or Satellite attitude and orbit at time, photon flight time t or and DEM are obtained by deriving and calculating the rigorous geometric model; point or (x or ,y or ) Photon return time after instrument calibration Photon flight time after instrument correction Single photon point or The flight time is t or ; After atmospheric delay error correction, single photon point or The actual laser straight line measurement distance is ΔL or For single photon point or Atmospheric delay error; Photon flight time after atmospheric delay error correction Based on satellite attitude and orbit, photon return time Photon flight time The DEM is updated by rigorous geometric model calculation or ,y or , get the corrected point cloud data and record point r (x r ,y r ).

7. The method for correcting sea surface ranging errors of a spaceborne single-photon laser radar according to claim 1, characterized in that: In S7, the atmospheric delay error ΔL is: Among them, ρ wa =1.0×10 3 kg / m 3 , is the density of water; P w is the water vapor partial pressure, P surf is the surface atmospheric pressure; m d (ε) is the stem mapping function, m w (v) is the wet term mapping function; M d 、M w are the molecular weight of dry air and water vapor respectively, R=8314.51J / (kmol*K) is the gas constant, g0=9.80665m / s 2 is the average gravitational acceleration in the atmosphere; k1 and k2 are functions of the laser center wavelength λ and are calculated according to the empirical formula given by Owens:

8. The method for correcting sea surface ranging errors of a space-borne single-photon laser radar according to claim 1, characterized in that: In S9, the processed point cloud data point r (x r ,y r ) to denoise and obtain valid photon point cloud data, including: S91: According to the processed point cloud data point r (x r ,y r ), press point r middle y r Take the value and mark the vertical layer as z p , where z is the vertical distance between the centers of each layer, p is the subscript of each layer, p∈[1,M], there are M height layers in total, and the number of photons in each height layer is recorded as num p ; S92: Set the coarse denoising threshold Among them, e∈(1,M / 3) and e is a positive integer, i is the subscript of each altitude layer, traverse each altitude layer, if the corresponding altitude layer num p >th is established, corresponding to the altitude layer z p The data in the middle is the photon point cloud data of the sea surface after rough denoising, denoted as point rc (x rc ,y rc ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction respectively, rc takes the value of 1, 2…Tm, representing the point cloud data sequence, and the total number of point clouds after coarse denoising is Tm; S93: Based on the photon point cloud data point after rough denoising rc Medium rc The value is taken and the segments along the track are marked as w l , where w is the distance along the track from the center of each segment, l is the subscript of each segment, l∈[1,V], there are V segments along the track, and the number of photons extracted from each segment is recorded as lum l , the median extraction result of the photon height of each segment is recorded as zh l ; If the segment w l There is lum l = 0, then the median photon height of the segment is Where Vn≤V, Vn is the number of segments with non-zero photon counts in the along-track segment, li is the segment subscript, and the average number of photons in the segments with non-zero photon counts S94: Each segment zh l By subtracting adjacent values, we can get the median change function of photon height along the track. l∈[1,V-1], for the difference diff(w l ) is smoothed to obtain the sea surface fluctuation function slope(w l ); S95: For the continuous point cloud data of the sea surface, take the neighborhood, set the neighborhood radius eps, and set the neighborhood to an ellipse. The major and minor axes of the neighborhood are e respectively. a =5·eps,e b =0.4·eps, the direction of the neighborhood long axis is based on the sea surface wave function slope(w l ) is determined; according to the number of photons in each segment lum l , and the long axis e of the neighborhood a and short axis e b Determine the core point threshold minpts; S96: According to the core point threshold minpts and the sea surface fluctuation function slope (w l ), and the DBSCAN algorithm with dynamic parameter optimization is used to denoise the photon point cloud data point rc (x rc ,y rc ) for fine denoising, and the point cloud data after fine denoising is recorded as point rf (x rf ,y rf ), where x and y represent the distance information of the point cloud point along the track and in the vertical direction respectively, rf takes the value of 1, 2…Tn, representing the point cloud data sequence, and the total number of point clouds after fine denoising is Tn, (x rf ,y rf ) is the point cloud effective measurement distance screening result.

9. The method for correcting sea surface ranging errors of a space-borne single-photon laser radar according to claim 8, characterized in that: The calculation method of the core point threshold minpts is as follows: Among them, the segment step length Δw = w l:1 -w l , when the track is evenly segmented, Δw does not change with the value of l; minpts(w l ) is for different track segments w l Set the core point threshold, is the average number of photons in the segments where the number of photons is non-zero; e a is the long axis of the neighborhood, e b is the minor axis of the neighborhood.

10. The method for correcting sea surface ranging errors of a space-borne single-photon laser radar according to claim 9, characterized in that: The DBSCAN algorithm for dynamic parameter optimization is as follows: Point cloud after rough denoising rc (x rc ,y rc ) is the center, rc takes the value of 1, 2…Tm, and the point cloud point is determined. rt (x rt ,y rt ),(rt≠rc,and rt∈[1,Tm]) whether the condition is satisfied: Among them, w l Point cloud point rc The track segment to which it belongs, x′ rt and y′ rt are the distance calculation process quantities along the major axis and minor axis of the ellipse neighborhood, respectively. If f(rc,rt)<1, it is considered that point rt At point rc Otherwise, point rt Not at point rc In the neighborhood of rc (x rc ,y rc ) as the center, the point that meets the above conditions rt The number of point clouds rum(rc), if rum(rc)>minpts(w l ), point rc is valid point cloud data; otherwise, point rc Invalid point cloud data will be removed during fine denoising.

Images