A method for inverting water optical parameters using spaceborne single-photon lidar

Abstract

The present invention relates to a method for inverting water body optical parameters of a satellite-borne single-photon laser radar, and belongs to the field of laser remote sensing technology. First, the effective signal of the satellite-borne single-photon laser radar when flying above the water surface and measuring is extracted, and the afterpulse effect of the single-photon laser radar is eliminated by a deconvolution algorithm. By establishing a theoretical model of water body backscattering, the water body backscattering coefficient and the diffuse attenuation coefficient are calculated by using the least squares fitting method and the water body backscattering laser pulse point cloud obtained after eliminating the afterpulse effect, and the water body backscattering coefficient is calculated according to the water body parameter backscattering model. The present invention can quickly and accurately calculate the main optical parameters of the seawater subsurface layer of the local water body through the point cloud data obtained when the satellite-borne single-photon laser radar flies over the ocean, which plays an important role in effectively monitoring water pollution, water eutrophication, water quality degradation, etc.

Description

Technical Field

[0001] The invention belongs to the technical field of laser remote sensing, and in particular relates to a method for inverting water body optical parameters of a satellite-borne single-photon laser radar. Background Art

[0002] For the study of the ocean subsurface, in-situ measurements by shipborne and buoy methods can provide high-precision results, but they are expensive, have low detection efficiency, and are difficult to implement over a large area; satellite-borne passive optical remote sensing methods are restricted by solar radiation and satellite orbits, making it difficult to work in high-latitude areas and unable to provide effective vertical profile distribution. Active detection lidar remote sensing benefits from the active emission of laser pulses and can work all day long. Satellite-borne lidar has the advantages of wide coverage and has gradually become a reliable tool for research. The extreme detection capability and detection efficiency of traditional linear lidar are limited, while single-photon detectors such as Geiger-mode Avalanche Photodiode (GM-APD) and photomultiplier tube (PMT) can respond to the energy of a single photon, which is about 1000 times higher than the detection sensitivity of traditional linear detectors, greatly improving the performance of lidar.

[0003] In the inversion of optical parameters of the ocean subsurface, there are relatively complete theories, methods and experiments using traditional linear laser radar to obtain optical parameters of the ocean subsurface, that is, using the full waveform data obtained by the traditional linear laser radar, by establishing a physical or empirical model between the waveform parameters and the water body parameters, extracting and analyzing the main parameters of the measured waveform, and inverting the water body attenuation coefficient and scattering coefficient. However, the water body optical parameter inversion method of linear laser radar is not suitable for the single photon system. This is because: the single photon laser radar does not have the waveform recording function, and the result of each detection is a discrete photon event. It can only record the presence or absence of the signal, but cannot identify the size of the signal; and because the single photon detector has a high sensitivity and is greatly affected by the background light, the original data contains a large number of noise photons, while only the underwater backscattered signal photons are needed for the inversion of water body optical parameters; in addition, the single photon detector is affected by the afterpulse effect, which will bring additional noise that the traditional linear detector does not have, and needs to be corrected in practical applications. Summary of the invention

[0004] The purpose of the present invention is to provide a method for inverting the optical parameters of water bodies of a space-borne single-photon laser radar, so as to solve the problem in the prior art that the method for inverting the optical parameters of water bodies of a linear laser radar is not applicable to a single-photon system.

[0005] A method for inverting water body optical parameters of a space-borne single-photon laser radar, comprising:

[0006] S1. Acquire the effective signal of the spaceborne single-photon laser radar when it flies above the water surface and performs measurements;

[0007] S2. Eliminate the after-pulse effect of single-photon lidar using deconvolution algorithm;

[0008] S3. Based on the theoretical model of water backscattering, the water backscattering coefficient and water diffuse attenuation coefficient are obtained by least square fitting;

[0009] S4. Calculate the backscattering coefficient of the water body.

[0010] Preferably, the effective signal is a water body backscattered point cloud, that is, a point cloud formed by a water body backscattered laser pulse. A water surface reflection point cloud with a higher spatial point cloud density is extracted through a spatial point cloud density classification method. The water surface reflection point cloud is a point cloud formed by a water surface reflecting a laser pulse. Point clouds with higher height values ​​than the water surface reflection point cloud correspond to point clouds formed by atmospheric scattered laser pulses, and point clouds with lower height values ​​than the water surface reflection point cloud correspond to water body backscattered point clouds.

[0011] Preferably, step S1 includes the following sub-steps:

[0012] S1.1. Calculate the point cloud density within the unit space;

[0013] S1.2. Calculate the mean value μ and standard deviation σ within the spatial grid;

[0014] S1.3. Calculate the discrimination threshold TH between the water surface reflection point cloud and the atmosphere and water body backscattering point cloud;

[0015] S1.4. Calculate the average backscattered laser pulse point cloud sequence y(z) = [y(z1), y(z2), …, y(z n )], z represents the depth sequence, z1,z2,…,z n Represents different depths, and the statistical interval in the depth direction is Δz=z2-z1=z3-z2=z4-z3.

[0016] Preferably, step S1.1 includes: the length interval of each spatial grid along the satellite flight direction is Δl, and the interval in the height direction is Δh. The total length of the data acquired by the satellite-borne single-photon laser radar in the study area along the satellite flight direction is L, and the total length in the height direction is H. Then the number of grids along the satellite flight direction is n l and the number of grids in the height direction n h They can be expressed as: n l =ceil(L / Δl),nh =ceil(H / Δh), ceil represents the upward rounding function; the data in the study area are divided into spatial grids, and the spatial coordinate range F(i,j) of the i-th grid along the satellite flight direction and the j-th grid along the height direction satisfies:

[0017]

[0018] L start and L end are the starting coordinate and ending coordinate in the study area along the satellite flight direction, H min and H max They are the height direction recording range in the study area, that is, the minimum and maximum values ​​of the distance window of the lidar; the number of photon point clouds within each spatial grid coordinate range F(i,j) is counted, that is, the point cloud density within the unit space, recorded as N(i,j).

[0019] Preferably, in step S1.2,

[0020] Preferably, in step S1.3, TH = μ + e a ×σ,e a is the proportional coefficient. In the study area, all point clouds within the spatial range where the spatial grid density N(i, j) is greater than the discrimination threshold TH are classified as water surface reflection point clouds.

[0021] Preferably, in step S2, y(z) is a real water backscatter signal sequence x(z)=[x(z1),x(z2),…,x(z n )] is the convolution of the impulse response function H(z) of the single-photon lidar system, that is, y(z)=H(z)x(z), where z represents the depth sequence, and y(z)=H(z)x(z) is equivalent to:

[0022]

[0023] z1,z2,…,z n Represents different depths.

[0024] Preferably, the pulse response function H(z) of the single-photon laser radar system can be obtained by calculating the received signal of the laser pulse incident on the planar target, or by calculating the received signal corresponding to a portion of the laser pulse energy intercepted inside the single-photon laser radar system. When the pulse response function H(z) is known, the deconvolution x(z)=H -1 (z)y(z) obtains the real underwater backscatter signal sequence x(z) after eliminating the after-pulse effect, where H -1 (z) represents the inverse matrix of H(z).

[0025] Preferably, the water body backscatter signal sequence x(z) and the water body parameter backscatter theoretical function can be expressed as:

[0026]

[0027] exp is the exponential function, sec is the secant trigonometric function; β(π) represents the volume scattering coefficient when the angle is π (backward), K d represents the water diffuse attenuation coefficient, which is the two parameters to be calculated; E t is the single pulse laser emission energy, η is the comprehensive efficiency of the single photon lidar system, h p is the Planck constant, which is approximately equal to 6.63×10-34 J·s, v is the laser wavelength efficiency, R is the flight altitude of the satellite-borne single-photon laser radar, and A r is the effective area of ​​the receiving aperture of the laser radar telescope. The above parameters are known parameters for a given laser radar system; n w is the refractive index of water, θ is the refraction angle of the emitted laser pulse after entering the water, T a is the attenuation coefficient of the laser pulse when passing through the atmosphere, T w is the attenuation coefficient of the laser pulse when it passes through the water surface and enters the water body. The above parameters are usually constant values ​​under given conditions or can be obtained by calculating the laser radar system parameters;

[0028] In the above formula, there are two parameters to be substituted. The data in the backscatter signal sequence x(z) is usually greater than 10, and there is a data redundancy problem. In order to obtain more stable results and suppress the random errors in the data, the least squares fitting method is used to calculate the backscattering coefficient β(π) and the diffuse attenuation coefficient K d .

[0029] Preferably, the water body backscattering coefficient in step S4 is a more commonly used water body optical parameter, and the water body backscattering coefficient b is calculated using the backscattering coefficient β(π) obtained in step 3. b , assuming that the backscattering is isotropic, the backscattering coefficient b b Defined as the sum of the volume scattering coefficients in the hemisphere space behind the incident light: b b =2πβ(π).

[0030] Compared with the existing technology, the present invention establishes a method for inverting the optical parameters of water bodies through theoretical modeling. It only needs to give the system hardware parameters and measurement data of the single-photon lidar, combined with the signal statistical model and the water parameter backscattering model, to accurately invert the optical parameters of the sub-surface water body of the measured water area, which plays an important role in effectively monitoring water pollution, eutrophication, water quality degradation, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is a technical flow chart of the method of the present invention.

[0032] Figure 2 This is the distribution map of photon point cloud obtained by satellite-borne single-photon lidar in nighttime measurements in open ocean areas. The horizontal axis is the relative cumulative distance along the satellite's flight orbit, and the vertical axis is the altitude direction. The ellipsoid elevation datum is used, and the units are all meters. The grid in the figure is the spatial grid division of the photon point cloud. The spatial grid size is Δl in length interval along the satellite's flight direction and Δh in altitude interval. The right picture is a partial enlargement of the left picture.

[0033] Figure 3 is the water depth value z calculated by using the point cloud data collected by ICESat-2 / ATLAS in the Sahara Desert at night. i The corresponding impulse response function h(z i ) between the curves. DETAILED DESCRIPTION

[0034] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention is described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0035] The present invention uses the system parameters and acquired data of the world's first satellite-borne single-photon laser radar ATLAS (Advance Topographic Laser Altimeter System) carried by the US ICESat-2 (Ice, Cloud, and land Elevation Satellite-2) as an example. ICESat-2 / ATLAS has 6 laser beams, a laser pulse repetition frequency of 10kHz, a laser sampling interval of 0.7 meters along the track, and a spot diameter of 12 meters. The multi-beam sampling mode makes the spatial point cloud more dense, which will provide higher authenticity for rough and uneven terrain. The ATL03 data product provided by ICESat-2 / ATLAS gives the point cloud data acquired along the satellite flight trajectory. The point cloud data contains the signal points and noise points acquired by ICESat-2 / ATLAS. Each point cloud has a unique three-dimensional coordinate (longitude, latitude and altitude), and the coordinate system is based on the WGS84 ellipsoid. The sample data of the embodiment selects the point cloud data acquired by the ICESat-2 / ATLAS single-photon lidar flying over the Indian Ocean at 20:38:55 (local time) on December 13, 2018.

[0036] A method for inverting water body optical parameters of a space-borne single-photon laser radar, comprising:

[0037] S1. Acquire the effective signal of the spaceborne single-photon laser radar when it flies above the water surface and performs measurements;

[0038] S2. Eliminate the after-pulse effect of single-photon lidar using deconvolution algorithm;

[0039] S3. Based on the theoretical model of water backscattering, the water backscattering coefficient and water diffuse attenuation coefficient are obtained by least square fitting;

[0040] S4. Calculate the backscattering coefficient of the water body.

[0041] Preferably, the effective signal is a water body backscattered point cloud, that is, a point cloud formed by a water body backscattered laser pulse. A water surface reflection point cloud with a higher spatial point cloud density is extracted through a spatial point cloud density classification method. The water surface reflection point cloud is a point cloud formed by a water surface reflecting a laser pulse. Point clouds with higher height values ​​than the water surface reflection point cloud correspond to point clouds formed by atmospheric scattered laser pulses, and point clouds with lower height values ​​than the water surface reflection point cloud correspond to water body backscattered point clouds.

[0042] Preferably, step S1 includes the following sub-steps:

[0043] S1.1. Calculate the point cloud density within the unit space;

[0044] S1.2. Calculate the mean value μ and standard deviation σ within the spatial grid;

[0045] S1.3. Calculate the discrimination threshold TH between the water surface reflection point cloud and the atmosphere and water body backscattering point cloud;

[0046] S1.4. Calculate the average backscattered laser pulse point cloud sequence y(z) = [y(z1), y(z2), …, y(z n )], z represents the depth sequence, z1,z2,…,z n Represents different depths, and the statistical interval in the depth direction is Δz=z2-z1=z3-z2=z4-z3.

[0047] When detecting a target, single-photon laser radar cannot record waveforms like traditional linear detectors, but it can record the three-dimensional coordinates of each photon point cloud. Figure 2The distribution map of photon point clouds obtained by satellite-borne single-photon laser radar in nighttime measurements in open ocean areas is shown; the horizontal axis is the relative cumulative distance along the satellite flight orbit, and the vertical axis is the height direction, using the ellipsoid elevation datum, and the units are all meters. When detecting water targets, the photon signals received by the satellite-borne single-photon laser radar are mainly formed by point clouds formed by atmospheric scattered laser pulses, point clouds formed by water surface reflected laser pulses, point clouds formed by water body backscattered laser pulses, point clouds formed by water bottom reflected laser pulses, point clouds formed by atmospheric scattered solar background radiation, point clouds formed by water surface reflected solar background radiation, point clouds formed by water body backscattered solar background radiation, point clouds formed by water bottom reflected solar background radiation, point clouds formed by detector dark counts, and point clouds formed by detector residual pulses. When obtaining water body optical parameters through single-photon laser radar, only the point cloud formed by water body backscattered laser pulses belongs to the signal, and all point clouds related to solar background radiation belong to noise, and the point clouds formed by atmospheric scattering, water surface reflection and water bottom reflected laser pulses belong to noise, and the point clouds formed by detector dark counts and laser radar system residual pulses also belong to noise. Among them, the influence of detector dark count is very small and can be ignored; the influence of the residual pulse of the lidar system is coupled into the point cloud formed by the backscattered laser pulse of the water body, and it is necessary to extract the point cloud formed by the backscattered laser pulse of the water body first.

[0048] In view of the above problems, in order to better extract the effective water body backscatter signal point cloud, the present invention adopts the following methods to process and screen the data: (1) select data measured at night to exclude all influences related to solar background radiation; (2) select data measured in the open sea away from islands, reefs and coastal zones to exclude the influence of underwater reflected laser pulses; (3) through the spatial point cloud density classification method, extract the point cloud formed by the water surface reflected laser pulse with higher spatial point cloud density, the point cloud with higher height values ​​than these water surface point clouds corresponds to the point cloud formed by the atmospheric scattered laser pulse, and the point cloud with lower height values ​​than these water surface point clouds corresponds to the point cloud formed by the water body backscattered laser pulse.

[0049] Preferably, step S1.1 includes: the length interval of each spatial grid along the satellite flight direction is Δl, and the interval in the height direction is Δh. The total length of the data acquired by the satellite-borne single-photon laser radar in the study area along the satellite flight direction is L, and the total length in the height direction is H. Then the number of grids along the satellite flight direction is n l and the number of grids in the height direction n h They can be expressed as: n l =ceil(L / Δl),n h =ceil(H / Δh), ceil represents the upward rounding function; the data in the study area are divided into spatial grids, and the spatial coordinate range F(i,j) of the i-th grid along the satellite flight direction and the j-th grid along the height direction satisfies:

[0050]

[0051] L start and L end are the starting coordinate and ending coordinate in the study area along the satellite flight direction, H min and H max They are the height direction recording range in the study area, that is, the minimum and maximum values ​​of the distance window of the lidar; the number of photon point clouds within each spatial grid coordinate range F(i,j) is counted, that is, the point cloud density within the unit space, recorded as N(i,j).

[0052] Preferably, in step S1.2,

[0053] Preferably, in step S1.3, since the water surface reflection coefficient is much larger than the backscattering coefficient of the atmosphere and water body, the point cloud in the grid with higher spatial point cloud density corresponds to the water surface reflection point cloud, and based on this, the discrimination threshold TH=μ+e between the water surface point cloud and the backscattering point cloud of the atmosphere and water body is calculated. a ×σ,e a is the proportional coefficient. In the study area, all point clouds within the spatial range where the spatial grid density N(i, j) is greater than the discrimination threshold TH are classified as water surface reflection point clouds.

[0054] Preferably, in step S1.4, since the point cloud signal formed by the backscattered laser pulse of the water body is weak (the number of point clouds is small) and a single measurement has a certain randomness, after accumulating multiple consecutive pulses in the direction of the track (satellite flight), the backscattered laser pulse point cloud sequence y(z)=[y(z1),y(z2),…,y(z n )], z represents the depth sequence, z1,z2,…,z n Represents different depths, and the depth direction statistical interval Δz=z2-z1=z3-z2=z4-z3. Along the satellite flight direction, the depth zero point at different positions is equal to the lower limit value of the spatial grid coordinate classified as the water surface point cloud at the current position.

[0055] Preferably, in step S2, y(z) is a real water backscatter signal sequence x(z)=[x(z1),x(z2),…,x(z n )] is the convolution of the impulse response function H(z) of the single-photon lidar system, that is, y(z)=H(z)x(z), where z represents the depth sequence, and y(z)=H(z)x(z) is equivalent to:

[0056]

[0057] z1,z2,…,z n Represents different depths.

[0058] Preferably, the pulse response function H(z) of the single-photon laser radar system can be obtained by calculating the received signal of its laser pulse incident on a planar target, or by calculating the received signal corresponding to a portion of the laser pulse energy intercepted inside the single-photon laser radar system. In this embodiment, there are two methods to obtain the pulse response function H(z) of ICESat-2 / ATLAS. First, the internal TEP (Transmitter Echo Path) module of the ICESat-2 / ATLAS system can be used to intercept part of the laser pulse energy internally to form a corresponding received signal to obtain H(z). In addition, the night measurement data obtained by ICESat-2 / ATLAS in the Sahara Desert area can be selected because: (1) the data obtained by night measurement excludes all influences related to solar background radiation; (2) the laser pulses with a wavelength of 532nm emitted by ICESat-2 / ATLAS cannot penetrate the desert surface; (3) the surface cover of the Sahara Desert (such as vegetation) is very small, and the surface is relatively flat in terms of the size of the 12m laser spot. The point cloud data obtained by ICESat-2 / ATLAS in the Sahara Desert at night are accumulated along the satellite flight track and then layered along the depth direction, with each layer Δz = 0.15m. The first layer z1 corresponds to the height with the largest number of point clouds, and its depth z1 = 0. The number of point clouds in this layer is h′(z1). The number of point clouds in the i-th layer is counted as h′(z i ), the number of each layered point cloud is normalized to obtain h(z i ), that is, h(z i )=h′(z i ) / h′(z1). This defines the contribution of the real backscatter signal sequence x(z) to the backscatter laser pulse point cloud sequence y(z) observed by the single-photon lidar in percentage form. Figure 3 The water depth z in the embodiment is shown i The corresponding impulse response function h(z i ). When the impulse response function H(z) is known, by deconvolution x(z) = H -1 (z)y(z) obtains the real underwater backscatter signal sequence x(z) after eliminating the after-pulse effect, where H -1 (z) represents the inverse matrix of H(z).

[0059] Preferably, the water body backscatter signal sequence x(z) and the water body parameter backscatter theoretical function can be expressed as:

[0060]

[0061] exp is the exponential function, sec is the secant trigonometric function; β(π) represents the volume scattering coefficient when the angle is π (backward), K d represents the water diffuse attenuation coefficient, which is the two parameters to be calculated; E t is the single pulse laser emission energy, η is the comprehensive efficiency of the single photon lidar system, h p is the Planck constant, which is approximately equal to 6.63×10-34 J·s, v is the laser wavelength efficiency, R is the flight altitude of the satellite-borne single-photon laser radar, and A r is the effective area of ​​the receiving aperture of the laser radar telescope. The above parameters are known parameters for a given laser radar system; the single pulse emission energy E t =93μJ, the overall efficiency of the single-photon lidar system η = 0.0615, and the laser wavelength efficiency v = 3e 8 / 532e -9 , the effective area of ​​the receiving aperture of the lidar telescope is A r =0.41m 2 , the flight altitude of the satellite-borne single-photon laser radar is R = 500km. w is the refractive index of water, usually 1.334; θ is the refraction angle of the emitted laser pulse after entering the water body. ICESat-2 / ATLAS is near-nadir incidence, θ≈0°; T a is the attenuation coefficient of the laser pulse when passing through the atmosphere, which can be obtained in real time by ICESat-2 / ATLAS. In the embodiment, T a =0.88; T w is the attenuation coefficient of the laser pulse entering the water body through the water surface. For the 532nm laser wavelength of ICESat-2 / ATLAS, T w ≈0.98. The above parameters are usually constant values ​​under given conditions or can be obtained by calculating the parameters of the lidar system;

[0062] In the above formula, there are two parameters to be substituted. The data in the backscatter signal sequence x(z) is usually greater than 10, and there is a data redundancy problem. In order to obtain more stable results and suppress the random errors in the data, the least squares fitting method is used to calculate the backscattering coefficient β(π) and the diffuse attenuation coefficient K d The data in the backscatter signal sequence x(z) is usually greater than 10 (i.e., n value). In the embodiment, the water depth range is 10m, the water depth direction interval Δz=0.15m, n=67, and there is a data redundancy problem. In order to obtain more stable results and suppress the random errors in the data, the least squares fitting method is used to calculate the backscatter coefficient β(π) and the diffuse attenuation coefficient K d When performing the least squares fitting, the backscattering coefficient β(π) and the diffuse attenuation coefficient K need to be givend The initial values ​​of the fitting can be selected from typical water optical parameters as the initial values. The 532nm wavelength light beam corresponds to the water absorption coefficient a, scattering coefficient b and backscattering coefficient b of different types of waters. b Table 1. Estimated diffuse attenuation coefficient K d :K d =a+4.18b b [1-0.52exp(-10.8a)].

[0063] Preferably, the water body backscattering coefficient in step S4 is a more commonly used water body optical parameter, and the water body backscattering coefficient b is calculated using the backscattering coefficient β(π) obtained in step 3. b , assuming that the backscattering is isotropic, the backscattering coefficient b b Defined as the sum of the volume scattering coefficients in the hemisphere space behind the incident light: b b =2πβ(π).

[0064] Table 1 Typical optical parameter values ​​of three types of water bodies (including water absorption coefficient, scattering coefficient, backscattering coefficient reference values)

[0065] Clean ocean water Coastal waters Turbid water in port <![CDATA[Absorption coefficient a (m -1 )]]> 0.114 0.179 0.366 <![CDATA[Scattering coefficient b (m -1 )]]> 0.037 0.398 1.824 <![CDATA[Backscattering coefficient b b (m -1 )]]> 0.00163 0.00052 0.0365

[0066] According to the above steps, in the embodiment, the backscatter coefficient b of the Indian Ocean waters to be measured is finally calculated. b =0.00247m -1 , diffuse attenuation coefficient K d =0.1047m -1 In fact, the backscatter coefficient b obtained by MODIS in this water area is b =0.00252m -1 , diffuse attenuation coefficient K d =0.1036m -1 , which is very consistent with the results calculated by the present invention, and also proves the feasibility and accuracy of the method of the present invention. Therefore, the present invention can quickly and accurately calculate the main optical parameters of the local water body through the point cloud data obtained by the satellite-borne single-photon laser radar when flying over the ocean, which serves as the basis for further analyzing various components and environmental parameters in the water body.

[0067] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, a person skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some or all of the technical features may be replaced by equivalents, and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for inverting water body optical parameters using a spaceborne single-photon laser radar, characterized in that: include: S1. Acquire the effective signal of the spaceborne single-photon laser radar when it flies above the water surface and performs measurements; S2. Eliminate the after-pulse effect of single-photon lidar using deconvolution algorithm; S3. Based on the theoretical model of water backscattering, the water backscattering coefficient and water diffuse attenuation coefficient are obtained by least square fitting; S4. Calculate the backscattering coefficient of water body; The effective signal is a water body backscattering point cloud, that is, a point cloud formed by a water body backscattering laser pulse. A water surface reflection point cloud with a higher spatial point cloud density is extracted by a spatial point cloud density classification method. The water surface reflection point cloud is a point cloud formed by a water surface reflecting a laser pulse. A point cloud with a higher height value than the water surface reflection point cloud corresponds to a point cloud formed by an atmospheric scattered laser pulse, and a point cloud with a lower height value than the water surface reflection point cloud corresponds to a water body backscattering point cloud. Step S1 includes the following sub-steps: S1.

1. Calculate the point cloud density within the unit space; S1.

2. Calculate the mean value μ and standard deviation σ within the spatial grid; S1.

3. Calculate the discrimination threshold TH between the water surface reflection point cloud and the atmosphere and water body backscattering point cloud; S1.

4. Calculate the average backscattered laser pulse point cloud sequence y(z) = [y(z1), y(z2), …, y(z n )], z represents the depth sequence, z1,z2,…,z n Represents different depths, and the statistical interval in the depth direction is Δz=z2-z1=z3-z2=z4-z3; Step S1.1 includes: the length interval of each spatial grid along the satellite flight direction is Δl, and the interval in the height direction is Δh. The total length of the data acquired by the satellite-borne single-photon lidar in the study area along the satellite flight direction is L, and the total length in the height direction is H. Then the number of grids along the satellite flight direction is n l and the number of grids in the height direction n h They can be expressed as: n l =ceil(L / Δl),n h =ceil(H / Δh), ceil represents the upward rounding function; the data in the study area are divided into spatial grids, and the spatial coordinate range F(i,j) of the i-th grid along the satellite flight direction and the j-th grid along the height direction satisfies: L start and L end are the starting coordinate and ending coordinate in the study area along the direction of satellite flight, H min and H max are the height direction recording range in the study area, that is, the minimum and maximum values ​​of the range window of the laser radar; the number of photon point clouds within each spatial grid coordinate range F(i,j) is counted, that is, the point cloud density within the unit space, recorded as N(i,j); In step S1.2, In step S1.3, TH = μ + e a ×σ,e a is the proportional coefficient. In the study area, all point clouds within the spatial range where the spatial grid density N(i, j) is greater than the discrimination threshold TH are classified as water surface reflection point clouds; In step S2, y(z) is the real water backscatter signal sequence x(z) = [x(z1), x(z2), ..., x(z n )] is the convolution of the impulse response function H(z) of the single-photon lidar system, that is, y(z)=H(z)*x(z), z represents the depth sequence, y(z)=H(z)*x(z) is equivalent to: z1,z2,…,z n Represents different depths; The pulse response function H(z) of the single-photon laser radar system can be obtained by calculating the received signal of the laser pulse incident on the planar target, or by calculating the received signal corresponding to a portion of the laser pulse energy intercepted inside the single-photon laser radar system. When the pulse response function H(z) is known, the deconvolution x(z) = H -1 (z)*y(z) to obtain the real underwater backscatter signal sequence x(z) after eliminating the after-pulse effect, where H -1 (z) represents the inverse matrix of H(z); The water body backscatter signal sequence x(z) and the water body parameter backscatter theoretical function can be expressed as: exp is the exponential function, sec is the secant trigonometric function; β(π) represents the volume scattering coefficient when the angle is π (backward), K d represents the water diffuse attenuation coefficient, which is the two parameters to be calculated; E t is the single pulse laser emission energy, η is the comprehensive efficiency of the single photon lidar system, h p Planck's constant is approximately equal to 6.63×10 -34 J·s, v is the laser wavelength efficiency, R is the flight altitude of the satellite-borne single-photon laser radar, A r is the effective area of ​​the receiving aperture of the laser radar telescope. The above parameters are known parameters for a given laser radar system; n w is the refractive index of water, θ is the refraction angle of the emitted laser pulse after entering the water, T a is the attenuation coefficient of the laser pulse when passing through the atmosphere, T w is the attenuation coefficient of the laser pulse when it passes through the water surface and enters the water body. The above parameters are usually constant values ​​under given conditions or can be obtained by calculating the laser radar system parameters; In the above formula, there are two parameters to be substituted. The data in the backscatter signal sequence x(z) is usually greater than 10, and there is a data redundancy problem. In order to obtain more stable results and suppress the random errors in the data, the least squares fitting method is used to calculate the backscattering coefficient β(π) and the diffuse attenuation coefficient K d ; The water body backscattering coefficient in step S4 is a more commonly used water body optical parameter. The water body backscattering coefficient b is calculated using the backscattering coefficient β(π) obtained in step 3. b , assuming that the backscattering is isotropic, the backscattering coefficient b b Defined as the sum of the volume scattering coefficients in the hemisphere space behind the incident light: b b =2πβ(π).

Images