A method for correcting sea surface ranging errors of a spaceborne single-photon lidar

Abstract

A method for correcting sea surface ranging errors of a spaceborne single-photon lidar includes: (1) correcting the internal clock error, transmission pulse delay, and pulse detection error of the single-photon lidar; (2) extracting the time centroid compensation system drift by accumulating multiple photon events from a static target to generate a photon histogram; (3) converting the laser elevation data to the geopotential height, interpolating with NCEP related data, calculating the atmospheric pressure at the laser foot point, and completing the NCEP pressure data correction; (4) using the CFA2.2 model to construct an atmospheric delay mapping function adapted to different laser pointing angles, calculating the dry and wet atmospheric delay terms, and correcting the laser measurement distance; (5) completing coarse denoising of the sea surface point cloud data and extracting the sea surface undulation function based on the density distribution characteristics of the sea surface point cloud; (6) constructing a DBSCAN denoising algorithm model with dynamic parameter optimization, completing fine denoising of the point cloud data, and realizing the extraction of effective ranging values ​​of the sea surface point cloud.

Description

Technical Field

[0001] This invention relates to a method for correcting sea surface ranging errors of a spaceborne single-photon lidar, belonging to the field of laser optics. Background Technology

[0002] Photon counting detection systems possess unique advantages in sea surface detection due to their micro-pulse, multi-beam, high sensitivity, and high repetition rate characteristics. The sea surface ranging errors of spaceborne single-photon lidar include transmit / receive link errors, photon detection probability errors, and atmospheric refraction errors. Currently, the atmospheric refraction error mapping function of spaceborne laser altimetry systems does not consider the influence of atmospheric refraction. As the laser pointing angle towards the nadir increases, the delay error of this method increases exponentially, and the influence of atmospheric refraction cannot be ignored. The photon detection probability error is usually calculated using static targets in a laboratory setting. In actual on-orbit dynamic ranging, the introduction of various noises affects the calculation of laser ranging values, requiring error correction based on sea surface point cloud characteristics.

[0003] Currently, the main methods for improving the ranging accuracy of spaceborne single-photon lidar focus on the refined processing of point cloud data. This includes denoising single photons and clustering signal and noisy photons. Commonly used unsupervised classification methods include DBSCAN, LOF, OPTICS, and PSO algorithms, while supervised classification methods include random forests, support vector machines, and conditional random fields. Unsupervised classification methods cluster based on point cloud density, and the key is to determine the clustering parameters related to the target characteristics. Affected by factors such as sea breezes and atmosphere, sea surface single-photon data contains a large number of noisy photons, resulting in significant fluctuations in the sea surface point cloud signal. Furthermore, the point cloud density of signal photons varies along the orbital direction. Extracting and correcting the effective point cloud information from the sea surface is fundamental to maximizing the application effectiveness of lidar products. Summary of the Invention

[0004] The technical problem solved by this application is to overcome the shortcomings of the prior art and provide a method for correcting the ranging error of a spaceborne single-photon lidar on the sea surface, thereby correcting the ranging error of the original point cloud data of the sea surface. This overcomes the problem in the prior art that "when performing dynamic ranging in orbit, the sea surface is affected by factors such as sea wind and atmosphere, resulting in large fluctuations in the point cloud signal and a certain change in the point cloud density along the orbital direction, and the introduction of various noises causes ranging errors".

[0005] This patent proposes a method for correcting sea surface ranging errors using a spaceborne single-photon lidar, comprehensively considering transceiver link errors, photon detection probability errors, atmospheric refraction errors, and sea surface point cloud ranging characteristics. An atmospheric refraction mapping function model is constructed to reduce delay calculation errors when the laser's pointing angle to the sky is large. A dynamically optimized DBSCAN denoising algorithm is designed to adapt to sea surface ranging value calculation.

[0006] The technical solution provided in this application is as follows:

[0007] A method for correcting sea surface ranging errors using a spaceborne single-photon lidar includes the following steps:

[0008] Step 1: Let the uncorrected single-photon flight time be t. Based on the internal clock error Δt1, transmission pulse delay Δt2, and pulse detection error Δt3 of the single-photon lidar, the corrected single-photon flight time after transmit / receive link error correction is calculated. It can be represented as

[0009] Step 2: Use the average value t of 200 consecutive single-photon detection results from the same target. mean Compared with the true value t tar Calculate the system drift error Δt4=t mean -t tar Based on the calculation results in step one, the photon flight time after instrument correction is:

[0010] Step 3: Interpolate the 1°×1° latitude and longitude grid sampling data provided by NCEP. The conversion relationship between geopotential height H and laser footprint elevation h is as follows: Where φ is the geographical latitude, and R e =6371009m is the Earth's radius, g0 = 9.80665m / s 2 Let g be the average gravitational acceleration in the atmosphere, and g be a constant. eq = 9.7803267715 m / s 2 k e =0.001931851353, e 2 =0.00669438002290. In the atmospheric NCEP, the upper and lower layers of the two standard atmospheric pressure layers adjacent to the geopotential height H are denoted as T0, T1, and T2, respectively, based on the temperature, relative humidity, and geopotential height of the upper layer. H0, the temperature, relative humidity, and potential height of the lower layer are denoted as T1, T2, and T3, respectively. H1, temperature and relative humidity vary linearly with geopotential height across different standard atmospheric pressure layers, i.e. Among them, T and R h These are the temperature and relative humidity of the standard pressure level, respectively.

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

[0012]

[0013] In the formula, the constant P b =1000hPa, E s (δ) is a Chebyshev polynomial with 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 interpolated standard pressure layer, and the temperature constant T0 is the temperature constant. max =648K,T min =273K.

[0014] Step 5: Rewrite the atmospheric pressure hydrostatic equation in the form of potential height, and construct the relationship between atmospheric pressure P and potential height H. gt The relational equation is as follows:

[0015]

[0016] in, Let R represent the reciprocal of the compressibility of dry air and the reciprocal of the compressibility of water vapor, respectively. R = 8314.51 J / (kmol*K) is the gas constant. d M w These are the molecular weights of dry air and water vapor, respectively, M. d =28.9644 g / mol, M w =18.0152 g / mol, T is the temperature of the standard barysphere after interpolation, P w The partial pressure of water vapor is given. The fourth-order Runge-Kutta formula is used to solve the above equation, yielding the atmospheric pressure P at the laser trail point. surf .

[0017] Step Six: Calculate the atmospheric delay mapping function. Use the CFA2.2 model (a 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] Where the subscripts j = d and w are used to distinguish between dry and wet terms, ε is the laser beam elevation angle, when the laser beam pointing to the nadir angle θ = 90 - ε < 10°, a j b j cj The value is much less than 1, so the above formula can be simplified to: As the laser pointing angle θ increases, 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 rate of temperature decay with elevation, here taken as ξ = -6.5℃ / km, A j (j=d,w) are the dry and wet height parameters, A d =40.136 + 0.14872·T c A w =11,T c Temperature, in °C, is represented by T. c =T-273.16, P w For water vapor partial pressure, P surf This refers to the atmospheric pressure at the Earth's surface.

[0020] The single-photon lidar of this patent emits a multi-beam laser 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 improvement is more significant as the nadir angle increases.

[0021] Step 7: Decompose the atmospheric delay into the dry term delay ΔL d and wet term delay ΔL w Atmospheric delay error is Where, ρ wa =1.0×103 kg / m 3 P is the density of water. w P surf m d (ε), m w (ε) is obtained from the above calculation, M d M w R and g0 are the known constants mentioned above, and k1 and k2 are functions of the laser center wavelength λ. Calculations are performed using 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] Based on the instrument-corrected photon flight time The measured distance L after the laser beam is refracted by the atmosphere is obtained. c is the speed of light;

[0025] The actual laser linear measurement distance is Based on actual laser linear measurement distance Obtain the photon flight time after atmospheric delay error correction.

[0026] Step 8: Calculate the input sea surface laser point cloud data according to the photon flight time after atmospheric delay error correction for each single photon. And the actual laser linear measurement distance Photon time-of-flight and laser measurement distance corrections are performed, and the processed point cloud data is denoted as point. r (x r ,y r ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, and r takes values ​​of 1, 2...T. r , represents a point cloud data sequence, with a total number of points T. r Press point r middle y r The value is determined by labeling the vertical layers 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 a total of M height layers, and the number of photons in each height layer is denoted as num. p ;

[0027] Step 9: Set the coarse noise reduction threshold Where e∈(1,M / 3) and e is a positive integer, i is the index of each height layer, and the process iterates through each height layer. If the num of the corresponding height layer is... p >th is true, corresponding to the height layer z p The data in the middle is the coarsely denoised sea surface photon point cloud data, denoted as point. rc (x rc ,y rc ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, rc takes values ​​of 1, 2...Tm, representing the point cloud data sequence, and the total number of points in the coarsely denoised cloud is Tm;

[0028] Step 10: Based on the photon point cloud data after coarse denoising in Step 9... rc Chinese x rc The value is determined by dividing the track into segments and marking them as w. l Where w is the distance along the track from the center of each segment, l is the index of each segment, l∈[1,V], there are a total of V segments along the track, and the number of photons extracted from each segment is denoted as lum. l The median photon height extracted from each segment is denoted as zh. l If segmented w l There is lum l If = 0 is true, then the median photon height of that segment is taken as . Where Vn≤V, Vn is the number of segments with non-zero photon counts along the orbital direction, and li is the segment index, representing the average photon count of the non-zero segments.

[0029]

[0030] Step 11: Divide each segment into zh l By subtracting adjacent values, we obtain the function of variation of the median photon height along the orbital direction. l∈[1,V-1], for the difference diff(w) l The sea surface wave function slope(w) is obtained by smoothing the surface. l ).

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

[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 ) after rough denoising. The point cloud data after fine denoising 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 value of rf is 1, 2... Tn, representing the point cloud data sequence. The total number of points in the point cloud after fine denoising is Tn. (x rf ,y rf ) is the screening result of the effective measurement distance of the point cloud.

[0033] In 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] In 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] In 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 high and low emission pulse detection thresholds, 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 the Gaussian fitting method and a single detection threshold for the four moment points:

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

[0041] Among them, t pos1 t pos2 These represent the times corresponding to the peak values ​​of the emitted pulses obtained using Gaussian fitting and single threshold methods, respectively.

[0042] In the above-mentioned method for correcting sea surface ranging errors of a spaceborne single-photon lidar, step five involves using the fourth-order Runge-Kutta formula to calculate the atmospheric pressure P at the laser endpoint. surf The method is as follows:

[0043]

[0044] Set the initial value P0 = P w The number of iterations is N kr Ls is the corresponding integration step size, kr1, kr2, kr3, kr4 are the Runge-Kutta formula parameters, and H... gti P is the potential height at the i-th iteration. 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 using a spaceborne single-photon lidar, step eight involves the following method for correcting sea surface laser point cloud data:

[0046] point or (x or ,y or () represents the information recorded in the original single-photon event. The earliest received single photon in the extracted point cloud data is defined as point1, point... or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or Tt1 and Tt1 are points or The photon return time at point 1; vt is the satellite's velocity relative to the ground, y or By point or Satellite attitude and orbit, photon flight time t or The DEM was obtained through rigorous geometric model derivation and calculation;

[0047] point or (x or ,y or Photon return time after instrument calibration Time of flight of photons after instrument correction single photon point orFlight time is denoted as t or After atmospheric delay error correction, single photon point or The actual laser linear measurement distance is ΔL or For single photon point or Atmospheric delay error; photon flight time corrected for atmospheric delay error Based on satellite attitude and orbit, photon return time Photon flight time And DEM updated by rigorous geometric model calculation x or ,y or Obtain the corrected point cloud data and record it as a point. r (x r ,y r ).

[0048] In the above-mentioned method for correcting sea surface ranging errors using a spaceborne single-photon lidar, the calculation method for the core point threshold minpts in step twelve is as follows:

[0049]

[0050] Where, Δw=w l:1 -w l When the track is divided into uniform segments, Δw does not change with the value of l, minpts(w l For different track sections w l Set the core point threshold, The average number of photons in a non-zero segment.

[0051] In the above-mentioned method for correcting sea surface ranging errors using a spaceborne single-photon lidar, the DBSCAN algorithm for dynamic parameter optimization in step thirteen is as follows:

[0052] Point cloud after coarse denoising rc (x rc ,y rc Centered on ), rc takes values ​​of 1, 2...Tm, to determine the point cloud. rt (x rt ,y rt Does ), (rt≠rc, and rt∈[1,Tm]) satisfy the condition:

[0053]

[0054] Among them, w l For point cloud rc The corresponding segment along the track, x′ rt and y′ rtThese are the distances calculated along the major and minor axes of the ellipse's neighborhood, respectively. If f(rc,rt) < 1, then point is considered... rt At point rc Within the neighborhood of; otherwise, the point is considered to be outside the neighborhood of. rt Not at point rc Within the neighborhood of [point]. Statistics are based on [point]. rc (x rc ,y rc Centered on a point that satisfies the above conditions rt The number of point clouds is rum(rc). If rum(rc) > minpts(w l ), point rc This is valid point cloud data; otherwise, point... rc Invalid point cloud data is removed during fine denoising.

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

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

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

[0058] (3) This invention extracts density parameters along the track direction and combines them with the sea surface wave function to design a DBSCAN denoising algorithm with dynamic parameter optimization. This enables adaptive adjustment of neighborhood direction and discrimination threshold, effectively avoiding the problem of partial signal point cloud loss in the point cloud density denoising algorithm extraction results caused by uneven distribution of point cloud density along the track due to sea wind and fog and sea surface fluctuations. It is more suitable for sea surface ranging value calculation. Attached Figure Description

[0059] Figure 1 A schematic diagram illustrating the ranging error analysis of a spaceborne single-photon lidar facing the sea surface, provided as an embodiment of the present invention;

[0060] Figure 2 This is a flowchart illustrating a method for correcting sea surface ranging errors using a spaceborne single-photon lidar, as provided in an embodiment of the present invention. Detailed Implementation

[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 in conjunction with 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 an 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 time delay from the internal clock to the emission pulse, and the emission pulse time detection error. The photon event detection probability model analyzes the static target system drift error, and 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 the rough denoising of the dynamic target point cloud and the 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 an 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 time 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 time 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 n is an integer;

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

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

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

[0070] Among them, t pos1 t pos2 These represent the times corresponding to the peak values ​​of the emitted pulses obtained using Gaussian fitting and single threshold methods, respectively.

[0071] Step 2: Use the average value t of 200 consecutive single-photon detection results from the same target. mean Compared with the true value t tar Calculate the system drift error Δt4=t mean -t tar Based on the calculation results in step one, the photon flight time after instrument correction is:

[0072] Step 3: Interpolate the 1°×1° latitude and longitude grid sampling data provided by NCEP. The conversion relationship between geopotential height H and laser footprint elevation h is as follows: Where φ is the geographical latitude, and R e =6371009m is the Earth's radius, g0 = 9.80665m / s 2 Let g be the average gravitational acceleration in the atmosphere, and g be a constant. eq = 9.7803267715 m / s 2 k e =0.001931851353, e 2 =0.00669438002290. The upper and lower layers of the two standard atmospheric pressure layers adjacent to H in the atmospheric NCEP are denoted as T0, T1, and T2 respectively, based on the temperature, relative humidity, and geopotential height of the upper layer. H0, the temperature, relative humidity, and potential height of the lower layer are denoted as T1, T2, and T3, respectively. H1, temperature and relative humidity vary linearly with geopotential height across different standard atmospheric pressure layers, i.e.

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

[0074]

[0075] In the formula, the constant P b =1000hPa, E s (δ) is a Chebyshev polynomial with 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 interpolated standard pressure layer, and the temperature constant T0 is the temperature constant. max =648K,T min =273K.

[0076] Step 5: Rewrite the atmospheric pressure hydrostatic equation in the form of potential height, and construct the relationship between atmospheric pressure P and potential height H. gt The relational equation is as follows:

[0077]

[0078] in, Let R represent the reciprocal of the compressibility of dry air and the reciprocal of the compressibility of water vapor, respectively. R = 8314.51 J / (kmol*K) is the gas constant. d M w These are the molecular weights of dry air and water vapor, respectively, M. d =28.9644 g / mol, M w =18.0152 g / mol, T is the temperature of the standard barysphere after interpolation, P w The partial pressure of water vapor is given. The atmospheric pressure P at the laser trail point is obtained using the fourth-order Runge-Kutta formula. surf The method is as follows:

[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 size, kr1, kr2, kr3, kr4 are the Runge-Kutta formula parameters, H gti P is the potential height at the i-th iteration. i The potential height is H gti Atmospheric pressure at that time, take

[0081] Step Six: 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] Where the subscripts j = d and w are used to distinguish between dry and wet terms, ε is the laser beam elevation angle, when the laser beam pointing to the nadir angle θ = 90 - ε < 10°, a j b j c j The value is much less than 1, so the above formula can be simplified to: As the laser pointing angle θ increases, 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 rate of temperature decay with elevation, here taken as ξ = -6.5℃ / km, A j (j=d,w) are the dry and wet height parameters, A d =40.136 + 0.14872·T c A w =11,T c Temperature, in °C, is represented by T. c =T-273.16, P w For water vapor partial pressure, P surf This refers to the atmospheric pressure at the Earth's surface.

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

[0085]

[0086] When λ = 1064 nm, k1 = 0.80277 K / Pa, k2 = 0.66388 K / Pa; when λ = 532 nm, k1 = 0.77493 K / Pa, k2 = 64873 K / Pa, based on the photon flight time after instrument calibration. Obtain the measurement distance after laser refraction through the atmosphere c is the speed of light, and the actual linear distance measured by laser is... Based on actual laser linear measurement distance Obtain the photon flight time after atmospheric delay error correction.

[0087] Step 8: Calculate the input sea surface laser point cloud data according to the photon flight time after atmospheric delay error correction for each single photon. And the actual laser linear measurement distance Perform photon flight time and laser measurement distance correction. or (x or ,y or () represents the information recorded in the original single-photon event. The earliest received single photon in the extracted point cloud data is defined as point1, point... or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or Tt1 and Tt1 are points or The photon return time at point 1; vt is the satellite's velocity relative to the ground, y or By point or Satellite attitude and orbit, photon flight time tor The DEM was obtained through rigorous geometric model derivation and calculation;

[0088] point or (x or ,y or Photon return time after instrument calibration Time of flight of photons after instrument correction single photon point or Flight time is denoted as t or After atmospheric delay error correction, single photon point or The actual laser linear measurement distance is ΔL or For single photon point or Atmospheric delay error; photon flight time corrected for atmospheric delay error Based on satellite attitude and orbit, photon return time Photon flight time And DEM updated by rigorous geometric model calculation x or ,y or Obtain the corrected point cloud data and record it as a point. r (x r ,y r ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, and r takes values ​​of 1, 2...T. r , represents a point cloud data sequence, with a total number of points T. r Press point r middle y r The value is determined by labeling the vertical layers 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 a total of M height layers, and the number of photons in each height layer is denoted as num. p ;

[0089] Step 9: Set the coarse noise reduction threshold Where e∈(1,M / 3) and e is a positive integer, i is the index of each height layer, and the process iterates through each height layer. If the num of the corresponding height layer is... p >th is true, corresponding to the height layer z p The data in the middle is the coarsely denoised sea surface photon point cloud data, denoted as point. rc (x rc ,y rc ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, rc takes values ​​of 1, 2...Tm, representing the point cloud data sequence, and the total number of points in the coarsely denoised cloud is Tm;

[0090] Step 10: Based on the photon point cloud data after coarse denoising in Step 9... rc Chinese x rc The value is determined by dividing the track into segments and marking them as w. l Where w is the distance along the track from the center of each segment, l is the index of each segment, l∈[1,V], there are a total of V segments along the track, and the number of photons extracted from each segment is denoted as lum. l The median photon height extracted from each segment is denoted as zh. l If segmented w l There is lum l If = 0 is true, then the median photon height of that segment is taken as . Where Vn≤V, Vn is the number of segments with non-zero photon counts along the orbital direction, and li is the segment index, representing the average photon count of the non-zero segments.

[0091]

[0092] Step 11: Divide each segment into zh l By subtracting adjacent values, we obtain the function of variation of the median photon height along the orbital direction. l∈[1,V-1], for the difference diff(w) l The sea surface wave function slope(w) is obtained by smoothing the surface. l ).

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

[0094]

[0095] Where, Δw=w l:1 -w l When the track is divided into uniform segments, Δw does not change with the value of l, minpts(w l For different track sections w l Set the core point threshold, The average number of photons in a non-zero segment.

[0096] Step 13: Substitute the above parameter values ​​and use the dynamically optimized DBSCAN algorithm to process the coarsely denoised photon point cloud data. rc (xrc ,y rc Fine denoising is performed, and the point cloud after coarse denoising is used. rc (x rc ,y rc Centered on ), rc takes values ​​of 1, 2...Tm, to determine the point cloud. rt (x rt ,y rt Does ), (rt≠rc, and rt∈[1,Tm]) satisfy the condition:

[0097]

[0098] Among them, w l For point cloud rc The corresponding segment along the track, x′ rt and y′ rt These are the distances calculated along the major and minor axes of the ellipse's neighborhood, respectively. If f(rc,rt) < 1, then point is considered... rt At point rc Within the neighborhood of; otherwise, the point is considered to be outside the neighborhood of. rt Not at point rc Within the neighborhood of [point]. Statistics are based on [point]. rc (x rc ,y rc Centered on a point that satisfies the above conditions rt The number of point clouds is rum(rc). If rum(rc) > minpts(w l ), point rc This is valid point cloud data; otherwise, point... rc Invalid point cloud data is removed during fine denoising.

[0099] The point cloud data after fine denoising is denoted as point rf (x rf ,y rf ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, rf takes values ​​of 1, 2...Tn, representing the point cloud data sequence, and the total number of points in the finely denoised cloud is Tn, (x rf ,y rf This is the result of filtering the effective measurement distance of the point cloud.

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

[0101] The present application has been described in detail above with reference to specific embodiments and exemplary examples; however, these descriptions should not be construed as limiting the present application. Those skilled in the art will understand that various equivalent substitutions, modifications, or improvements can be made to the technical solutions and implementation methods of the present application without departing from the spirit and scope of the present application, and all such modifications and improvements fall within the scope of the present application. The scope of protection of the present application is determined by the appended claims.

Claims

1. A method for correcting sea surface ranging errors using a spaceborne single-photon lidar, characterized in that, include: S1: The single-photon flight time is denoted as t. Based on 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 correction for transmit / receive link errors is obtained. S2: The average value t of the detection results obtained by accumulating multiple consecutive single photons emitted from the same target. mean Compared with the true value t tar The difference is the system drift error Δt4; the photon flight time 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: Based on the temperature T and relative humidity R of the standard pressure level h The partial pressure of water vapor P was calculated. w ; S5: Rewrite the atmospheric pressure hydrostatic equation in the form of potential height, and construct the relationship between atmospheric pressure P and potential height H. gt Relationship equation; Solve the above equation to obtain the atmospheric pressure P at the laser footpoint. surf ; S6: Based on the partial pressure of water vapor P w Atmospheric pressure P at the laser footprint surf The atmospheric delay interference mapping function m is determined by combining the mapping function model for calculating atmospheric delay values. d (ε) and wetted term mapping function m w (ε); S7: Decompose atmospheric delay into dry term delay ΔL d and wet term delay ΔL w According to the term mapping function m d (ε), wetted term mapping function m w (ε), Delayed term ΔL d and wet term delay ΔL w Calculate the atmospheric delay error ΔL; Based on the instrument-corrected photon flight time Obtain the measured distance L of the laser beam after atmospheric refraction; c is the speed of light; The actual laser linear measurement distance is Based on actual laser linear measurement distance Obtain the photon flight time after atmospheric delay error correction. S8: Based on the photon point cloud data of the sea surface obtained by single-photon lidar, the photon flight time of each data point in the photon point cloud data is calculated according to the photon flight time after atmospheric delay error correction. And the actual laser linear measurement distance Photon time-of-flight and laser measurement distance corrections are performed to obtain processed point cloud data. r (x r ,y r ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, and r takes values ​​of 1, 2...T. r , represents a point cloud data sequence, with a total number of points T. r ; S9: Process the point cloud data. r (x r ,y r Denoising is performed to obtain effective photon point cloud data.

2. The method for correcting sea surface ranging errors of a spaceborne single-photon lidar 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 spaceborne single-photon lidar according to claim 1, characterized in that: The transmission pulse delay Δt2=Δt o +Δt c , where Δt o Δt is the time delay between the internally generated transmitted signal and the laser's reception of the transmitted signal. c This refers to the time delay between the laser receiving the transmitted signal and the pulse being emitted.

4. The method for correcting sea surface ranging errors of a spaceborne single-photon lidar according to claim 1, characterized in that, The method for calculating the pulse detection error Δt3 is as follows: Set two emission pulse detection thresholds, one for high and one for low pulses, and record the times when the rising and falling edges of the emission pulses reach the two detection thresholds respectively. Calculate the difference in peak emission pulse values ​​at four time points obtained using Gaussian fitting and a single detection threshold: Δt3=t pos1 -t pos2 = inside ma Among them, t pos1 t pos2 These represent the times corresponding to the peak values ​​of the emitted pulses obtained using Gaussian fitting and single threshold methods, respectively.

5. The method for correcting sea surface ranging errors of a spaceborne single-photon lidar according to claim 1, characterized in that, In S6, based on the partial pressure of water vapor P w Atmospheric pressure P at the laser footprint surf The atmospheric delay interference mapping function m is determined by combining the mapping function model for calculating atmospheric delay values. d (ε) and wetted term mapping function m w (ε), include: The mapping function model for calculating the atmospheric delay value is the CFA2.2 model; the subscripts j = d, w are used to distinguish between dry and wet terms; ε is the laser beam elevation angle, when the laser pointing to the nadir angle θ = 90 - ε < 10°, a j b j c j The value is much less than 1, so the above formula can be simplified to: As the laser pointing angle θ towards the nadir 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 rate of temperature decay with elevation, ξ = -6.5℃ / km, A j (j=d,w) are the dry and wet height parameters, A d =40.136 + 0.14872·T c A w =11,T c Temperature, in °C, is represented by T. c =T-273.16; P w For water vapor partial pressure, P surf This refers to the atmospheric pressure at the Earth's surface.

6. The method for correcting sea surface ranging errors of a spaceborne single-photon lidar according to claim 1, characterized in that, In step S8, based on the photon point cloud data of the sea surface obtained by the single-photon lidar, each data point in the photon point cloud data is processed according to the photon flight time after atmospheric delay error correction. And the actual laser linear measurement distance Laser measurement distance correction is performed to obtain processed point cloud data. r (x r ,y r ),include: point or (x or ,y or () represents the information recorded in the original single-photon event. The earliest received single photon in the extracted point cloud data is defined as point1, point... or The distance relative to point1 along the flight direction is x or =(Tt or -Tt1)·vt, where Tt or Tt1 and Tt1 are points or The photon return time at point 1; vt is the satellite's velocity relative to the ground, y or By point or Satellite attitude and orbit, photon flight time t or The DEM was obtained through rigorous geometric model derivation and calculation; point or (x or ,y or Photon return time after instrument calibration Time of flight of photons after instrument correction single photon point or Flight time is denoted as t or After atmospheric delay error correction, single photon point or The actual laser linear measurement distance is ΔL or For single photon point or Atmospheric delay error; photon flight time corrected for atmospheric delay error Based on satellite attitude and orbit, photon return time Photon flight time And DEM updated by rigorous geometric model calculation x or ,y or Obtain the corrected point cloud data and record it as a point. r (x r ,y r ).

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

8. The method for correcting sea surface ranging errors of a spaceborne single-photon lidar according to claim 1, characterized in that, In step S9, the processed point cloud data is... r (x r ,y r Denoising is performed to obtain effective photon point cloud data, including: S91: Based on the processed point cloud data r (x r ,y r ), press point r middle y r The value is determined by labeling the vertical layers 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 a total of M height layers, and the number of photons in each height layer is denoted as num. p ; S92: Set coarse noise reduction threshold Where e∈(1,M / 3) and e is a positive integer, i is the index of each height layer, and the process iterates through each height layer. If the num of the corresponding height layer is... p >th is true, corresponding to the height layer z p The data in the middle is the coarsely denoised sea surface photon point cloud data, denoted as point. rc (x rc ,y rc ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, rc takes values ​​of 1, 2...Tm, representing the point cloud data sequence, and the total number of points in the coarsely denoised cloud is Tm; S93: Based on the coarsely denoised photon point cloud data rc Chinese x rc The value is determined by dividing the track into segments and marking them as w. l Where w is the distance along the track from the center of each segment, l is the index of each segment, l∈[1,V], there are a total of V segments along the track, and the number of photons extracted from each segment is denoted as lum. l The median photon height extracted from each segment is denoted as zh. l If segmented w l There is lum l If = 0 is true, then the median photon height of that segment is taken as . Where Vn≤V, Vn is the number of segments with non-zero photon counts along the orbital direction, and li is the segment index, representing the average photon count of the segments with non-zero photon counts. S94: Translate each segment zh l By subtracting adjacent values, we obtain the function of variation of the median photon height along the orbital direction. l∈[1,V-1], for the difference diff(w) l The sea surface wave function slope(w) is obtained by smoothing the surface. l ); S95: For continuous point cloud data of the sea surface, take the neighborhood, set the neighborhood radius eps, and set the neighborhood as an ellipse with the major axis and minor axis of the neighborhood being e. a =5·eps,e b =0.4·eps, the major axis direction of the neighborhood is based on the sea surface wave function slope(w) l Determine; based on the number of photons in each segment (lum). l and the major axis e of the neighborhood a and short axis e b Determine the core point threshold minpts; S96: Based on the core point threshold minpts and the sea surface undulation function slope(w) l The DBSCAN algorithm with dynamic parameter optimization is used to process the coarsely denoised photon point cloud data. rc (x rc ,y rc Perform fine denoising, and denote the denoised point cloud data as a point. rf (x rf ,y rf ), where x and y represent the distance information of the point cloud along the track and vertical direction, respectively, rf takes values ​​of 1, 2...Tn, representing the point cloud data sequence, and the total number of points in the finely denoised cloud is Tn, (x rf ,y rf This is the result of filtering the effective measurement distance of the point cloud.

9. A method for correcting sea surface ranging errors using a spaceborne single-photon lidar according to claim 8, characterized in that, The core point threshold minpts is calculated as follows: Where, the piecewise step size Δw = w l:1 -w l When the track is divided into uniform segments, Δw does not change with the value of l; minpts(w l For different track sections w l Set the core point threshold, e is the average photon count for segments with non-zero photon counts. a e is the major axis of the neighborhood. b It is the short axis of the neighborhood.

10. A method for correcting sea surface ranging errors using a spaceborne single-photon lidar according to claim 9, characterized in that, The dynamically optimized DBSCAN algorithm is as follows: Point cloud after coarse denoising rc (x rc ,y rc Centered on ), rc takes values ​​of 1, 2...Tm, to determine the point cloud. rt (x rt ,y rt Does ), (rt≠rc, and rt∈[1,Tm]) satisfy the condition: Among them, w l For point cloud rc The corresponding segment along the track, x′ rt and y′ rt These are the distances calculated along the major and minor axes of the ellipse's neighborhood, respectively. If f(rc,rt) < 1, then point is considered... rt At point rc Within the neighborhood of; otherwise, the point is considered to be outside the neighborhood of. rt Not at point rc Within the neighborhood; statistics are based on points rc (x rc ,y rc Centered on a point that satisfies the above conditions rt The number of point clouds is rum(rc). If rum(rc) > minpts(w l ), point rc This is valid point cloud data; otherwise, point cloud data is invalid. rc Invalid point cloud data is removed during fine denoising.

Images