CN115825920B - ICESat-2 photon denoising method considering glacier morphology - Google Patents

ICESat-2 photon denoising method considering glacier morphology Download PDF

Info

Publication number
CN115825920B
CN115825920B CN202310092563.7A CN202310092563A CN115825920B CN 115825920 B CN115825920 B CN 115825920B CN 202310092563 A CN202310092563 A CN 202310092563A CN 115825920 B CN115825920 B CN 115825920B
Authority
CN
China
Prior art keywords
photon
photons
density
neighborhood
ransac
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310092563.7A
Other languages
Chinese (zh)
Other versions
CN115825920A (en
Inventor
黄荣刚
常瑞杰
江利明
张文典
汪汉胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Precision Measurement Science and Technology Innovation of CAS
Original Assignee
Institute of Precision Measurement Science and Technology Innovation of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Precision Measurement Science and Technology Innovation of CAS filed Critical Institute of Precision Measurement Science and Technology Innovation of CAS
Priority to CN202310092563.7A priority Critical patent/CN115825920B/en
Publication of CN115825920A publication Critical patent/CN115825920A/en
Application granted granted Critical
Publication of CN115825920B publication Critical patent/CN115825920B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

The invention discloses an ICESat-2 photon denoising method considering glacier morphology, which comprises the following steps of: the multiscale RANSAC determines the optimal photon ellipse neighborhood; an estimate of the weighted density comprising an estimate of the continuity density, the similarity density, and the blend density; determination of adaptive thresholds along track segments. The novel photon denoising method provided by the invention better solves the problem that noise photons attached near the ground surface are difficult to separate due to similar density as signal photons caused by factors such as different morphological dimensions, direction change and the like, and can effectively improve the extraction quality of the signal photons, thereby serving for detecting the elevation of the glacier surface.

Description

ICESat-2 photon denoising method considering glacier morphology
Technical Field
The invention belongs to the field of satellite-borne laser radar photon point cloud processing, and particularly relates to an ICESat-2 photon denoising method considering glacier morphology.
Background
The laser altimeter satellite ICESat-2 (The Ice, cloud and Land Elevation Satellite-2) is provided with an advanced topography laser altimeter system (Advanced Topography Laser Altimeter System, ATLAS), and a brand new micropulse multi-beam photon counting laser radar is used, so that a three-dimensional high-precision and high-density photon point cloud can be captured. One of the main scientific objectives of this task is to monitor the surface morphology of the frozen collar and its changes. ATL03 is a global geographic positioning photon data product, and under clear atmospheric conditions, the resolution along the orbit can reach 0.7 meter, thus providing unprecedented opportunities for measuring the surface morphology of fine glaciers, such as small fissures, glaciers, and the like. However, photon counting lidar has higher sensitivity, is susceptible to solar background noise, atmospheric scattering and instrument dark counting, resulting in ATL03 data products containing a large number of noise photons, and noise photons are randomly and widely distributed in space, and are difficult to reject from echo signal photons. Therefore, there is a need to further develop denoising methods to retrieve glacier surface photons with high resolution from ATL03 data products.
The existing denoising method basically assumes that the target surface photon density is larger than the background noise density, so as to extract the target surface photon. These methods have higher quality in areas of strong signal-to-noise and simple morphology, while robustness and quality in complex morphology (e.g., surface fissures) and low signal-to-noise treatments remain to be improved. In addition, existing denoising algorithms cannot effectively filter out noise photons near signal photons. From the raw photon data distribution characteristics, we found that: in addition to photon density, there are other differences in local spatial distribution, such as glacier surface spatial continuity, density scale difference, direction change, etc., and deep excavation of the device can help to improve photon denoising quality and adaptability.
Disclosure of Invention
The invention aims to provide an ICESat-2 photon denoising method considering glacier morphology aiming at the problem of difficult photon denoising in a complex morphology region in the prior art. The invention can remove noise photons attached near the ground surface to the maximum extent, can more completely retain signal photons of complex morphology areas, is beneficial to providing high-resolution surface features in the research of glacier change process, provides surface height for undisturbed ice surface, and provides a basis for evaluating surface elevation change and quality change in the range of the ice cover.
In order to solve the technical problems, the invention adopts the following technical scheme:
an ICESat-2 photon denoising method considering glacier morphology comprises the following steps:
step 1, calculating the K neighbor distance of each photon;
step 2, selecting a current central photon, and performing multi-scale RANSAC linear fitting by using a neighborhood based on a circle;
step 3, determining the optimal circle neighborhood scale of the current center photon and the major axis size and direction of the ellipse neighborhood;
step 4, calculating the continuity density of the current photon center point;
step 5, calculating the similarity density of the current photon center point;
step 6, calculating a mixing weighted density according to the continuity density and the similarity density;
step 7, dividing the along-track photons into disjoint along-track sections;
step 8, carrying out histogram statistics on photon elevation of each rail-along section, and extracting an elevation section of the candidate signal photons;
step 9, determining an adaptive density threshold according to the normalized mixed weighted density average value in the section;
and step 10, separating the signal photons from the noise photons according to the adaptive density threshold.
Step 2 as described above comprises the steps of:
step 2.1, setting the multi-scale set as
Figure SMS_1
Wherein->
Figure SMS_2
For each scale, each scale corresponds to the radius of a circle, and each photon is used as the current center photon to search photons in the circle neighborhood under each scale;
step 2.2, preliminary classification is carried out on the neighborhood photons of each scale of the current center photon by using the Gaussian mixture model and taking the K neighbor distance of the photons as a classification basis, so as to obtain a possible signal photon set under each scale
Figure SMS_3
；
Step 2.3, performing RANSAC linear fitting on possible signal photons of the current center photon under each scale by using a RANSAC algorithm to obtain the direction of the current center photon under each scale, a RANSAC fitting line and an interior point set
Figure SMS_4
。
Step 3 as described above comprises the steps of:
optimal circular neighborhood scale
Figure SMS_5
A circle neighborhood scale corresponding to the maximum value of the ratio of the number of the interior points obtained by RANSAC fitting under each scale in the step 2 to the number of possible signal photons;
if the current central photon is obtained in step 2Is determined to be effective in RANSAC linear fitting, and the major axis of the ellipse is
Figure SMS_6
The major axis direction of the ellipse is the direction of the current central photon under the optimal circle neighborhood scale, and the minor axis of the ellipse is a set fixed value;
if the current center photon is not in the interior point set obtained in step 2, it is determined that the RANSAC linear fitting is not effective, the elliptical long axis direction is defined as the horizontal direction, the elliptical long axis size is a set fixed value, and the elliptical short axis is a set fixed value.
Step 4 as described above comprises the steps of:
when the RANSAC linear fitting is effective, assigning a continuity weight according to the distance from photons in the elliptical neighborhood to the RANSAC fitting line;
when the RANSAC linear fit is not valid, the continuity weight of the photons is assigned to zero,
and counting the sum of the continuity weights of all photon points in the elliptic neighborhood as the continuity density of the current photon center point.
Step 5 as described above comprises the steps of:
step 5.1, assigning similarity weight to photons in the elliptical neighborhood of the current center photon:
searching for photons within an elliptical neighborhood
Figure SMS_7
The nearest photon; calculating photons in the vicinity of the ellipse and +.>
Figure SMS_8
The average value of the distances between the nearest photons is taken as photon spacing +.>
Figure SMS_9
The method comprises the steps of carrying out a first treatment on the surface of the Evaluating the difference of photon spacing between photons in the elliptical neighborhood and the current center photon, and distributing a similarity weight for photons in the elliptical neighborhood;
and 5.2, counting the sum of the similarity weights of all photon points in the elliptic neighborhood to be used as the similarity density of the current photon center point.
The mixing weight density in step 6 as described above is obtained based on the following formula:
Figure SMS_11
wherein: />
Figure SMS_13
And->
Figure SMS_17
Continuity density coefficient and similarity density coefficient, respectively,/->
Figure SMS_12
For the continuity density>
Figure SMS_15
For similarity density, ++>
Figure SMS_16
Oval major axis of oval neighborhood for current center photon, +.>
Figure SMS_18
Elliptical short axis of elliptical neighborhood for current center photon, +.>
Figure SMS_10
Is the current central photon +.>
Figure SMS_14
Is used to normalize the mixing weighted density.
The altitude interval of the candidate signal photon in the step 8 is as described above
Figure SMS_19
Wherein->
Figure SMS_20
For the minimum elevation value of the maximum frequency interval of the photon elevation histogram,/for the maximum frequency interval of the photon elevation histogram>
Figure SMS_21
The maximum elevation value of the maximum frequency interval of the photon elevation histogram,
Figure SMS_22
is an extension value.
The adaptive density threshold in step 9 as described above is based on the following formula:
Figure SMS_23
wherein: />
Figure SMS_24
Is an adaptive density threshold->
Figure SMS_25
For the normalized mixed weighted density average within the segment, +.>
Figure SMS_26
For density correction value +.>
Figure SMS_27
Effective fitting ratio for RANSAC linear fit within a segment, +.>
Figure SMS_28
An effective fit ratio threshold for RANSAC linear fit.
The continuity weight of photons within the elliptical neighborhood in step 4 as described above is obtained based on the following formula:
Figure SMS_30
wherein: />
Figure SMS_32
For continuity weight of photons within an elliptical neighborhood, < +.>
Figure SMS_34
For the first Gaussian kernel function standard deviation, +.>
Figure SMS_31
For photon->
Figure SMS_33
Distance to RANSAC fit line, +.>
Figure SMS_35
Is an effective fitting index of RANSAC linear fitting, if effective fitting is possible, then +.>
Figure SMS_36
If the fitting is invalid, then +.>
Figure SMS_29
。
The similarity weight for photons within the elliptical neighborhood in step 5 as described above is obtained based on the following formula:
Figure SMS_37
wherein (1)>
Figure SMS_38
For similarity weight of photons in an elliptical neighborhood, +.>
Figure SMS_39
For the second Gaussian kernel function standard deviation, +.>
Figure SMS_40
Is the difference in the spatial distribution of photons within the elliptical neighborhood from the current center photon.
Compared with the prior art, the invention has the following beneficial effects:
(1) High precision: according to the invention, the optimal elliptic neighborhood of the photon point is determined through the multiscale RANSAC, and the estimation of the continuity density, the similarity density and the mixed density is carried out by adopting the one-dimensional Gaussian kernel, so that the problem that the noise photon density attached near the ground surface is difficult to separate due to the fact that the noise photon density is similar to the signal photon density caused by factors such as different morphological scales and direction changes of glaciers is solved, the extraction quality of the signal photons can be effectively improved, and the detection of the glacier surface elevation is further serviced;
(2) Practicality: the invention can realize high-precision photon denoising of different glacier landforms, can acquire the surface gradient change of most scenes, and provides a basis for extracting the three-dimensional characteristics of the surface cracks.
Drawings
FIG. 1 is a schematic flow chart of the present invention.
FIG. 2 (a) is a graph showing the result of photon denoising in region one of the complex morphology of the jacobian glacier surface, where
Figure SMS_42
And->
Figure SMS_45
The original data diagram of the weak wave beam and the original data diagram of the strong wave beam are respectively +.>
Figure SMS_47
And->
Figure SMS_43
The method comprises the steps of respectively carrying out a weak beam mixed weighted density map and a strong beam mixed weighted density map, < >>
Figure SMS_44
And->
Figure SMS_46
Respectively a low-beam denoising result diagram and a high-beam denoising result diagram, < + >>
Figure SMS_48
And->
Figure SMS_41
Respectively carrying out a low-beam denoising result precision evaluation graph and a high-beam denoising result precision evaluation graph;
FIG. 2 (b) is a graph showing the result of photon denoising in region two of a complex morphology of the jacobian glacier surface, where
Figure SMS_50
And
Figure SMS_53
respectively the original numbers of the weak wave beamsData map and strong beam raw data map, +.>
Figure SMS_54
And->
Figure SMS_51
The method comprises the steps of respectively carrying out a weak beam mixed weighted density map and a strong beam mixed weighted density map, < >>
Figure SMS_52
And->
Figure SMS_55
Respectively a low-beam denoising result diagram and a high-beam denoising result diagram, < + >>
Figure SMS_56
And->
Figure SMS_49
Respectively carrying out a low-beam denoising result precision evaluation graph and a high-beam denoising result precision evaluation graph;
FIG. 2 (c) is a graph of the photon denoising result of the sloped region of the jacobian glacier surface, where
Figure SMS_59
And->
Figure SMS_61
The original data diagram of the weak wave beam and the original data diagram of the strong wave beam are respectively +.>
Figure SMS_63
And->
Figure SMS_58
The method comprises the steps of respectively carrying out a weak beam mixed weighted density map and a strong beam mixed weighted density map, < >>
Figure SMS_60
And->
Figure SMS_62
A low-beam denoising result diagram and a high-beam denoising result diagram respectively,
Figure SMS_64
and->
Figure SMS_57
Respectively carrying out a low-beam denoising result precision evaluation graph and a high-beam denoising result precision evaluation graph;
FIG. 2 (d) is a graph of the photon denoising result for a flat region of the jacobian glacier surface, where
Figure SMS_66
And->
Figure SMS_68
The original data diagram of the weak wave beam and the original data diagram of the strong wave beam are respectively +.>
Figure SMS_71
And->
Figure SMS_67
The method comprises the steps of respectively carrying out a weak beam mixed weighted density map and a strong beam mixed weighted density map, < >>
Figure SMS_69
And->
Figure SMS_70
Respectively a low-beam denoising result diagram and a high-beam denoising result diagram, < + >>
Figure SMS_72
And->
Figure SMS_65
And respectively carrying out a low-beam denoising result precision evaluation graph and a high-beam denoising result precision evaluation graph.
FIG. 3 (a) is a graph showing the comparison of the denoising result of the complex morphology region one with the denoising result of the RANSAC algorithm, the modified DBSCAN algorithm and the KNN algorithm, respectively, wherein
Figure SMS_75
And->
Figure SMS_76
Respectively a weak wave beam RANSAC denoising result diagram and a strong wave beam RANSAC denoising result diagram, < ->
Figure SMS_78
And->
Figure SMS_74
DBSCAN denoising result graph and DBSCAN denoising result graph with improved weak wave beam and improved strong wave beam respectively,>
Figure SMS_77
and->
Figure SMS_79
Respectively a weak wave beam KNN denoising result diagram and a strong wave beam KNN denoising result diagram, < + >>
Figure SMS_80
And->
Figure SMS_73
The denoising result diagram of the invention of the weak wave Shu Ben and the denoising result diagram of the invention of the strong wave beam are respectively;
FIG. 3 (b) is a graph showing the comparison of the denoising result of the complex morphology region II with the denoising results of the RANSAC algorithm, the modified DBSCAN algorithm and the KNN algorithm, respectively, wherein
Figure SMS_82
And->
Figure SMS_85
Respectively a weak wave beam RANSAC denoising result diagram and a strong wave beam RANSAC denoising result diagram, < ->
Figure SMS_86
And->
Figure SMS_83
DBSCAN denoising result graph and DBSCAN denoising result graph with improved weak wave beam and improved strong wave beam respectively,>
Figure SMS_84
and->
Figure SMS_87
Respectively a weak wave beam KNN denoising result diagram and a strong wave beam KNN denoising result diagram, < + >>
Figure SMS_88
And->
Figure SMS_81
The denoising result diagram of the invention of the weak wave Shu Ben and the denoising result diagram of the invention of the strong wave beam are respectively;
FIG. 3 (c) is a graph comparing the denoising result of the ramp region with the denoising results of the RANSAC algorithm, the modified DBSCAN algorithm and the KNN algorithm, respectively, wherein
Figure SMS_90
And->
Figure SMS_93
Respectively a weak wave beam RANSAC denoising result diagram and a strong wave beam RANSAC denoising result diagram, < ->
Figure SMS_95
And->
Figure SMS_91
DBSCAN denoising result graph and DBSCAN denoising result graph with improved weak wave beam and improved strong wave beam respectively,>
Figure SMS_92
and->
Figure SMS_94
Respectively a weak wave beam KNN denoising result diagram and a strong wave beam KNN denoising result diagram, < + >>
Figure SMS_96
And->
Figure SMS_89
The denoising result diagram of the invention of the weak wave Shu Ben and the denoising result diagram of the invention of the strong wave beam are respectively;
FIG. 3 (d) is a planeComparing the denoising result of the slow region with the denoising results of the RANSAC algorithm, the improved DBSCAN algorithm and the KNN algorithm respectively, wherein
Figure SMS_99
And->
Figure SMS_100
Respectively a weak wave beam RANSAC denoising result diagram and a strong wave beam RANSAC denoising result diagram, < ->
Figure SMS_102
And->
Figure SMS_98
DBSCAN denoising result graph and DBSCAN denoising result graph with improved weak wave beam and improved strong wave beam respectively,>
Figure SMS_101
and->
Figure SMS_103
Respectively a weak wave beam KNN denoising result diagram and a strong wave beam KNN denoising result diagram, < + >>
Figure SMS_104
And->
Figure SMS_97
A denoising result graph of the invention for the weak wave Shu Ben and a denoising result graph of the invention for the strong wave beam respectively.
FIG. 4 is a graph comparing the denoising result of the present invention with ATL03 confidence photons, standard land ice orbital height product ATL 06. Wherein, (a) and (c) are respectively a complex morphology region-weak beam ATL03 confidence photon graph and a strong beam ATL03 confidence photon graph, (b) and (d) are respectively a complex morphology region-weak beam inventive denoising result and ATL06 contrast graph and a strong beam inventive denoising result and ATL06 contrast graph, (e) and (g) are respectively a complex morphology region-weak beam ATL03 confidence photon graph and a strong beam ATL03 confidence photon graph, and (f) and (h) are respectively a complex morphology region-weak wave Shu Ben inventive denoising result and ATL06 contrast graph and a strong beam inventive denoising result and ATL06 contrast graph.
Detailed Description
The present invention will be further described in detail below in conjunction with the following examples, for the purpose of facilitating understanding and practicing the present invention by those of ordinary skill in the art, it being understood that the examples described herein are for the purpose of illustration and explanation only and are not intended to limit the invention.
Examples
An ICESat-2 photon denoising method considering glacier morphology specifically comprises the following steps:
step one: the K-nearest neighbor distance of the photon is calculated. For ICESat-2 ATL03 photon data, a photon spatial index was constructed using KdTrae, and the K nearest neighbor distance (KNN) for each photon was calculated.
Figure SMS_105
Formula (1)
Wherein:
Figure SMS_108
representing +.>
Figure SMS_112
The nearest +.>
Figure SMS_113
Total distance of individual photons, +.>
Figure SMS_109
And
Figure SMS_111
respectively represent +.>
Figure SMS_116
Individual photons and->
Figure SMS_117
Along-track distance of individual photons, < >>
Figure SMS_106
And->
Figure SMS_110
Respectively represent +.>
Figure SMS_114
Individual photons and->
Figure SMS_115
Photon elevation of individual photons. />
Figure SMS_107
Is to calculate->
Figure SMS_118
The important parameter of (a) is to select the number of neighboring photons, which is essentially a smoothing factor that acts to smooth the spatial heterogeneity inside the same class of photons. If->
Figure SMS_119
If the distance between the signal photon and the noise photon is too small, the smoothing effect is not obvious, the distance sum of the signal photon and the noise photon is relatively similar, and the signal photon and the noise photon are difficult to distinguish; if->
Figure SMS_120
Too large, the difference between the noise photons and the signal photons is excessively smoothed, so that the photon denoising is incomplete.
Step two: and (5) performing multiscale RANSAC fitting. After the K-nearest neighbor distance of each photon is obtained, a multi-scale RANSAC linear fit is performed, and since a circle-based neighborhood is used, each photon is as follows.
First, a multi-scale set is set as
Figure SMS_121
Wherein->
Figure SMS_122
For each dimension, and each dimension corresponds to a radius of a circle. And taking each photon in the photon data as a current center photon, and searching photons in a circle adjacent region under each scale. />
Then, in order to reduce the influence of background noise on the fitting result, preliminary classification is carried out on the neighborhood photons of each scale of the current center photon by using a Gaussian Mixture Model (GMM) and taking the K neighbor distance of the photons as a classification basis so as to obtain a possible signal photon set under each scale
Figure SMS_123
。
Finally, carrying out RANSAC linear fitting on possible signal photons of the current center photon under each scale by using a RANSAC algorithm to obtain the direction of the current center photon under each scale, a RANSAC fitting line and an interior point set
Figure SMS_124
。
Step three: the optimal circular neighborhood scale of the current center photon and the major axis size and direction of the elliptical neighborhood are determined. Optimal circular neighborhood scale
Figure SMS_125
And (3) for the circle neighborhood scale corresponding to the maximum value of the ratio of the number of internal points obtained by RANSAC fitting to the number of possible signal photons in each scale in the second step, the larger the ratio is, the larger the probability that the internal point obtained by RANSAC fitting is a signal photon in the scale is, and the better the fitting result is. Considering that signal photons are densely distributed in the horizontal direction, in order to restrain influence of noise photons in the vertical direction on subsequent signal photon density estimation, a circular neighborhood of a current center photon is transformed into an elliptical neighborhood.
If the current central photon is in the inner point set obtained in the second step, the RANSAC linear fitting is judged to be effective, and the major axis of the ellipse is
Figure SMS_126
The elliptical long axis direction is the direction of the current center photon under the optimal circular neighborhood scale.
If the current center photon is not in the interior point set obtained in the second step, it is determined that the RANSAC linear fitting is not valid, in which case the elliptical long axis direction is defined as the horizontal direction and the elliptical long axis size is set to a set fixed value (e.g., 20 m). Further, the minor axis of the ellipse is set to a fixed value regardless of whether the center photon is the set of interior points obtained in step two.
Figure SMS_127
Formula (2)
Wherein:
Figure SMS_130
is->
Figure SMS_132
Inner points fitting on scale +.>
Figure SMS_133
Quantity of->
Figure SMS_129
Is->
Figure SMS_131
Possible signal photons +.derived from the Gaussian mixture model on the scale>
Figure SMS_134
Quantity of->
Figure SMS_136
For the scale number>
Figure SMS_128
For the scale with the highest fitting ratio, i.e. +.>
Figure SMS_135
The largest corresponding scale, argmax, is a function of the function parameters.
Step four: a continuity density estimate of the current photon center point. The elliptic neighborhood of the current central photon is determined in a self-adaptive manner according to the step three, when the RANSAC linear fitting is effective, the possible signal photons of the local surface in the elliptic neighborhood can be considered to be smooth, continuous and linear, and the continuity weight is distributed according to the distance from the photons in the elliptic neighborhood to the RANSAC fitting line by utilizing a one-dimensional Gaussian kernel function, and the smaller the distance is, the larger the weight is. Furthermore, when the RANSAC linear fit is not effective, the local surface is not smooth or present, and the continuity weight of the photons is directly assigned to zero as in equation (3).
Figure SMS_137
Formula (3)
Wherein:
Figure SMS_138
continuity weight indicating to which photons within the elliptical neighborhood are assigned, +.>
Figure SMS_139
For the first Gaussian kernel function standard deviation, +.>
Figure SMS_140
For photon->
Figure SMS_141
Distance to RANSAC fit line, +.>
Figure SMS_142
Is an effective fitting index of RANSAC linear fitting, if effective fitting is possible, then +.>
Figure SMS_143
If the fitting is invalid, then +.>
Figure SMS_144
。/>
After assigning a continuity weight to all photons in the elliptical neighborhood, the continuity-based weighted density may be calculated according to equation (4), i.e., for each current center photon, the sum of the continuity weights of all photon points in the elliptical neighborhood is counted.
Figure SMS_145
Formula (4)
Wherein: photons (photon)
Figure SMS_146
From the current central photon->
Figure SMS_147
Photon set in the vicinity of ellipse->
Figure SMS_148
，/>
Figure SMS_149
Is photon point in oval neighborhood->
Figure SMS_150
Is/are continuity weights of->
Figure SMS_151
For the current photon center point->
Figure SMS_152
Is a continuous density of (c).
Step five: and estimating the similarity density of the current photon center point. Within the elliptical neighborhood of the current center photon, the properties of the photons may be non-uniform, being a mixture of signal photons and noise photons. In this way, the density estimation of signal photons and nearby noise photons is severely affected and their density differences become smaller. Thus, assigning different weights to different classes of photons in the density estimation will improve the density difference between the signal and nearby noise photons. In general, the spatial distribution of photons having the same class of properties is approximately uniform in a localized region. That is, the photon spacing of one signal photon is about the same as the photon spacing of an adjacent signal photon, and the same rule should be followed by adjacent noise. Therefore, the difference between the photon spacing in the elliptical neighborhood of the current center photon and the photon spacing of the current center photon is used to evaluate the weight using the photon spacing as a similarity index. For a current central photon
Figure SMS_153
Assigning photons within its elliptical neighborhoodThe steps of similarity weighting are as follows. 1) Is a photon in the vicinity of an ellipse->
Figure SMS_154
Search->
Figure SMS_155
The nearest photon; 2) Calculating photon +.>
Figure SMS_156
And to this
Figure SMS_157
The average value of the distances between the nearest photons is taken as photon spacing +.>
Figure SMS_158
The method comprises the steps of carrying out a first treatment on the surface of the 3) Evaluating photons in the vicinity of the ellipse and the current center photon +.>
Figure SMS_159
And the difference of photon spacing between the two photons is used for distributing a similarity weight for photons in the oval neighborhood through a one-dimensional Gaussian kernel function, as shown in a formula (5).
Figure SMS_160
Formula (5)
Wherein:
Figure SMS_161
similarity weight indicating that photons within the neighborhood of ellipse are assigned +.>
Figure SMS_162
A second Gaussian kernel function standard deviation, wherein the second Gaussian kernel function standard deviation is only numerically different from the first Gaussian kernel function standard deviation, +.>
Figure SMS_163
Is a photon in the vicinity of an ellipse->
Figure SMS_164
And at present inHeart photon->
Figure SMS_165
I.e. the difference in photon spacing.
After assigning similarity weights to all photon points in the ellipse neighborhood, the current center photon can be calculated according to equation (6)
Figure SMS_166
A weighted density based on similarity.
Figure SMS_167
Formula (6)
Wherein: photons (photon)
Figure SMS_168
From the current central photon->
Figure SMS_169
Photon set in the vicinity of ellipse->
Figure SMS_170
，/>
Figure SMS_171
Is photon point in oval neighborhood->
Figure SMS_172
Similarity weight of->
Figure SMS_173
For the current photon center point->
Figure SMS_174
Is a similarity density of (c).
Step six: and (5) mixing weighted density estimation. The continuity density reflects the shape characteristics of the local surface, while the similarity density expresses the spatial distribution difference between photon points. Thus, the two types of weighted densities are combined as a mixed weighted density. The mixed weighting density is advantageous for amplifying the difference between signal photons and noise photons. Furthermore, since the major axis size of the elliptical neighborhood may be different for photons in different regions, the hybrid weighting density should be normalized.
Figure SMS_175
Formula (7)
Wherein:
Figure SMS_178
and->
Figure SMS_179
Continuity density coefficient and similarity density coefficient, respectively,/->
Figure SMS_181
For the current central photon->
Figure SMS_177
Oval major axis of oval neighborhood, +.>
Figure SMS_180
For the current central photon->
Figure SMS_182
Elliptic minor axis of elliptic neighborhood, +.>
Figure SMS_183
Is the current central photon +.>
Figure SMS_176
Is used to normalize the mixing weighted density.
Step seven: segmentation of the along-track photons. To overcome the effect of the region density difference between signal photons and noise photons on threshold selection, the off-track photons are divided into disjoint off-track segments. In this way, the signal photon density, the noise photon density and the difference between the two densities will be kept well consistent within one along-track section. In addition, the local surface shape will be relatively simple and the elevation of the signal photons will be more concentrated.
Step eight: extraction of candidate signal photons. Rail-following lightAnd after the sub-segments are segmented, carrying out histogram statistics on the photon elevation of each along-track segment so as to extract the altitude interval of the candidate signal photons. For an along-track segment, a histogram based on photon elevation is first generated, where the abscissa is along-track distance and the ordinate is photon elevation. Then, a photon elevation interval is selected
Figure SMS_184
As candidate signal photons, wherein +.>
Figure SMS_185
For the minimum elevation value of the maximum frequency interval of the photon elevation histogram,/for the maximum frequency interval of the photon elevation histogram>
Figure SMS_186
The maximum elevation value of the maximum frequency interval of the photon elevation histogram is the interval with the maximum photon quantity. In order to ensure that all signal photons are contained in the altitude interval, an expansion value is set>
Figure SMS_187
Expanding the altitude interval, wherein the altitude interval of the final candidate signal photon is
Figure SMS_188
。
Step nine: determination of an adaptive density threshold. In general, there is a density difference between the signal photons and the noise photons in each of the along-track segments, the density threshold
Figure SMS_190
May be set as the average of all photon densities within the along-track segment. However, the density difference between signal photons and noise photons is different for different along-track segments. To improve->
Figure SMS_194
Is added to the density threshold value by a correction value +.>
Figure SMS_197
To accommodate along-track segments of different data quality, as shown in equation (8). In the present invention, the data is of good quality when most of the photons and photons within their elliptical neighbors are able to meet the RANSAC effective fit in step three. Therefore, the effective fitting ratio +.>
Figure SMS_191
For describing how many photons in the along-track segment can be effectively fit to photons within its elliptical neighborhood. When->
Figure SMS_192
Greater than threshold->
Figure SMS_195
In this case, it is shown that the data quality of the track-following section is better, and in order to remove the noise around the signal photons better, the density threshold can be set to a higher value, so that a correction value +_ is added>
Figure SMS_198
. Conversely, when +.>
Figure SMS_189
Less than threshold->
Figure SMS_193
Data quality is not considered good and the signal to noise ratio may be low. Therefore, in order to ensure adequate extraction of signal photons, a lower density threshold should be set, the density threshold is subtracted by a correction value +.>
Figure SMS_196
。
Figure SMS_199
Formula (8)
Wherein:
Figure SMS_200
is an adaptive density threshold->
Figure SMS_201
For the normalized mixed weighted density average within a segment,
Figure SMS_202
for density correction value +.>
Figure SMS_203
Effective fitting ratio for RANSAC linear fit within a segment, +.>
Figure SMS_204
And determining an effective fit ratio threshold of the RANSAC linear fit of the data quality condition for the effective fit ratio.
Step ten: separation of signal photons from noise photons. After the self-adaptive density threshold value of each along-track section is obtained, judging photons with photon density in each along-track section larger than the self-adaptive density threshold value as signal photons; photons having a photon density less than the adaptive density threshold are determined to be noise photons.
In order to illustrate the effectiveness of the invention, the strong and weak signal data of the complex morphology area, the slope area and the gentle area of the surface of the jacobian glacier are selected for testing. Table 1 is a parameter setting specification table.
Table 1 parameter setting explanatory table
Figure SMS_205
The invention provides a denoising result with higher precision, can process data with different qualities, and can reach more than 0.9 in F1 fraction of various topographic structures. The overall evaluation of the denoising results for the different regions is shown in table 2.
Table 2 the denoising result accuracy evaluation table of the present invention
Figure SMS_206
Compared with the prior art, the method has better rejection effect on noise photons attached near the ground surface, which provides more accurate high-resolution ground surface photons, thereby acquiring more accurate ground height information. Tables 3-5 are precision evaluation tables of denoising results based on the RANSAC algorithm, the modified DBSCAN algorithm, and the KNN algorithm, respectively.
Table 3 precision evaluation table of denoising result based on RANSAC algorithm
Figure SMS_207
/>
Table 4 precision evaluation table of denoising results based on improved DBSCAN algorithm
Figure SMS_208
Table 5 precision evaluation table of denoising result based on KNN algorithm
Figure SMS_209
The specific embodiments described herein are offered by way of example only to illustrate the spirit of the invention. Those skilled in the art may make various modifications or additions to the described embodiments or substitutions thereof without departing from the spirit of the invention or exceeding the scope of the invention as defined in the accompanying claims.

Claims (4)

1. An ICESat-2 photon denoising method considering glacier morphology, which is characterized by comprising the following steps of:
step 1, calculating the K neighbor distance of each photon;
step 2, selecting a current central photon, and performing multi-scale RANSAC linear fitting by using a neighborhood based on a circle;
step 3, determining the optimal circle neighborhood scale of the current center photon and the major axis size and direction of the ellipse neighborhood;
step 4, calculating the continuity density of the current photon center point;
step 5, calculating the similarity density of the current photon center point;
step 6, calculating a mixing weighted density according to the continuity density and the similarity density;
step 7, dividing the along-track photons into disjoint along-track sections;
step 8, carrying out histogram statistics on photon elevation of each rail-along section, and extracting an elevation section of the candidate signal photons;
step 9, determining an adaptive density threshold according to the normalized mixed weighted density average value in the section;
step 10, separating the signal photons from noise photons according to the adaptive density threshold,
the step 2 comprises the following steps:
step 2.1, setting the multi-scale set as
Figure QLYQS_1
Wherein->
Figure QLYQS_2
For each scale, each scale corresponds to the radius of a circle, and each photon is used as the current center photon to search photons in the circle neighborhood under each scale;
step 2.2, preliminary classification is carried out on the neighborhood photons of each scale of the current center photon by using the Gaussian mixture model and taking the K neighbor distance of the photons as a classification basis, so as to obtain a possible signal photon set under each scale
Figure QLYQS_3
；
Step 2.3, performing RANSAC linear fitting on possible signal photons of the current center photon under each scale by using a RANSAC algorithm to obtain the direction of the current center photon under each scale, a RANSAC fitting line and an interior point set
Figure QLYQS_4
，
The step 3 comprises the following steps:
optimal circular neighborhood scale
Figure QLYQS_5
A circle neighborhood scale corresponding to the maximum value of the ratio of the number of the interior points obtained by RANSAC fitting under each scale in the step 2 to the number of possible signal photons;
if the current center photon is in the inner point set obtained in the step 2, the RANSAC linear fitting is judged to be effective, and the ellipse major axis is as follows
Figure QLYQS_6
The major axis direction of the ellipse is the direction of the current central photon under the optimal circle neighborhood scale, and the minor axis of the ellipse is a set fixed value;
if the current central photon is not in the interior point set obtained in the step 2, the RANSAC linear fitting is judged to be invalid, the elliptical long axis direction is defined as the horizontal direction, the elliptical long axis size is a set fixed value, the elliptical short axis is a set fixed value,
the step 4 comprises the following steps:
when the RANSAC linear fitting is effective, assigning a continuity weight according to the distance from photons in the elliptical neighborhood to the RANSAC fitting line;
when the RANSAC linear fit is not valid, the continuity weight of the photons is assigned to zero,
the sum of the continuity weights of all photon points in the elliptical neighborhood is counted as the continuity density of the current photon center point,
the continuity weight of photons in the elliptical neighborhood in the step 4 is obtained based on the following formula:
Figure QLYQS_8
wherein: />
Figure QLYQS_10
For continuity weight of photons within an elliptical neighborhood, < +.>
Figure QLYQS_12
Is the firstStandard deviation of Gaussian kernel function->
Figure QLYQS_9
For photon->
Figure QLYQS_11
Distance to RANSAC fit line, +.>
Figure QLYQS_13
Is an effective fitting index of RANSAC linear fitting, if effective fitting is possible, then +.>
Figure QLYQS_14
If the fitting is invalid, then +.>
Figure QLYQS_7
，
Said step 5 comprises the steps of:
step 5.1, assigning similarity weight to photons in the elliptical neighborhood of the current center photon:
searching for photons within an elliptical neighborhood
Figure QLYQS_15
The nearest photon; calculating photons in the vicinity of the ellipse and +.>
Figure QLYQS_16
The average value of the distances between the nearest photons is taken as photon spacing +.>
Figure QLYQS_17
The method comprises the steps of carrying out a first treatment on the surface of the Evaluating the difference of photon spacing between photons in the elliptical neighborhood and the current center photon, and distributing a similarity weight for photons in the elliptical neighborhood;
step 5.2, counting the sum of the similarity weights of all photon points in the elliptic neighborhood as the similarity density of the current photon center point,
the similarity weight of photons in the ellipse neighborhood in the step 5 is obtained based on the following formula:
Figure QLYQS_18
wherein (1)>
Figure QLYQS_19
For similarity weight of photons in an elliptical neighborhood, +.>
Figure QLYQS_20
For the second Gaussian kernel function standard deviation, +.>
Figure QLYQS_21
Is the difference in the spatial distribution of photons within the elliptical neighborhood from the current center photon.
2. The method of denoising the ICESat-2 photons taking account of glacier morphology according to claim 1, wherein the mixing weight density in step 6 is obtained based on the following formula:
Figure QLYQS_24
wherein: />
Figure QLYQS_27
And->
Figure QLYQS_29
Continuity density coefficient and similarity density coefficient, respectively,/->
Figure QLYQS_23
For the continuity density>
Figure QLYQS_26
For similarity density, ++>
Figure QLYQS_28
Oval major axis of oval neighborhood for current center photon, +.>
Figure QLYQS_30
Elliptical short axis of elliptical neighborhood for current center photon, +.>
Figure QLYQS_22
Is the current central photon +.>
Figure QLYQS_25
Is used to normalize the mixing weighted density.
3. The method of denoising the ICESat-2 photons taking into account glacier morphology according to claim 2, wherein the altitude interval of the candidate signal photons in step 8 is
Figure QLYQS_31
Wherein->
Figure QLYQS_32
For the minimum elevation value of the maximum frequency interval of the photon elevation histogram,/for the maximum frequency interval of the photon elevation histogram>
Figure QLYQS_33
For the maximum elevation value of the maximum frequency interval of the photon elevation histogram,/for the maximum frequency interval of the photon elevation histogram>
Figure QLYQS_34
Is an extension value.
4. A method of denoising ICESat-2 photons, according to claim 3, wherein the adaptive density threshold in step 9 is based on the following formula:
Figure QLYQS_35
wherein: />
Figure QLYQS_36
Is an adaptive density threshold->
Figure QLYQS_37
For the normalized mixed weighted density average within the segment, +.>
Figure QLYQS_38
For density correction value +.>
Figure QLYQS_39
Effective fitting ratio for RANSAC linear fit within a segment, +.>
Figure QLYQS_40
An effective fit ratio threshold for RANSAC linear fit. />
CN202310092563.7A 2023-02-10 2023-02-10 ICESat-2 photon denoising method considering glacier morphology Active CN115825920B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310092563.7A CN115825920B (en) 2023-02-10 2023-02-10 ICESat-2 photon denoising method considering glacier morphology

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310092563.7A CN115825920B (en) 2023-02-10 2023-02-10 ICESat-2 photon denoising method considering glacier morphology

Publications (2)

Publication Number Publication Date
CN115825920A CN115825920A (en) 2023-03-21
CN115825920B true CN115825920B (en) 2023-05-05

Family

ID=85520954

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310092563.7A Active CN115825920B (en) 2023-02-10 2023-02-10 ICESat-2 photon denoising method considering glacier morphology

Country Status (1)

Country Link
CN (1) CN115825920B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116243273B (en) * 2023-05-09 2023-09-15 中国地质大学（武汉） Photon counting laser radar data filtering method for vegetation canopy extraction

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5520786B2 (en) * 2010-11-11 2014-06-11 株式会社パスコ Laser density distribution estimation apparatus, laser density distribution estimation method, and program
CN109799494B (en) * 2017-11-17 2020-09-29 中国林业科学研究院资源信息研究所 Satellite-borne photon counting laser radar data rapid denoising and filtering method
CN111665517B (en) * 2020-05-29 2022-11-18 同济大学 Density statistics-based single photon laser height finding data denoising method and device
CN112986964B (en) * 2021-02-26 2023-03-31 北京空间机电研究所 Photon counting laser point cloud self-adaptive denoising method based on noise neighborhood density
CN114355367A (en) * 2022-01-10 2022-04-15 中国人民解放军61540部队 Method for measuring shallow sea water depth based on satellite-borne single photon laser radar data
CN114779215A (en) * 2022-04-21 2022-07-22 昆明理工大学 Data denoising method for spaceborne photon counting laser radar in planting coverage area
CN115222949B (en) * 2022-09-21 2022-12-06 自然资源部第一海洋研究所 Shallow sea area photon denoising method based on laser satellite data

Also Published As

Publication number Publication date
CN115825920A (en) 2023-03-21

Similar Documents

Publication Publication Date Title
CN108957453B (en) High-precision moving target imaging and identifying method based on multi-target tracking
CN111830506B (en) K-means clustering algorithm-based sea surface wind speed method
US7787657B2 (en) SAR ATR treeline extended operating condition
CN108681525B (en) Road surface point cloud intensity enhancing method based on vehicle-mounted laser scanning data
CN101727662B (en) SAR image nonlocal mean value speckle filtering method
US7116265B2 (en) Recognition algorithm for the unknown target rejection based on shape statistics obtained from orthogonal distance function
CN111665517B (en) Density statistics-based single photon laser height finding data denoising method and device
CN103971127B (en) Forward-looking radar imaging sea-surface target key point detection and recognition method
CN104050477B (en) Infrared image vehicle detection method based on auxiliary road information and significance detection
CN115825920B (en) ICESat-2 photon denoising method considering glacier morphology
CN108171193B (en) Polarized SAR (synthetic aperture radar) ship target detection method based on super-pixel local information measurement
CN110379007B (en) Three-dimensional highway curve reconstruction method based on vehicle-mounted mobile laser scanning point cloud
CN110516646B (en) Superglacial moraine covering type glacier identification method combining polarization decomposition and topographic features
CN109100697B (en) Target condensation method based on ground monitoring radar system
CN113034569A (en) Point cloud data-based ship overrun early warning method and system
CN107678025B (en) Sea wave height calculation method and device, storage medium and processor
CN108983194B (en) Target extraction and condensation method based on ground monitoring radar system
CN113256990B (en) Method and system for collecting road vehicle information by radar based on clustering algorithm
CN113466827B (en) Denoising method based on improved local sparse algorithm
CN113269889A (en) Self-adaptive point cloud target clustering method based on elliptical domain
CN111080647B (en) SAR image segmentation method based on adaptive sliding window filtering and FCM
CN116165635A (en) Denoising method for photon cloud data of different beams under daytime condition of multistage filtering algorithm
CN115436966A (en) Batch extraction method for laser radar reference water depth control points
CN113281716B (en) Photon counting laser radar data denoising method
CN114779235A (en) Road boundary detection method and system based on vehicle-mounted millimeter wave radar

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant