A method for optical stereoscopic image information assisting in denoising of spaceborne laser altimetry data

Abstract

The application relates to a method for removing noise from spaceborne laser altimetry data by using optical stereo image information. Firstly, dense matching point clouds are generated by using optical stereo image dense matching, and noise contained in the dense matching point clouds is removed; then, according to the along-track distance information, longitude and latitude information and height information of the laser altimetry data, along-track terrain slope and along-track terrain height information are calculated by using the dense point cloud data; then, coarse registration is performed by using the optimal terrain height extracted from the laser altimetry data and the along-track terrain height data extracted from the dense matching point clouds, and the dense matching point cloud coordinates are updated by using the coarse registration parameters; finally, the slope and height information extracted from each data segment are used to improve the density-based clustering algorithm, so that the noise points in the laser altimetry data are effectively removed.

Description

Technical Field

[0001] This invention relates to a method for the composite processing of spaceborne laser altimetry data and optical stereo images, belonging to the field of remote sensing technology; in particular, it relates to a method for noise removal from spaceborne laser altimetry data using satellite optical stereo image information. Background Technology

[0002] With the launch and use of various photon-counting laser altimetry satellites such as ICESAT-2 / ATLAS, high-precision three-dimensional coordinate data of the global Earth's surface can be obtained quickly.

[0003] However, due to the influence of various factors such as solar background noise, laser altimetry data contains a large amount of noise, especially daytime laser altimetry data, which greatly interferes with the extraction of three-dimensional surface information from laser altimetry data, and has a more serious impact on the information extraction of weak beam laser altimetry data.

[0004] Optical stereo imaging can enable rapid localization of targets of interest. Although the localization accuracy is low when using the original orientation parameters of optical images, the relative accuracy of the localization data is relatively stable.

[0005] Therefore, the surface slope and elevation information estimated by optical stereo images can be used to assist in noise removal from laser altimetry data. Summary of the Invention

[0006] The purpose of this invention is to provide a method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information, in order to solve the problem of signal photon extraction from laser altimetry data.

[0007] To achieve the above objectives, the present invention includes:

[0008] The technical solution of the present invention, a method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information, includes the following steps:

[0009] 1) Divide the laser altimetry data into several segments according to a certain distance along the track, and fit the centerline of the corresponding data segment based on the laser altimetry data in each data segment;

[0010] 2) Establish a buffer zone of a certain width with the center line as the center, extract the three-dimensional coordinate data of the densely matched point cloud in the buffer zone, and use the densely matched point cloud data in this buffer zone to fit the terrain plane.

[0011] The densely matched point cloud data is generated by densely matching the original orientation parameters of the optical stereo image of the same area.

[0012] 3) Calculate the terrain slope along the track of the terrain profile line where the corresponding data segment is located; the terrain profile line is the intersection line of the vertical plane where the center line of the corresponding data segment is located and the terrain plane obtained by fitting the densely matched point cloud.

[0013] 4) Using a density-based clustering algorithm, the laser altimeter cloud data within each data segment is iterated and the density value of each point based on the optimal search direction of the corresponding data segment is calculated. Points with density values ​​less than the classification density threshold within the corresponding data segment are designated as noise points. The optimal search direction is determined by using the terrain slope along the track of the data segment corresponding to the point to be determined, as well as the terrain slope along the track of several data segments before and after it, as the neighborhood search direction, and then adding a horizontal search direction. The maximum density value calculated by each search direction is taken as the density value of that point, and the search direction corresponding to the maximum density value is the optimal search direction for that point. The classification density threshold is calculated based on the density values ​​of all points within the corresponding data segment.

[0014] Laser altimetry technology boasts high precision, geographical independence, and all-weather operation, enabling rapid acquisition of high-precision three-dimensional coordinates of the Earth's surface globally. However, it also suffers from drawbacks such as susceptibility to solar background noise, the ability to acquire only topographic profile data, lack of texture features, and difficulty in extracting surface signal points. This invention proposes a method for denoising spaceborne laser altimetry data using optical stereo imagery. This method fully utilizes the slope and elevation information from optical stereo imagery to effectively eliminate noise in laser altimetry data, achieving a superior denoising effect.

[0015] Specifically, considering that the dense matching point cloud of optical stereo imagery has good relative positioning accuracy and can estimate the terrain slope and elevation more accurately, this invention uses a segmented fitting method on the surface profile based on laser altimetry data to map the corresponding data segments to the laser altimetry data based on the track-side terrain slope and elevation obtained from the optical stereo imagery. Noise identification is then performed on the laser altimetry data based on point cloud density, terrain slope, and elevation, thus realizing noise removal of spaceborne laser altimetry data assisted by optical stereo imagery information.

[0016] Furthermore, step 3) further calculates the terrain elevation along the track of the topographic profile line where each data segment is located. The terrain elevation along the track is the elevation calculated at the center point of the corresponding data segment's terrain profile line. The best-fit slope and best terrain elevation for the corresponding data segment's laser altimetry data are then calculated using the terrain slope along the track. Before step 4), the following coarse registration steps are performed between the laser altimetry data and the densely matched point cloud data:

[0017] Extract the track-side terrain elevation and optimal terrain elevation for each data segment, calculate the standard deviation of the elevation difference for each data segment, and remove elevation values ​​exceeding the set standard deviation. The elevation difference is the difference between the optimal terrain elevation estimated by the laser altimetry data for the corresponding data segment and the track-side terrain elevation estimated by the densely matched point cloud data. After removing the excess elevation values, fit the elevation values ​​of the densely matched point cloud data for each data segment into a curve, and use interpolation to ensure that there is a corresponding elevation value for each data segment. Translate the elevation fitting curve of the matched point cloud according to a set step size, and calculate the standard deviation of the elevation difference for each data segment after each translation. Take the translation distance with the smallest standard deviation of elevation difference as the coarse registration parameter.

[0018] After correcting the densely matched point cloud using the coarse registration parameters, the terrain slope along the track is refitted for use in the density-based clustering algorithm in step 4). Before being used to correct the elevation of the densely matched point cloud, the coarse registration parameters need to be converted from track distance and elevation to three-dimensional coordinates.

[0019] Optical stereo imagery inherently suffers from low original positioning accuracy, and the accuracy of corresponding densely matched point clouds is also limited. The slope of the terrain along the track obtained from fitting the densely matched point clouds for each data segment does not accurately correspond to the potential signal points of the laser altimetry data for the same data segment in the actual terrain, exhibiting certain deviations in both the planar and elevation directions. The best-fit slope data for each data segment obtained from the terrain slope data along the track is not accurate. Considering that the elevation and slope information estimated by both laser altimetry data and optical stereo imagery can reflect the actual terrain well, coarse registration of laser altimetry data and densely matched point clouds based on the principle of minimizing the standard deviation of elevation difference can effectively improve the accuracy of slope and elevation estimation for each data segment of laser altimetry data, thereby improving the denoising quality of the laser altimetry data.

[0020] Furthermore, due to noise interference, a few estimated optimal terrain elevations may be located in noisy areas with high density, constituting pseudo-elevation data. These pseudo-elevation data typically have either high or low elevation values. Therefore, abrupt pseudo-elevation data was removed before coarse registration, improving the effectiveness of the coarse registration. The pseudo-elevation data removal method involves removing the optimal terrain elevations and track-side terrain elevations of the corresponding data segments from all data segments where the elevation difference exceeds three times the standard deviation.

[0021] According to previous analysis, the terrain obtained by dense matching point cloud fitting and the potential signal points of laser altimetry point cloud have a certain offset in spatial distribution. Therefore, when using the track-side terrain slope fitted by dense matching point cloud fitting to obtain the best fitting slope for the potential signal points of laser altimetry point cloud in each data segment, in addition to considering the track-side terrain slope of the current data segment, the track-side terrain slopes of several data segments before and after should also be added.

[0022] Further, in step 3), the optimal fitting slope and optimal terrain elevation of the corresponding laser altimetry data segment are calculated using the track-side terrain slope as follows:

[0023] For each data segment, several vertical partitioning windows are established in the elevation direction with a set step size. The length of the window is equal to the length of the data segment, and the height is a preset window height. The set step size is less than the window height to ensure a certain degree of overlap between windows. The track-side terrain slope of the corresponding data segment, the data segment before the corresponding data segment, and the data segment after the corresponding data segment are used as the rotation direction of the vertical partitioning windows. The laser altimeter cloud density of each vertical partitioning window is calculated. The track-side terrain slope corresponding to the vertical partitioning window with the highest laser altimeter cloud density is taken as the best-fit slope of the corresponding data segment. The elevation of the center point of the vertical partitioning window with the highest laser altimeter cloud density is taken as the optimal terrain elevation.

[0024] The higher the density of laser altimetry cloud points, the better it can reflect the true potential surface signals. When statistically analyzing the potential surface laser altimetry cloud data along the terrain slope direction, it has the highest data density. Therefore, the maximum value of the laser cloud density obtained by extracting the terrain slope along the track using the current data segment and several data segments before and after it can reflect the true slope of this area.

[0025] Furthermore, the vertical segmentation window is rotated so that its extension direction is consistent with the slope direction of the corresponding data segment along the track terrain.

[0026] Due to differences in factors such as surface reflectivity, atmospheric environment, altitude, and solar incidence angle, noise density varies significantly along the track direction. Therefore, the classification density threshold cannot be set as a global threshold; instead, a local classification density threshold needs to be set for each data segment.

[0027] Further, in step 4), the maximum inter-class variance of all laser altimeter cloud densities within the corresponding data segment is used as the classification density threshold. To balance the noise removal quality of strong and weak laser altimeter data, the classification density threshold should also be limited based on the maximum noise density and the maximum signal point density within the corresponding data segment. When the classification density threshold is greater than the maximum signal density, the surface signal density is taken as the classification density threshold for the current data segment. When the classification density threshold is less than the maximum noise density, the average of the maximum noise density and the classification density threshold is taken as the final classification density threshold. The maximum signal point density is the maximum value of the laser altimeter cloud density within a certain range above and below the track-side terrain elevation of the corresponding data segment. The maximum noise density is calculated as follows: using the track-side terrain elevation of the corresponding data segment as a benchmark, a range above a certain height of this elevation is selected as the upper statistical window, and a range below a certain height of this elevation is selected as the lower statistical window. The laser altimeter cloud densities in the upper and lower statistical windows are calculated respectively, and the maximum value of all density values ​​in the upper and lower statistical windows is taken as the maximum noise density of the current data segment.

[0028] Furthermore, in step 4), the density-based clustering algorithm performs an elliptical neighborhood search along the optimal search direction within the corresponding data segment.

[0029] The elliptical neighborhood search method along the optimal search direction can better reflect the true density of potential signal points and help to accurately identify noise points.

[0030] Furthermore, the size of the major semi-axis of the elliptical neighborhood is the diameter of the laser spot used for spaceborne laser altimetry.

[0031] Using the diameter of the laser spot as the semi-major axis of the search neighborhood can better reflect the details of terrain changes, adapt to situations where the laser spot lengthens and the signal density decreases when the terrain slope is large, and also avoids introducing too much noise when calculating local density.

[0032] Furthermore, in step 3), the slope of the terrain along the track in each data segment is calculated through the following steps:

[0033] Using the linear equation y=kx+b and the plane equation Z=Ax+By+C, calculate the equation of the intersection between the vertical plane containing the centerline and the terrain plane in the corresponding data segment. Based on the intersection equation, further calculate the terrain slope along the track for that data segment.

[0034] The slope of the tangent at the midpoint of the data segment corresponding to the intersection equation is obtained by geometric calculation and used as the slope of the terrain along the track for that data segment. The scheme is simple, effective and easy to implement. Attached Figure Description

[0035] Figure 1This is a schematic diagram of the terrain slope along the track, which is a segmented fitting of densely matched point cloud from optical stereo images.

[0036] Figure 2 This is a schematic diagram of laser altimetry data and densely matched point cloud data used for segmented fitting of the terrain slope along the track;

[0037] Figure 3(a) is a schematic diagram of the terrain slope along the track for the seven data segments in the embodiment;

[0038] Figure 3(b) is a schematic diagram of the vertical partitioning window of the third data segment counting forward from the current data segment in the embodiment;

[0039] Figure 3(c) is a schematic diagram of the vertical partitioning window of the second data segment counting forward from the current data segment in the embodiment;

[0040] Figure 3(d) is a schematic diagram of the vertical partitioning window of the first data segment counting backwards from the current data segment in the embodiment;

[0041] Figure 3(e) is a schematic diagram of the vertical partitioning window of the current data segment in the embodiment;

[0042] Figure 3(f) is a schematic diagram of the vertical partitioning window of the first data segment counting backwards from the current data segment in the embodiment;

[0043] Figure 3(g) is a schematic diagram of the vertical partitioning window of the second data segment counting backwards from the current data segment in the embodiment;

[0044] Figure 3(h) is a schematic diagram of the vertical partitioning window of the third data segment counting backwards from the current data segment in the embodiment;

[0045] Figure 3(i) is a schematic diagram of the vertical partition window with the maximum point cloud density of laser altimetry in the embodiment;

[0046] Figure 4 This is a flowchart of the method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to the present invention. Detailed Implementation

[0047] The present invention will now be described in further detail with reference to the accompanying drawings.

[0048] This invention proposes a method for denoising spaceborne laser altimetry data assisted by optical stereo imagery. By utilizing densely matched point clouds generated from optical stereo images as auxiliary information, noise removal from the laser altimetry data is achieved.

[0049] 1. The method for noise removal in spaceborne laser altimetry data according to this invention is as follows:

[0050] Noise removal from laser altimetry data is aided by optical stereo imagery. First, densely matched point clouds are generated using optical stereo imagery, and residual noise within these point clouds is removed. Then, based on the orbital distance (orbit of the onboard laser satellite), latitude, longitude, and elevation information of the laser altimetry data, slope and elevation information are extracted from the dense point cloud data in segments. Finally, the extracted slope and elevation information from each data segment is used to improve the density-based clustering algorithm, effectively removing noisy points from the laser altimetry data.

[0051] like Figure 4 As shown, this is specifically achieved through steps 1.1-1.5:

[0052] Step 1.1 Generate densely matched point clouds from optical stereo images:

[0053] Using the original orientation parameters of the optical stereo image, dense matching is performed to generate a densely matched point cloud, and noise in the densely matched point cloud is filtered out.

[0054] Step 1.2 Dense matching of point cloud segmentation for fitting the terrain slope along the track:

[0055] First, the laser altimetry data is divided into segments according to a certain distance along the track. Based on the distance along the track, latitude, longitude, and elevation information of each segment, a centerline is fitted to each segment. Then, a buffer zone of a certain width is established based on the fitted centerline. The 3D coordinate data of the densely matched point cloud within the buffer zone is extracted, and the densely matched point cloud data within this buffer zone is used to fit a terrain plane, such as... Figure 1 As shown in the figure, the point cloud distribution of laser altimetry data is shown. One of the areas with a high density of point cloud is the potential surface signal point area 20, which describes the surface profile. The point cloud distributed more evenly above and below it is the noise above and below the ground. The continuous curve in the figure is the track slope 10 fitted according to the densely matched point cloud data. The gray levels of the point cloud and the terrain slope in the figure, which vary according to a certain track distance, represent different data segments. Finally, the equation (3) of the intersection line (i.e., the track slope 10) of the vertical plane where the center line is located and the spatial plane fitted by the densely matched point cloud is calculated using formulas (1) and (2), and the track slope of the corresponding data segment is calculated according to formula (4).

[0056] y = kx + b(1)

[0057] Z = Ax + By + c(2)

[0058]

[0059]

[0060] Specifically, using the latitude, longitude, and elevation coordinates of the point cloud data within each data segment, a straight line is fitted using the least squares method, which serves as the centerline of the data segment, expressed by formula (1). Using this centerline as the center, a buffer zone of a certain width is taken on both sides, and densely matched point clouds within the buffer zone are extracted, such as... Figure 2 As shown, the densely distributed gray dots in the strip-shaped area of ​​part A are laser altimetry data, and the black dots a distributed at intervals along the strip direction are points on the center line of the laser altimetry data fitting. The small gray dots distributed throughout the figure are densely matched point cloud data, and the black dots distributed in the strip-shaped area of ​​part A are densely matched point clouds in the buffer zone. The plane is fitted according to the densely matched point cloud in the buffer zone, and expressed by formula (2). Formula (1) is a vertical plane in the three-dimensional environment. By combining formula (1) and formula (2), the surface slope line can be calculated and expressed by formula (3). Formula (4) can be used to further calculate the slope of the terrain along the track.

[0061] In this embodiment, each data segment obtained by segmentation is 40 meters long, and the segmentation step size is 10 meters, so that there is a certain overlap between the data segments. Figure 1 In the buffer zone, points with different gray values ​​are mixed to some extent, and the width of the buffer zone is 40 meters.

[0062] Step 1.3 Calculate the point cloud density of the laser altimetry data:

[0063] Because the noise level of laser altimetry data is affected by various factors such as terrain slope, solar background noise, land cover, and laser beam energy, the spatial distribution of noise data varies greatly. To more accurately remove noise and extract effective laser altimetry data, the slope (the terrain slope along the track for that data segment) and terrain elevation obtained from fitting the original optical stereo image in step 1.2 can be used as auxiliary information to statistically analyze the density distribution of laser altimetry data, which will then serve as input information for the next clustering algorithm. Furthermore, due to the limited accuracy of the orientation parameters in the original optical stereo image, there is a certain deviation between the actual positioning coordinates and the true coordinates. When calculating the point cloud density based on slope and elevation information, the differences in slope and elevation between the current data segment and several preceding and following data segments must be considered.

[0064] Theoretically, potential surface laser altimetry point cloud data has the highest data density when statistically analyzed along the terrain slope. Therefore, using the slope extension direction as the neighborhood statistical direction is crucial for accurately calculating the maximum point cloud data density within a given data segment. However, due to spatial differences between the terrain slope and elevation obtained from fitting the original optical stereo image in step 1.2 and the actual terrain (e.g., ...), ... Figure 1As shown, the track-side terrain slope 10 obtained by fitting the original optical stereo image is offset from the potential surface signal point area 20 with the maximum density of the laser point cloud in the figure. Therefore, when calculating the maximum density of the current data segment, the track-side slope of the laser altimetry data needs to be optimized.

[0065] Using the length and step size of the laser altimeter data along the track designed in step 1.2, compare and analyze the fitted slope values ​​at the current data segment and several data segments before and after. If the slope value changes significantly, it indicates that the terrain cover or terrain undulation has changed significantly; otherwise, the change is small.

[0066] First, for the current data segment in the vertical direction (i.e. Figure 1 A vertical partitioning window is established along the mid-elevation direction. The length of the window is the length of the current data segment, and the height is set to 2-5 meters. The window height should not be set too small, otherwise the computational workload will be large and the statistics will be inaccurate; it should also not be set too large, otherwise it may introduce more noise, resulting in an underestimation of the density statistics. In addition, to ensure the accuracy of the estimation as much as possible, the windows should have some overlap, that is, the step size of the vertical partitioning window is smaller than the window height. Furthermore, each quadrilateral vertical partitioning window is rotated so that the top and bottom sides of the window in the elevation direction are parallel to the slope direction of the current data segment along the track, and the extension of the vertical partitioning window of the data segment along the track direction is consistent with the extension direction of the slope of the current data segment, so that the laser altimeter cloud data is statistically analyzed along the terrain slope as much as possible.

[0067] Specifically, such as Figure 3(a) , 3(b) As shown in Figures 3(c), 3(d), 3(e), 3(f), 3(g), 3(h), and 3(i), the gray dots represent laser altimeter cloud points, with the horizontal axis representing the distance along the track and the vertical axis representing the elevation. Based on the track slope and elevation of the current data segment, windows are divided vertically. The rotation angle of the window is selected from the track slope of the current data segment, establishing a vertical partition window 30 for the current data segment as shown in Figure 3(e). Furthermore, the track slopes of several data segments before and after the current data segment are used as the rotation angles of the windows, and vertical partition windows are established in the same manner. In this embodiment, as shown... Figure 3(b) , 3(c) As shown in Figures 3(d), 3(f), 3(g), and 3(h), the terrain slope along the track for a total of 7 data segments (3 segments forward and 3 segments backward) is shown in Figure 3(a). The established vertical partitioning window is shown in Figure 3(a). Figure 3(b) , 3(c)Figures 3(d), 3(f), 3(g), and 3(h) are shown. The maximum laser altimeter cloud density of all vertical partition windows (30°, 40°, 50°, 60°, 70°, 80°, and 90°) in these seven rotation directions is statistically estimated. The window with the maximum laser altimeter cloud density is shown in Figure 3(i). The terrain slope along the track corresponding to the maximum density (calculated from the spatial plane of the densely matched point cloud fitting in step 1.2) is taken as the best fitting slope for the current data segment. The elevation position of the vertical partition window with the maximum density at this time is recorded as the best terrain elevation of the laser altimeter cloud in this data segment.

[0068] Step 1.4 Coarse registration of laser altimetry data with densely matched point clouds:

[0069] Based on the statistical data from step 1.3, the maximum density corresponds to the optimal fitted slope. However, since the accuracy of densely matched point clouds is not necessarily high, especially in terms of planar accuracy, even if the optimal fitted slope is selected, it may not accurately reflect the actual terrain undulations. This problem exists for both flat and undulating terrain areas. Nevertheless, the slope and elevation along the track fitted by the densely matched point cloud can basically reflect the actual terrain conditions. That is, the optimal fitted slope and elevation are mostly consistent with the slope and elevation estimated by the densely matched point cloud (e.g., Figure 1 As shown, the terrain slope 10 along the track obtained by fitting the original optical stereo image has a certain offset from the potential surface signal point area 20 with the highest density of laser point cloud in the figure. However, the terrain undulation reflected by the terrain slope 10 along the track is basically consistent with the actual terrain undulation. Based on this, coarse registration is performed between the laser altimetry data and the dense matching point cloud. After coarse registration, a density-based clustering algorithm is used to remove noise, which can improve the denoising effect.

[0070] Because there may be some slope anomalies and elevation anomalies during the statistical process, and in rare cases, the area with the highest density may not be a surface laser altimeter point, values ​​with large elevation differences are first filtered out as gross errors when calculating the registration parameters.

[0071] Coarse registration can be achieved according to steps 1.4.1-1.4.5:

[0072] 1.4.1 Extract the optimal terrain elevation at the midpoint of each data segment from the laser altimetry data.

[0073] The method for extracting the optimal terrain elevation of a data segment has been described in step 1.3, which is to take the elevation position of the vertical subdivision window with the highest density as the optimal terrain elevation of the laser altimeter cloud in that data segment.

[0074] 1.4.2 Extract the terrain elevation along the track from the densely matched point cloud data for each data segment.

[0075] Based on the data segment length and step size set in step 1.2, and the calculated terrain profile of each data segment, calculate the elevation value at the center point of each data segment, which is used as the terrain elevation along the track for the data segment.

[0076] 1.4.3 Elevation anomalies are removed as gross errors.

[0077] The elevation values ​​at the center of the data segments calculated in steps 1.4.1 and 1.4.2 may contain elevation anomalies. Calculate the standard deviation of the elevation difference at the center point of all data segments and remove data that exceeds 3 times the standard deviation.

[0078] 1.4.4 Perform curve fitting on the terrain elevation along the track corresponding to the densely matched data.

[0079] Based on the terrain elevation along the track at the center of each data segment obtained from the densely matched data in step 1.4.2, fit the terrain curve.

[0080] 1.4.5 Translate the fitted curve according to the set step size and calculate the coarse registration parameters.

[0081] The curve fitted in step 1.4.4 is translated by a certain step size along the track. The track-side terrain elevation at the center point of the data segment is interpolated from the curve after each translation. The standard deviation of the elevation difference between the optimal terrain elevation from the laser altimetry data at the center of each data segment and the track-side terrain elevation fitted by the densely matched point cloud is calculated. The track-side translation distance with the smallest standard deviation is taken as the horizontal coarse registration parameter between the laser altimetry data and the densely matched point cloud data. The elevation error with the smallest standard deviation is taken as the elevation direction coarse registration parameter between the laser altimetry data and the densely matched point cloud data. At this point, the horizontal coarse registration parameter is the horizontal translation parameter along the track. Further, based on the track azimuth, the horizontal translation parameter needs to be decomposed into translation parameters in the x and y directions of the plane projection coordinate system, and the plane and elevation values ​​of the densely matched point cloud are translated. Assume the horizontal translation parameter is Δs, and the track azimuth angle is... The orbital slope is taken as the average value k of the slope of the fitted centerline of each data segment. The orbital azimuth angle is calculated by formula (5), and the translation parameters Δx in the x-direction and Δy in the y-direction of the plane projection coordinate system are calculated by formulas (6) and (7). By adding the translation parameters Δx and Δy to the X and Y coordinates of the densely matched point cloud, coarse registration between the densely matched point cloud and the laser altimetry data can be achieved.

[0082]

[0083]

[0084]

[0085] Step 1.5 Remove noise using a density-based clustering algorithm:

[0086] Using the densely matched point cloud after coarse registration, the track-side terrain slope is fitted again according to step 1.2. The new track-side terrain slope and potential signal points in the laser altimetry data are more closely matched. Parameters such as elevation difference, slope, point cloud density statistically obtained from each vertical window, and noise density can be comprehensively utilized to design appropriate neighborhood search parameters and classification thresholds for noise point and signal point densities for noise removal.

[0087] Specifically, this is achieved through steps 1.5.1-1.5.4:

[0088] 1.5.1 Design the clustering algorithm search region parameters.

[0089] In density-based clustering algorithms, the shape and size of the search neighborhood are crucial to algorithm performance. Elliptical neighborhood search along the terrain slope typically reflects the true density of potential signal points better. The optimal fitting slope determined in step 1.3 can be used as the search direction for the elliptical neighborhood. To better determine the optimal search direction for each point, the optimal fitting slope of the data segment corresponding to the current point and several data segments before and after it are used as candidate neighborhood search directions. A horizontal direction is also added (suitable for features with flat tops, such as buildings), using the optimal fitting slope corresponding to the data segment with the highest point density as the neighborhood search direction for the current point.

[0090] After determining the optimal neighborhood search direction for each point, the selection of the length and semi-axis parameters of the elliptical search neighborhood is highly selective, as long as they are not too small (severing the connection between terrain points) or too large (including too many noise points). Since the detector can receive reflected photons from anywhere within the laser footprint, but the laser spot lengthens and signal point density decreases when the terrain slope is steep, the length of the major semi-axis needs to be appropriately increased to capture more surface signal points. This embodiment recommends using the laser spot diameter, i.e., 17.5 meters, as the major semi-axis of the search neighborhood. This not only reflects the details of terrain changes well but also adapts to the reduced surface signal point density when the terrain slope is steep. For the minor semi-axis, if the length is too large, excessive noise will be introduced when calculating the regional noise density, leading to inaccurate local density values. This invention recommends a minor semi-axis length of 3 meters, which takes into account the accuracy of the slope fitting of the dense matching point cloud and different surface cover types, while also avoiding the introduction of excessive noise when calculating local density.

[0091] 1.5.2 Estimate the classification threshold for noise point and signal point density.

[0092] Based on the search direction and size of the elliptical search neighborhood designed in step 1.5.1, the density values ​​of all points in each data segment are calculated. The maximum inter-class variance of the laser altimeter cloud density in each data segment is used as the classification density threshold for noise points and surface signal points in the corresponding data segment. In order to take into account the noise removal quality of laser altimeter data with strong and weak beams, the classification density threshold should also be limited according to the maximum noise density and the maximum signal point density in the corresponding data segment. When the calculated classification density threshold is greater than the maximum signal point density of the current data segment, the maximum signal point density is taken as the classification density threshold of the current data segment. When the classification density threshold is less than the maximum noise density, the average value of the maximum noise density and the classification density threshold is taken as the classification density threshold of the current data segment.

[0093] The statistical method for the maximum signal point density of each data segment is as follows: taking the track-side terrain elevation determined in steps 1.2 and 1.4.2 as the benchmark, the maximum value of the laser altimeter cloud density within a certain range above and below the track-side terrain elevation of the data segment is used as the maximum signal point density of the corresponding data segment.

[0094] Data outside a certain range above and below the track elevation are generally noisy point cloud data. The statistical method for the maximum noise density of each data segment is as follows: taking the track elevation of each data segment determined in steps 1.2 and 1.4.2 as the benchmark, select the range above a certain height of the elevation as the upper statistical window, select the range below a certain height of the elevation as the lower statistical window, and take the maximum density value of all points in the upper and lower statistical windows as the maximum noise density of the current data segment.

[0095] In this embodiment, the vertical height range is 50 meters.

[0096] 1.5.3 Noise removal using density-based clustering algorithms.

[0097] The algorithm iterates through each point in the laser altimeter cloud data, calculates its density value based on the best-fit slope, removes points with a density less than a threshold as noise points, and retains points with a density greater than the threshold as surface signal points.

[0098] 1.5.4 Residual noise removal.

[0099] Because the distribution of noise points in laser altimetry data is random, in some cases a few noise points may be distributed in clusters, which may also be greater than the density threshold. In this case, it is necessary to further use the 3σ criterion and the absolute elevation value of the points to remove residual noise points in the surface signal points.

[0100] Specifically, the extracted surface signal points are first divided into segments of a certain length. In each data segment, points that are higher or lower than the terrain elevation calculated in steps 1.2 and 1.4.2 by a certain height along the track, such as 50 meters, are removed. Then, the mean and standard deviation of the elevation of the remaining surface signal points in each data segment are calculated. Points whose difference between the signal point elevation and the mean elevation exceeds 3 times the standard deviation are marked as noise points and removed from the results.

[0101] In this embodiment, each data segment obtained by segmentation is 40 meters long.

[0102] Laser altimetry technology boasts high precision, geographical independence, and all-weather operation, enabling rapid acquisition of high-precision three-dimensional coordinates of the Earth's surface globally. However, it also suffers from drawbacks such as susceptibility to solar background noise, the ability to acquire only topographic profile data, lack of texture features, and difficulty in extracting surface signal points. This invention proposes a method for denoising spaceborne laser altimetry data using optical stereo imagery. This method fully utilizes the slope and elevation information from optical stereo imagery to effectively eliminate noise in laser altimetry data, achieving a superior denoising effect.

Claims

1. A method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information, characterized in that, Includes the following steps: 1) Divide the laser altimetry data into several segments according to a certain distance along the track, and fit the centerline of the corresponding data segment based on the laser altimetry data in each data segment; 2) Establish a buffer zone of a certain width with the center line as the center, extract the three-dimensional coordinate data of the densely matched point cloud in the buffer zone, and use the densely matched point cloud in this buffer zone to fit the terrain plane. The densely matched point cloud is generated by densely matching the original orientation parameters of the optical stereo image of the same area. 3) Calculate the terrain slope along the track of the terrain profile line where each data segment is located; the terrain profile line is the intersection line of the vertical plane where the center line of the corresponding data segment is located and the terrain plane obtained by fitting the densely matched point cloud; 4) Using a density-based clustering algorithm, the laser altimeter cloud data within each data segment is iterated and the density value of each point based on the optimal search direction of the corresponding data segment is calculated. Points with density values ​​less than the classification density threshold within the corresponding data segment are designated as noise points. The optimal search direction is determined by using the slope of the terrain along the track in the data segment corresponding to the point to be determined, as well as the slope of the terrain along the track in several data segments before and after it, as the neighborhood search direction. A horizontal search direction is also added. The maximum density value calculated from each search direction is taken as the density value of that point, and the search direction corresponding to the maximum density value is taken as the optimal search direction for that point. The classification density threshold is calculated based on the density values ​​of all points within the corresponding data segment.

2. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 1, characterized in that, Step 3) further calculates the elevation along the track of the terrain profile line where each data segment is located. The elevation along the track is the elevation calculated at the center point of the corresponding data segment's terrain profile line. The best-fit slope and best terrain elevation for the corresponding data segment's laser altimetry data are then calculated using the terrain slope along the track. Before step 4), the following coarse registration steps are performed between the laser altimetry data and the densely matched point cloud data: Using the track-side terrain elevation and optimal terrain elevation calculated in step 3), the standard deviation of the elevation difference for each data segment is calculated, and elevation values ​​exceeding the set standard deviation are removed. The elevation difference is the difference between the optimal terrain elevation and the track-side terrain elevation for the corresponding data segment. After removing the excess elevation values, the track-side terrain elevation values ​​of the densely matched point cloud data at each data segment are fitted into an elevation curve, and interpolation is used to ensure that there is a corresponding elevation value at each data segment. The elevation curve is translated according to a set step size, and the standard deviation of the elevation difference for each data segment is calculated after each translation. The translation distance with the smallest standard deviation of elevation difference is taken as the coarse registration parameter. After the dense matching point cloud is corrected by applying the coarse registration parameters, the track-side terrain slope of each data segment is refitted for use in the density-based clustering algorithm in step 4). The coarse registration parameters are converted from track-side distance and elevation to three-dimensional coordinates before being used to correct the dense matching point cloud.

3. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 2, characterized in that, In step 3), the optimal fitting slope and optimal terrain elevation for the corresponding data segment are calculated based on the track-side terrain slope using the following steps: For each data segment, several vertical partitioning windows are established in the elevation direction with a set step size. The length of the window is equal to the length of the data segment, and the height is a preset window height. The set step size is less than the window height to ensure a certain degree of overlap between windows. The track-side terrain slope of the corresponding data segment, the data segment before the corresponding data segment, and the data segment after the corresponding data segment are used as the rotation direction of the vertical partitioning windows. The laser altimeter cloud density of each vertical partitioning window is calculated. The track-side terrain slope corresponding to the vertical partitioning window with the highest laser altimeter cloud density is taken as the best-fit slope of the corresponding data segment. The elevation of the center point of the vertical partitioning window with the highest laser altimeter cloud density is taken as the optimal terrain elevation.

4. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 3, characterized in that, The vertical segmentation windows are rotated so that their extension direction is consistent with the slope direction of the terrain along the track for the corresponding data segment, the data segment before the corresponding data segment, and the data segment after the corresponding data segment.

5. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 2, characterized in that, In step 4), the maximum inter-class variance of all laser altimeter cloud densities within the corresponding data segment is used as the classification density threshold. To balance the noise removal quality of strong and weak laser altimeter data, the classification density threshold should also be limited based on the maximum noise density and the maximum signal density within the corresponding data segment. When the classification density threshold is greater than the maximum signal density, the maximum signal density is taken as the classification density threshold for the current data segment. When the classification density threshold is less than the maximum noise density, the average of the maximum noise density and the classification density threshold is taken as the final classification density threshold. The maximum signal density is the maximum value of the laser altimeter cloud density within a certain range above and below the track-side terrain elevation of the corresponding data segment. The maximum noise density is calculated as follows: using the track-side terrain elevation of the corresponding data segment as a benchmark, a range above a certain height of this elevation is selected as the upper statistical window, and a range below a certain height of this elevation is selected as the lower statistical window. The laser altimeter cloud densities in the upper and lower statistical windows are calculated respectively, and the maximum value of all density values ​​in the upper and lower statistical windows is taken as the maximum noise density of the current data segment.

6. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 1, characterized in that, In step 4), the density-based clustering algorithm performs an elliptical neighborhood search along the optimal search direction within the corresponding data segment.

7. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 6, characterized in that, The size of the major semi-axis of the elliptical search neighborhood is the diameter of the laser spot used for spaceborne laser altimetry.

8. The method for denoising spaceborne laser altimetry data assisted by optical stereoscopic image information according to claim 1, characterized in that, In step 3), the slope of the terrain along the track in each data segment is calculated through the following steps: Using the linear equation y=kx+b and the plane equation Z=Ax+By+C, calculate the equation of the intersection between the vertical plane containing the centerline and the terrain plane in the corresponding data segment. Based on the intersection equation, further calculate the terrain slope along the track for that data segment.

Images