Ocean profile inversion method based on spread matching for spaceborne photon counting lidar

Abstract

The ocean profile inversion method of spaceborne photon counting lidar based on broadening matching includes land / ocean photon point cloud quality control; photon point cloud elevation detrending; generating system pulse response function; establishing water profile signal; post-pulse correction; ocean optical parameter inversion, and converting the corrected water profile signal into S c ( z ) is combined with the lidar equation, and the water diffuse attenuation coefficient is calculated based on the consistency assumption of the Collis method that the optical parameters of the subsurface ocean do not change significantly with depth. k d Inversion. Aiming at the problem of post-pulse in the ATLAS system, the present invention provides a post-pulse correction method based on the broadening matching theory. Under the premise of fully considering the effect of sea surface on signal broadening, by matching the system pulse response function with the same pulse width, the problem that the traditional post-pulse correction method is only applicable to some sea areas is solved, and the robustness of the ocean data processing method of space-borne photon counting lidar is effectively improved.

Description

Technical Field

[0001] The invention relates to a satellite-borne photon counting laser radar ocean data processing method based on signal broadening matching, belonging to the technical field of satellite-borne laser radar. Background Art

[0002] The study of changes in the optical properties of ocean water and their influencing mechanisms is of great significance for establishing high-precision water color remote sensing inversion methods, using optical means to monitor the dynamic changes of suspended particles, predicting changes in underwater visibility, and using optical means to study the impact of short-time-scale ocean and weather processes on water bodies. Satellite-borne lidar can achieve all-day, high-resolution global observations, providing a large amount of global observation data for multiple research fields such as remote sensing mapping and marine science. In order to monitor changes in sea and land ice elevation and measure the height of global vegetation canopy, NASA launched the first satellite ICESat-2 equipped with the photon counting lidar ATLAS system in November 2018. The ATLAS transmitting system adopts a design of 532nm micro-pulses, 10KHz high repetition rate, and three strong and three weak six beams. The receiving system uses a high-sensitivity photomultiplier tube as a single-photon detector to achieve the acquisition of high-temporal and spatial resolution global observation point cloud data. Since 532nm wavelength laser can penetrate the subsurface of seawater, the ATLAS system can also obtain water profile signals. Its receiving system has a maximum resolution of 50ps, which makes it possible to invert vertical high-resolution water optical parameters. However, the water profile signals directly obtained by the ATLAS system have obvious after-pulse phenomena, which restricts the application of observation data.

[0003] In current research, post-pulse correction is mainly performed by deconvolution. The correction method first assumes that the water body signal and the land signal are subjected to the same modulation effect, uses desert data to establish a system impulse response function, and performs post-pulse correction on the water body signal. The corrected water body profile signal can be used for ocean optical parameter inversion, and the relative difference between the inversion result and the MODIS data can be reduced to about 10%.

[0004] Since the algorithm needs to be based on the assumption of the same modulation of land and sea, good correction and inversion results can only be obtained in some sea areas when the algorithm is applied. The stability and robustness of the algorithm when applied in different sea areas around the world are still insufficient, and the correction results are sensitive to changes in sea surface wind speed. The algorithm needs to ensure that more accurate post-pulse correction can be performed in most sea areas to achieve high-precision ocean optical parameter inversion using space-borne lidar data on a global scale. Summary of the invention

[0005] The purpose of the present invention is to provide a method for inverting ocean profile of a space-borne photon counting lidar based on broadening matching. Based on the broadening effect of the land surface and the ocean surface on the signal, a system pulse response function with the same broadening is matched to the water profile signal under sea surface broadening, thereby solving the problem of insufficient robustness in traditional post-pulse correction methods. This method can effectively reduce the influence of post-pulses on the signal of the space-borne lidar ATLAS system, and obtain high-precision water profile signals and ocean optical parameters.

[0006] The present invention comprises the following steps:

[0007] Step 1, land / ocean photon point cloud quality control:

[0008] Use ATL09 cloud data and satellite optical images (such as Sentinel-2 optical images) that match the ATL03 data in time and space to determine the locations of clouds and other non-target types of surfaces on the surface or ocean in the ATL03 data, and remove the photon point cloud in the area according to the corresponding time and geographical location;

[0009] Step 2, photon point cloud elevation detrending:

[0010] The cumulative histogram of photon events from the previous step is obtained by accumulating along the track. The elevation position of the maximum number of photon events within the cumulative number is used as the surface or sea surface elevation for detrending. The elevation values ​​of all photon events within the cumulative number are subtracted from the surface or sea surface elevation to eliminate the influence of the surface or sea surface fluctuations.

[0011] Step 3: Generate the system impulse response function:

[0012] Calculate the surface slope and surface roughness along the track direction, and then obtain the surface width. Use the surface width as the classification label of the photon point cloud, establish a land photon point cloud library with different widths, and obtain the system pulse response function with different pulse widths through the accumulation method;

[0013] Step 4: Establish water profile signal:

[0014] The water body point cloud after elevation detrending is accumulated along the track, and the accumulated result is the water body profile signal;

[0015] Step 5, post pulse correction:

[0016] The sea surface slope and sea surface roughness along the track are calculated according to the sea surface elevation value in the water profile signal, and the sea surface broadening of the current water profile signal is obtained accordingly. The system pulse response function with the same pulse width is matched for deconvolution, and the pulse is corrected to obtain the real water profile signal. S c ( z );

[0017] Step 6, inversion of ocean optical parameters:

[0018] The corrected water profile signal S c ( z ) is combined with the lidar equation, and the water diffuse attenuation coefficient is calculated based on the consistency assumption of the Collis method that the optical parameters of the subsurface ocean do not change significantly with depth. k d Inversion.

[0019] In step 1, the water body photon point cloud is quality controlled. The quality control method is similar to that of the land photon point cloud, mainly cloud quality control and abnormal window quality control.

[0020] To eliminate the photon point clouds in the locations where there are clouds or other non-target types of surface above the ground or water, the ATL09 data is used as auxiliary data, the cloud tag positions of the ATL09 data are read, and the photon point clouds in the cloud-covered locations in the ATL03 photon point cloud are eliminated.

[0021] Photon point clouds in clouds and other non-target types of surface locations above water bodies are eliminated, and an abnormal windowing quality control method is adopted. The abnormal windowing quality control reads the maximum and minimum positions of the photon event ellipsoid height, calculates the elevation difference, sets a threshold based on the general windowing conditions in the sea area, and eliminates photon point clouds with elevation differences greater than the threshold.

[0022] In the detrending of the photon point cloud elevation in step 2, the specific processing flow of the photon point cloud is as follows:

[0023] First, the photon point cloud after quality control is divided into blocks according to a certain distance along the track and vertical height. The block division standard is to meet the minimum cumulative number of times that can identify the surface or sea surface, and the noise rate P n The number of photon events after vertical height binning is determined by the noise rate of ATL03. P s , after block division, the photon events in the small blocks are accumulated to obtain the photon event accumulation histogram;

[0024] Then, the elevation of the maximum position of the number of photon events in the cumulative histogram is used as the surface time slice, and the elevation value of the surface time slice is used as the surface elevation;

[0025] Finally, the photon point cloud elevation value within the current distance along the track is subtracted from the surface elevation value, the elevation of the surface position is set to 0, and the photon elevations of the remaining positions are set to the relative surface elevation after detrending. This can effectively reduce the impact of the difference in ground photon flight time on the system impulse response function when the surface elevation is inconsistent, and the algorithm data processing speed is faster.

[0026] The minimum cumulative number of times that can identify the land surface or the sea surface in step 2 is as follows:

[0027] Cumulative number of times at surface or sea position n The following standards must be met, where P s , P n are the detection probabilities of signal and noise photon events, respectively;

[0028] (1) .

[0029] The system impulse response function library is generated as described in step 3, as follows:

[0030] To calculate the surface slope and surface roughness in the along-track direction, first, the photon point cloud in a certain along-track direction is grouped into a group. A group contains several photon point cloud groups after elevation detrending in step 2. It is assumed that the broadening caused by the surface in a group is the same. When calculating the slope, the arc tangent value of the slope of the two elevation points before and after the target elevation point is used as the surface slope of the target point, and the mean of the slope in a group is used as the slope of the photon point cloud group; when calculating the roughness, the standard deviation of the surface elevation in a group is used as the surface roughness, and the surface broadening is calculated. Finally, the surface broadening is used as the classification label of the photon point cloud, and a land photon point cloud library with different broadenings is established. The noise rate is calculated based on the number of photon events above the surface and the cumulative length. The system pulse response function with different pulse widths is established using the noise-deducted photon point cloud library.

[0031] Step 4 describes the establishment of a water profile signal:

[0032] The detrending results of the photon point cloud that meet the consistency assumption of the optical properties of water bodies are accumulated, and the accumulated results of the underwater part are recorded as the water body profile signal before denoising. The average number of photon events above the sea surface is counted, and the noise rate is calculated considering the cumulative length along the track and the width of the time slice in the vertical direction. The number of photon events at each depth of the water body profile signal before denoising is subtracted from the average number of noise photon events to obtain the water body profile signal.

[0033] Aiming at the problem of post-pulse in the ATLAS system, the present invention provides a post-pulse correction method based on the broadening matching theory. Under the premise of fully considering the effect of sea surface on signal broadening, by matching the system pulse response function with the same pulse width, the problem that the traditional post-pulse correction method is only applicable to part of the sea area is solved, and the robustness of the ocean data processing method of space-borne photon counting lidar is effectively improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 is the expected effect of different surface roughness and slope on pulse width.

[0035] in, Figure 1 (a) is the surface pulse width under the corresponding conditions. Figure 1 (b) shows the sea surface pulse width under the corresponding conditions, and the contour lines represent the expected RMS pulse width of the received signal.

[0036] Figure 2 This is a comparison of the quality control of the land photon point cloud before and after.

[0037] in, Figure 2 (a) is the photon point cloud before quality control. Figure 2 (b) in the middle is the photon point cloud after quality control, and the box indicates the abnormal surface or the area with clouds above the surface.

[0038] Figure 3 It is the classification of land photon point cloud and the system impulse response function diagram.

[0039] in, Figure 3 (a) is the initial photon point cloud, and the surface broadening is marked with different grayscales; Figure 3 (b.1)-(b.3) are the local enlarged images of the photon point cloud at the positions of case1-case3 in (a); Figure 3 In the figure, (c.1)-(c.3) are the corresponding broadened system impulse response functions.

[0040] Figure 4 It is the water profile signal and deconvolution result diagram.

[0041] in, Figure 4 (a) is the water profile signal before correction. Figure 4 (b) is the corrected water profile signal.

[0042] Figure 5 This is a comparison chart of the inversion results of different correction methods.

[0043] in, Figure 5 (a)-(c) are the inversions of the Black Sea, Atlantic Ocean and Pacific Ocean respectively. k d result.

[0044] Figure 6 This is a statistical diagram of inversion errors in different sea areas.

[0045] Figure 7 Flowchart of the present invention. DETAILED DESCRIPTION

[0046] The target backscattered signal received by the spaceborne laser radar will be modulated by the atmosphere, target surface and receiving system during the transmission process. g ( t ) can be expressed as:

[0047] (2) .

[0048] In the formula, I ( t ) represents the target backscatter signal, h atm ( t ) represents the modulation effect of the atmosphere on the signal during transmission. h surf ( t ) represents the modulation effect of the target surface on the signal, h Rec ( t ) represents the modulation effect of the receiving system on the signal. Among them, the modulation effect of the atmosphere can cause the signal to attenuate and broaden, which is a nonlinear modulation, and it is necessary to accurately obtain the modulation function under different atmospheric conditions. However, due to the limitation of the data transmission bandwidth of the ATLAS system, the atmospheric modulation function cannot be accurately obtained, and it is necessary to select data under clear sky conditions to weaken the modulation effect of the atmosphere as much as possible.

[0049] In general, h surf ( t ) belongs to nonlinear modulation. Different sea surfaces have different modulation effects on the signal, but the modulation effect of the receiving system can be obtained through analysis. If the two modulation effects are regarded as a whole, equation (2) can be expressed as:

[0050] (3).

[0051] In the formula, h ’ Rec ( t ) represents the modulation function of the receiving system with the sea surface modulation effect superimposed. At this time, since the receiving system is a linear system, when the modulation function h ’ Rec ( t ) is represented by the receiving system modulation function with superimposed target surface modulation, the water body signal I ocean ( t ) can be obtained by using the deconvolution method.

[0052] According to Gardner's theory, the broadening of the reflected signal received by the space-borne laser radar from the ground or sea surface is different, which is mainly affected by the surface roughness and surface slope. The expected root mean square pulse width of the received signal on the ground and sea surface is expressed as:

[0053] (4),

[0054] (5).

[0055] In the formula, σ f is the RMS pulse width of the emitted light, σ h is the RMS pulse width of the impulse response, c is the speed of light, z is the satellite altitude, θ T is the laser divergence angle, ϕ is the nadir angle, θ is the surface slope, σ ξ is the root mean square of surface roughness, λ 3 is the skewness coefficient, S is the sea surface slope. For the ATLAS system, σ f , σ h , θ T , ϕ Determined by the system, z It is determined when the data location is selected. c is a constant, so the expected surface pulse width is mainly affected by the surface roughness Var( ξ ) and surface slope θ The expected sea surface pulse width is mainly affected by the sea surface roughness. σ ξ and sea surface slope S The expected effect of different surface roughness and slope on pulse width is as follows: Figure 1 shown.

[0056] Depend on Figure 1 It can be seen that under certain surface roughness and slope conditions, the land signal broadening may be the same as the water signal broadening, and different surface roughness and slope may also cause the same signal broadening. When the water signal broadening is known, the land signal broadening can be made consistent with the water signal broadening by limiting the surface roughness and slope. The surface conditions and signal broadening are not in a one-to-one correspondence. Photon point clouds with different surface roughness and slope but the same broadening can be accumulated to obtain the corresponding broadened system impulse response function.

[0057] In view of the above situation, the present invention provides a method for inverting ocean profiles using a spaceborne photon counting lidar based on broadening matching, such as Figure 7 As shown, the following steps are included.

[0058] Step 1, land / ocean photon point cloud quality control:

[0059] Step 1.1, the ATL09 data product provides cloud layer tags and gives the time of each cloud tag position. Use the cloud layer tags as ATL03 cloud tags and record the corresponding cloud layer position at that time. Find the corresponding cloud layer position in the ATL03 data and remove the photon point cloud in the ATL03 data that matches the cloud tag data time.

[0060] Step 1.2, the land study area used to establish the system impulse response function in the present invention may have abnormal areas such as surface water and surface cover. Sentinel-2 optical image data can effectively identify the location of surface water and cover. Select the Sentinel-2 optical image that matches the ATL03 data in time and space, identify and read the longitude and latitude range of the abnormal area, and mark it as the area to be eliminated. Mark the longitude and latitude range in the ATL03 data that matches the data in the area to be eliminated, and eliminate the photon point cloud in the area.

[0061] Step 2: Detrend the land photon point cloud elevation based on the quality control results of step 1:

[0062] Before accumulating the system impulse response function, the quality-controlled photon point cloud needs to be aligned with the ground surface and detrended in elevation. First, the quality-controlled photon point cloud is divided into blocks according to a certain along-track distance and vertical height, and the photon events in the small blocks are accumulated to obtain the cumulative histogram of photon events. The block results in the along-track direction and the vertical direction must meet the minimum accumulation times shown in formula (1). Then, the elevation of the position with the maximum number of photon events in the cumulative histogram is used as the ground surface elevation. Secondly, the elevation value of the photon point cloud within the current along-track distance is subtracted from the ground surface elevation value, and the elevation of the ground surface position is set to 0. The photon elevations of the remaining positions are all set to the relative ground surface elevation after detrending.

[0063] Step 3: Calculate the widening caused by the surface within a certain distance along the track, and establish the system impulse response function after classification:

[0064] To calculate the surface slope and surface roughness in the along-track direction, first, the photon point clouds in a certain along-track direction are grouped together. A group contains several photon point cloud groups after elevation detrending in step 2. It is assumed that the widening caused by the surface in a group is the same. When calculating the slope, the arc tangent value of the slopes of the two elevation points before and after the target elevation point is used as the surface slope of the target point, and the mean of the slopes in a group is used as the slope of the photon point cloud group; when calculating the roughness, the standard deviation of the surface elevation in a group is used as the surface roughness, and the surface slope and surface roughness are substituted into formula (3) to calculate the surface widening. Finally, the surface widening is used as the classification label of the photon point cloud, and a land photon point cloud library with different widenings is established. The photon point cloud library is used to establish a system pulse response function with different pulse widths. The classification of photon point clouds with different pulse widths is as follows: Figure 3 shown.

[0065] Step 4: Establish water profile signal:

[0066] Step 4.1, quality control of water body photon point cloud. The quality control method is similar to that of land photon point cloud, mainly cloud quality control and abnormal window quality control. When controlling cloud quality, ATL09 data is still used as auxiliary data, cloud tag positions are read, and photon point clouds at cloud cover positions in ATL03 photon point cloud are eliminated. When controlling abnormal window quality, the maximum and minimum positions of photon event ellipsoid height are read, the elevation difference is calculated, and the threshold is set according to the general window conditions in the sea area, and photon point clouds with elevation differences greater than the threshold are eliminated.

[0067] Step 4.2, detrending water elevation: First, the cumulative length of the water profile signal is set based on the assumption that the optical properties of the water body are consistent. Secondly, the photon point cloud within the cumulative length is divided into several time slices, and the photon point cloud within the time slice is accumulated in blocks according to a certain vertical resolution, and the maximum position in the photon event accumulation histogram is used as the sea surface elevation. Then, the photon event elevation within the time slice is subtracted from the sea surface elevation to align with the sea surface position.

[0068] Step 4.3, establish water profile signal: accumulate the detrending results of the photon point cloud within the range of the water body optical property consistency assumption, and record the accumulated results of the underwater part as the water body profile signal before denoising. Count the average number of photon events above the sea surface, and calculate the noise rate considering the cumulative length along the track and the vertical time slice width. Difference between the number of photon events at each depth of the water body profile signal before denoising and the average number of noise photon events to obtain the water body profile signal.

[0069] Step 5, post pulse correction and application:

[0070] Based on the water body elevation information in step 4.2, calculate the sea surface slope and sea surface roughness within the cumulative length of the water body profile signal, and substitute them into formula (4) to obtain the pulse width of the water body profile signal. Select the system pulse response function with the same pulse width as the water body profile signal, normalize the system pulse response function, and use the deconvolution method to perform post-pulse correction on the water body profile signal through the system pulse response function to obtain the real water body profile signal. Combine the real water body profile signal with the lidar equation to calculate the target ocean optical parameters. The water body profile signal and deconvolution results are shown in Figure 4 .

[0071] Step 6, inversion of ocean optical parameters:

[0072] The Collis method makes certain assumptions about the variation of water optical parameters with depth, as shown in formula (6). The logarithmic result of the processed water profile signal is fitted once, and -1 / 2 of the fitting result can be used as the water diffuse attenuation coefficient.

[0073] (6).

[0074] Example verification.

[0075] This experimental data uses desert and salt marsh data as the land research area, and uses ATL09 and Sentinel-2 optical image data to perform quality control on the land photon point cloud. Taking the salt marsh data as an example, the quality control results are as follows Figure 2 As shown. Accumulate the photon point clouds after quality control in several consecutive pulses. The cumulative length along the track must be greater than 7m, and the vertical resolution must be greater than 0.03m. Obtain the surface position, perform elevation detrending, and calculate the surface roughness and surface slope through the surface elevation within the cumulative length, and then calculate the pulse width. Classify the photon point clouds according to the pulse width, and establish a photon point cloud library. The photon point cloud library contains photon point clouds corresponding to multiple different pulse widths with a pulse width greater than 1.5ns. The photon point clouds within the same pulse width are accumulated as the system pulse response function under the current pulse width for subsequent post-pulse correction. Taking the desert data as an example, the photon point cloud classification results and their corresponding system pulse response functions with pulse widths of 1.6ns, 2.0ns and 2.8ns are shown in Figure 3 shown.

[0076] The post-pulse correction method based on pulse broadening matching was experimentally tested in the Black Sea (33°N, 43°E), the Atlantic Ocean (56°S, 28°E) and the Pacific Ocean (16°S, 43°E). First, the ATL09 data was used for cloud quality control. The photon point clouds that met the requirements of several pulses were accumulated. The accumulated length along the track was required to be greater than 7m, and the vertical resolution was required to be greater than 0.03m to obtain the sea surface position. After the elevation was detrended according to the sea surface position, the photon point clouds within the length that met the consistency assumption of the ocean optical properties were accumulated. The accumulated length along the track was generally not more than 20km to obtain the water profile signal. The sea surface in the Black Sea is relatively flat, with a total of 117 tracks of data, and a total of 1110 water profile signals available after quality control. The western Atlantic Ocean is located near the ocean current, with a single sea surface feature and stable sea conditions, with a total of 233 tracks of data, and a total of 1657 water profile signals available after quality control; the Pacific Ocean has strong winds and high waves, and the widening caused by the sea surface is the largest, with a total of 139 tracks of data, and a total of 965 water profile signals available after quality control. The sea surface roughness and sea surface slope are calculated by accumulating the sea surface elevation values ​​within the length, and then the pulse width is calculated. Match the system pulse response function that is the same as the sea surface widening, deconvolve the water profile signal, and use the corrected pulse to obtain the real water profile signal. Taking the Black Sea data as an example, the system pulse response function with a pulse width of 2.0ns, the water profile signal and its subsequent pulse correction results are shown in Figure 4Substitute the real water profile signal into the inversion formula to invert the water diffuse attenuation coefficient k d , and compared with MODIS data to evaluate the data processing results.

[0077] Figure 5 This is a comparison chart of the inversion results using the traditional post-pulse correction method and the post-pulse correction method based on signal broadening matching. MODIS data is selected as the reference truth value. The horizontal axis in the figure represents the along-track distance, and the vertical axis represents the inversion result. The dotted line represents the reference truth value provided by MODIS data, and the triangle line represents the inversion result obtained using the broadening matching method. k d The square line represents the inversion result obtained using the traditional method. k d As can be seen from the figure, in the data of the Black Sea, the Atlantic Ocean and the Pacific Ocean, the ocean optical parameters inverted using the broadening matching method are less different from the MODIS data.

[0078] Figure 6 The figure shows the statistical results of the inversion errors of data from multiple sea areas. The horizontal axis represents different sea areas, and the vertical axis represents the relative error distribution between the inversion results and MODIS. The dark box plot represents the inversion results using the widening matching method, and the light box plot represents the inversion results using the traditional method. The middle line of the box plot is the median value, the upper limit of the dotted line is the maximum value, the lower limit is the minimum value, the top of the box is the 75th percentile, and the bottom of the box is the 25th percentile. As can be seen from the figure, the inversion results of the widening matching method are closer to the reference true value provided by MODIS.

[0079] Based on the above verification results, it can be concluded that the method provided by the present invention can effectively reduce the correction error and improve the inversion accuracy of ocean optical parameters compared with traditional data processing methods.

[0080] The specific embodiments described herein are merely examples of the present invention. A person skilled in the art may make various modifications or additions to the specific implementation process described, or replace the specific implementation process with a similar method, but this will not deviate from the spirit of the present invention or exceed the scope defined by the appended claims.

Claims

1. A method for inverting ocean profiles using a spaceborne photon counting lidar based on broadening matching, comprising the following steps: Step 1, land / ocean photon point cloud quality control: Using ATL09 cloud data and satellite optical images that match the ATL03 data in time and space, determine the location of clouds and other non-target types of surfaces above the surface or water bodies in the ATL03 data, and remove the photon point cloud in the area according to the corresponding time and geographical location; Step 2, photon point cloud elevation detrending: The cumulative histogram of photon events from the previous step is obtained by accumulating along the track. The elevation position of the maximum number of photon events within the cumulative number is used as the surface or sea surface elevation for detrending. The elevation values ​​of all photon events within the cumulative number are subtracted from the surface or sea surface elevation to eliminate the influence of the surface or sea surface fluctuations. Step 3: Generate the system impulse response function: Calculate the surface slope and surface roughness along the track direction, and then obtain the surface width. Use the surface width as the classification label of the photon point cloud, establish a land photon point cloud library with different widths, and obtain the system pulse response function with different pulse widths through the accumulation method; Step 4: Establish water profile signal: The water body point cloud after elevation detrending is accumulated along the track, and the accumulated result is the water body profile signal; Step 5, post pulse correction: The sea surface slope and sea surface roughness along the track are calculated according to the sea surface elevation value in the water profile signal, and the sea surface broadening of the current water profile signal is obtained accordingly. The system pulse response function with the same pulse width is matched for deconvolution, and the pulse is corrected to obtain the real water profile signal. S c ( z ); Step 6, inversion of ocean optical parameters: The corrected water profile signal S c ( z ) is combined with the lidar equation, and the water diffuse attenuation coefficient is calculated based on the consistency assumption of the Collis method that the optical parameters of the subsurface ocean do not change significantly with depth. k d Inversion; In the detrending of the photon point cloud elevation in step 2, the specific processing flow of the photon point cloud is as follows: First, the photon point cloud after quality control is divided into blocks according to a certain distance along the track and vertical height. The block division standard is to meet the minimum cumulative number of times that can identify the surface or sea surface, and the noise rate P n The number of photon events after vertical height binning is determined by the noise rate of ATL03. P s , after block division, the photon events in the small blocks are accumulated to obtain the photon event accumulation histogram; Then, the elevation of the maximum position of the number of photon events in the cumulative histogram is used as the surface time slice, and the elevation value of the surface time slice is used as the surface elevation; Finally, the photon point cloud elevation value within the current along-track distance is subtracted from the surface elevation value, the elevation of the surface position is set to 0, and the photon elevations of the remaining positions are set to the detrended relative surface elevation.

2. The method as described in claim 1 is characterized in that in step 1, the photon point clouds in the clouds and other non-target surface positions above the surface or water bodies are eliminated by using ATL09 data as auxiliary data, reading the cloud tag position of the ATL09 data, and eliminating the photon point clouds in the cloud-covered positions in the ATL03 photon point cloud.

3. The method as described in claim 1 is characterized in that in step 1, the photon point clouds in clouds and other non-target types of surface positions above the water body are eliminated, and an abnormal windowing quality control method is adopted. The abnormal windowing quality control reads the maximum and minimum positions of the photon event ellipsoid height, calculates the elevation difference, sets a threshold according to the general windowing conditions in the sea area, and eliminates photon point clouds with elevation differences greater than the threshold.

4. The method according to claim 1, characterized in that The minimum cumulative number of times that can identify the land surface or the sea surface in step 2 is as follows: Cumulative number of times at surface or sea position n The following standards must be met, where P s , P n are the detection probabilities of signal and noise photon events, respectively; 。 5. The method according to claim 1, wherein the generating system impulse response function in step 3 is as follows: To calculate the surface slope and surface roughness in the along-track direction, first, the photon point clouds in a certain along-track direction are grouped together, and each group contains several photon point cloud groups after elevation detrending in step 2. It is assumed that the broadening caused by the surface in a group is the same; when calculating the slope, the arc tangent value of the slopes of the two elevation points before and after the target elevation point is taken as the surface slope of the target point, and the mean value of the slope in a group is taken as the slope of the photon point cloud group; when calculating the roughness, the standard deviation of the surface elevation in a group is taken as the surface roughness, and the surface broadening is calculated; finally, the surface broadening is used as the classification label of the photon point cloud, and a land photon point cloud library with different broadenings is established. The noise rate is calculated based on the number of photon events above the surface and the cumulative length, and the system pulse response function with different pulse widths is established using the noise-deducted photon point cloud library.

6. The method according to claim 1, characterized in that the step of establishing the water body profile signal in step 4 is as follows: The detrending results of the photon point cloud that meet the consistency assumption of the optical properties of water bodies are accumulated, and the accumulated results of the underwater part are recorded as the water body profile signal before denoising. The average number of photon events above the sea surface is counted, and the noise rate is calculated considering the cumulative length along the track and the width of the time slice in the vertical direction. The number of photon events at each depth of the water body profile signal before denoising is subtracted from the average number of noise photon events to obtain the water body profile signal.

Images