Ground surface identification method based on photon point cloud in polar sea ice observation scene

Abstract

The present application relates to the technical field of remote sensing data processing, and provides a kind of ground surface identification method based on photon point cloud in polar sea ice observation scene, comprising: obtaining single-photon laser altimetry data to construct two-dimensional photon profile;First scale window is divided, threshold is determined based on density gradient field, and low-density noise is filtered to obtain intermediate signal set;Elevation window is constructed in second scale window, and high-density and low-density tailing double-end outlier subsets outside the window are removed to obtain high-confidence signal set;Resampling elevation and sliding classification features are extracted;Filtering candidate centroid and optical high-resolution image registration and gating;Label consistency check is carried out using multi-temporal image, and spatial continuity is optimized;Physical verification is completed in combination with snow and ice radiation transfer model to construct training set;Sea ice, ice channel and thin ice identification results are output by inputting classifier.The present application improves the physical reliability of training set label, reduces the influence of optical label mismatch, and improves the accuracy of polar ground surface identification.

Description

Technical Field

[0001] This invention relates to the field of remote sensing data processing technology, and in particular to a method for surface identification based on photon point clouds in polar sea ice observation scenarios. Background Technology

[0002] Accurate identification of polar sea ice and interglacial channels is crucial for estimating sea ice mass balance, studying ice-sea interactions, and parameterizing climate models. The next-generation spaceborne single-photon lidar (ICESat-2) provides massive amounts of point cloud data for reconstructing polar surface elevation profiles with centimeter-level precision through high-frequency photon counting altimetry along its orbit. However, in polar sea ice observation scenarios, single-photon lidar altimetry point cloud data is easily affected by solar background photons, detector noise, and variations in surface reflectivity, making the separation of signal photons from background photons challenging. Furthermore, the limited differences in orbital elevation, local photon density, and optical image spectral characteristics among sea ice, thin ice, and interglacial channels reduce the separability between different surface types, posing several technical bottlenecks to directly converting photon point clouds into semantically labeled surface classification results.

[0003] In photonic cloud denoising, traditional methods typically rely on single-threshold density segmentation, which uniformly labels low-density photons below a certain statistical cutoff threshold as background noise. This single-boundary processing logic tends to fail under conditions of global increases in diffuse background noise during polar day and concurrent occurrences of transient high-energy artifacts from detectors: when the threshold is set too low, diffuse low-density ground object photons composed of remotely sensed reflected light fields are incorrectly retained; when the threshold is set too high, high-density isolated clusters caused by cosmic rays or detector dark current spikes are misjudged as strong surface echoes, forming artifact groups that are difficult to separate from the real surface signal in the density domain. Furthermore, relying solely on density thresholds for high-end truncation easily leads to the incorrect removal of real surface signal photons in high-density states. This problem has not been specifically addressed in the publicly available ICESat-2 data processing chain, resulting in the purity of the high-confidence signal set being cross-contaminated by outlier photon groups at both ends.

[0004] In terms of training set construction, supervised classification schemes for distinguishing sea ice from interglacial channels typically rely on high-quality physical labels to support the training set construction. Existing literature usually uses direct geographic coordinate matching to spatially associate laser beams with high-resolution optical image pixels to obtain classification labels. However, the narrow, long, and winding linear morphology of polar interglacial channels and the asymmetrical coverage of laser beam footprints along and across the track create a spatial mismatch. Furthermore, the solar geometry conditions at the observation time and cloud shadow interference cause spectral response shifts in the optical images. Schemes that rely solely on latitude and longitude nearest neighbor matching fail to consider the distribution differences of the two modal signals in the abstract feature space and lack physical consistency verification of the registration results. This can easily lead to the erroneous transfer of optical pixel labels with transiently changing spectral attributes to laser surface echoes within the corresponding time window, thereby contaminating the classifier's training benchmark.

[0005] Furthermore, after obtaining the surface classification results, the calculation of sea ice freeboard typically relies on the statistical percentile estimation of the elevation profile along the track, assuming a fixed low percentile as the sea surface reference. This global statistical assumption lacks physical anchor support in thick ice-covered areas where inter-ice channels are sparse or completely absent, and it also eliminates the anisotropy of local sea surface undulations and the dynamic smoothing effect within the channels, introducing systematic errors that deviate from the true sea surface height. Summary of the Invention

[0006] Therefore, it is necessary to address the problems of incomplete denoising of photon point clouds, easy contamination of training set labels, and low accuracy of surface identification in existing technologies under polar sea ice observation scenarios. This paper proposes a surface identification method based on photon point clouds for polar sea ice observation scenarios. By simultaneously suppressing diffuse reflection noise and transient artifacts through double-ended outlier subset removal, and by constructing a high-confidence training set using multi-layer progressive verification, high-precision surface identification is achieved.

[0007] This invention provides a surface identification method based on photon point clouds in polar sea ice observation scenarios, comprising: Acquire geophysically corrected polar single-photon laser altimetry point cloud data, calculate the orbital distance and residual elevation of each photon, and construct a two-dimensional photon profile of orbital distance-residual elevation. The two-dimensional photon profile is divided into first-scale windows along the track direction. Within any first-scale window, the density gradient field of the local density distribution of photons is calculated. The minimum value of the boundary of the connected graph corresponding to the change from negative to positive density gradient is extracted as the global density segmentation threshold. Photons with local density lower than the global density segmentation threshold are filtered out to obtain the intermediate signal photon set. For the intermediate signal photon set, a second-scale window is defined along the track direction to perform local analysis. An elevation window with positive and negative preset elevation thresholds is constructed based on the average elevation or modal elevation of the intermediate signal photon set. Double-ended outliers corresponding to high-density isolated tips of transient high-energy mutations are removed from the second-scale window, where the density fluctuation is at the bottom of the low-density tail and the local density exceeds the preset high-density threshold outside the elevation window. This is done to simultaneously suppress diffuse reflection noise and transient artifacts, and obtain a high-confidence signal set. Using the local photon density structure within the high-confidence signal set as a robust constraint, resampled elevations are extracted along the track direction using overlapping sliding windows, and a sliding classification feature set with the overlapping sliding windows as statistical units is generated. Within the feature space formed by the sliding classification feature set, based on the constraints of the local curvature of the local density field of each sliding window and the temporal persistence along the track, the local density maxima window is selected as the candidate cluster physical centroid. The candidate cluster physical centroid is registered with the spatiotemporally matched optical high-resolution remote sensing image in the real pixel space. The registration includes: projecting the photon feature set of the sliding window and the spectral texture feature set of the optical high-resolution remote sensing image pixels into a common feature subspace, and performing gated filtering based on the statistical distance and observation geometric angle difference in the common feature subspace. For the registration results that pass through the gated filtering, the label consistency is checked using the optical high-resolution remote sensing image at multiple time phases before and after the acquisition time of the polar single-photon laser altimeter cloud data, and conflict resolution is performed by combining the spatial label continuity of adjacent centroids along the track. The multispectral reflectance of the high-resolution optical remote sensing image pixels is input into the ice and snow radiative transfer model. The physical consistency of the registration results verified above is checked by using theoretical photon return intensity and theoretical elevation broadening. The physical centroids of the clusters that pass the verification are assigned physical labels to the entities to construct a training set. The training set is used as input to the classifier to perform supervised classification on all the overlapping sliding windows, and the surface identification results are output, which include three types: sea ice, interglacial waterways, and thin ice.

[0008] In one embodiment, the step of dividing the two-dimensional photonic profile into a first-scale window includes: A nonlinear adaptive compression transformation is applied to the orbital coordinates in the two-dimensional photon profile, and this transformation is expressed as: , In the formula, The coordinates along the track before transformation. These are the transformed coordinates along the track. Based on the scaling ratio, The background photon rate is based on neighborhood statistics. The preset environmental adjustment coefficient is used; and, the basic scaling ratio is set differently for the strong beam and weak beam of the single-photon altimeter system. .

[0009] In one embodiment, the step of removing the double-ended outlier subset corresponding to the high-density isolated tip of a transient high-energy mutation, where the density fluctuation within the second scale window is at the bottom of a low-density tail and outside the elevation window with a local density exceeding a preset high-density threshold, includes: Within the second scale window, the density value corresponding to the preset low percentile is used as the low-end cutoff threshold to remove photons with local densities lower than the low-end cutoff threshold, thereby eliminating diffuse reflection noise. Outside the elevation window, the density value corresponding to the preset high percentile is used as the preset high density threshold to remove photons with local densities higher than the preset high density threshold and spatially isolated patches, thus eliminating transient high-energy mutation artifacts.

[0010] In one embodiment, the step of extracting resampled elevations along the orbital direction using an overlapping sliding window, with the photon local density structure within the high-confidence signal set as a robust constraint, includes: Within each of the overlapping sliding windows, a Tukey dual-weight robust estimator is introduced to calculate the resampled elevation according to the following robust density-weighted formula. : , In the formula, For the first overlapping sliding window The residual elevation of a high-confidence signal photon. This represents the local density value of the photon. This refers to the Tukey robust weighting factor, which is updated iteratively based on photon elevation residuals. This represents the total number of high-confidence signal photons within the overlapping sliding window.

[0011] In one embodiment, the sliding classification feature set further includes: The track slope features are extracted from the first-order forward difference of the resampled elevations of adjacent overlapping sliding windows, and the residual elevation skewness features are extracted from the third-order central moments of the residual elevations of high-confidence signal photons within the overlapping sliding windows.

[0012] In one embodiment, during the local analysis performed within the defined second-scale window, the ratio of the number of signal photons to the total number of photons within the second-scale window is calculated as the signal photon proportion. When the signal photon proportion is higher than a preset saturation threshold, the detector is determined to be saturated, and a post-pulse denoising mechanism based on dynamic optimization of residual elevation is triggered within the second-scale window: the second derivative of the photon residual elevation distribution histogram within the second-scale window is calculated, the minimum point of the valley between the primary surface reflection peak and the secondary post-pulse peak is identified, the elevation value corresponding to the minimum point of the valley is established as the adaptive truncation threshold, and photons with residual elevation lower than the adaptive truncation threshold are determined to be post-pulse photons generated by detector saturation and are removed.

[0013] In one embodiment, projecting the photon feature set of the sliding window and the spectral texture feature set of the optical high-resolution remote sensing image pixels onto a common feature subspace includes: Canonical correlation analysis is performed on the photon feature set and the spectral texture feature set respectively to extract canonical correlation variables to form a common feature subspace; the gating filtering based on statistical distance in the common feature subspace includes calculating the Mahalanobis distance between the sliding window sample and the candidate image pixel in the common feature subspace, and weighting and fusing the Mahalanobis distance with the observation geometric angle difference to form a comprehensive matching cost, and retaining only the registration pairs whose comprehensive matching cost is lower than a preset gating threshold.

[0014] In one embodiment, the step of performing label consistency verification using multiple temporal phases of high-resolution optical remote sensing imagery and resolving conflicts by combining spatial label continuity of adjacent centroids along the track includes: For each candidate cluster physical centroid corresponding image pixel, the spectral label of the pixel is obtained at least one time phase before and after the acquisition time of the polar single-photon laser altimetry point cloud data, forming a label time series. If the number of consecutive time phases with the same label in the label time series is less than a preset consistency threshold, the registration result label is determined to be inconsistent and is removed. For centroids that pass the label consistency check, check their physical labels along the track direction and compare them with the labels of their immediate and adjacent centroids. If the label of a centroid is different from the labels of its immediate and adjacent centroids, then adjust the label of that centroid using a neighborhood-weighted smoothing strategy.

[0015] In one embodiment, the output of the surface identification result further includes: extracting the cluster physical centroids corresponding to the overlapping sliding windows identified as inter-ice channel types as channel anchor points; constructing a Gaussian decay weight function with the distance along the track from each channel anchor point to the current sea ice type window to be solved as the independent variable, and performing a distance-weighted moving average on the resampled elevations corresponding to the channel anchor points to obtain the local reference sea level height; for the overlapping sliding windows identified as sea ice types, differencing their resampled elevations with the local reference sea level heights at the corresponding locations to obtain the sea ice freeboard along the track.

[0016] In one embodiment, when there is no cluster physical centroid identified as inter-ice channel type within a certain track interval, the resampled elevation of the overlapping sliding window identified as thin ice type within that interval is extracted as the physical upper limit constraint of the local reference sea level. Combined with the local reference sea level of spatially adjacent intervals, the local reference sea level of that interval is estimated by one-dimensional spline interpolation with physical upper limit constraint. The thin ice is shallow sea ice with a thickness less than a preset threshold.

[0017] The aforementioned surface identification method based on photon point clouds in polar sea ice observation scenarios, in the noise removal stage, extracts the minimum value of the connected graph boundary based on the density gradient field within the first-scale window as the global density segmentation threshold and filters out low-density photons. Then, in the second-scale window, an elevation window is constructed based on the average elevation or modal elevation of the intermediate signal photon set. This removes double-ended outliers that correspond to the high-density isolated tips of transient high-energy mutations, such as those with density fluctuations at the bottom of the low-density tail and those located outside the elevation window with local densities exceeding the preset high-density threshold. This allows diffuse reflection noise and high-density isolated artifact clusters that are misjudged as strong surface echoes by traditional single-threshold methods to be suppressed simultaneously. At the same time, the elevation window constraint avoids the erroneous removal of real surface signal photons by relying solely on the density threshold. This overcomes the limitation of single-boundary density processing in scenarios where diffuse reflection noise rises across the entire polar day and transient high-energy artifacts occur concurrently, making it difficult to simultaneously preserve the fidelity of the noise and remove it. This results in a high-confidence signal set with improved purity. In the training set construction phase, after extracting resampled elevations and generating a sliding classification feature set using the local density structure of photons in the high-confidence signal set as a robust constraint, candidate cluster physical centroids are selected based on the local curvature of the density field and the temporal persistence constraints along the track. Subsequently, the sliding window photon feature set and the spectral texture feature set of optical high-resolution image pixels are projected into a common feature subspace and gated filtering is performed based on the statistical distance and observation geometric angle difference in this subspace. Label consistency verification is performed using multiple temporal phases of optical high-resolution image before and after the acquisition of altimetry data and conflict resolution is performed by combining the spatial label continuity of adjacent centroids along the track. Finally, the multispectral reflectance of optical pixels is input into the ice and snow radiation transmission model to perform physical consistency verification by broadening the theoretical photon return intensity and theoretical elevation. This improves the matching of laser and optical modes from direct mapping of geographic coordinates to robust physical registration in the feature space. Label mistransfer caused by the transient change of optical pixel spectral attributes is blocked by multi-temporal verification and physical consistency verification, thereby obtaining physically anchored physical labels to construct a high-confidence training set. In terms of surface identification, the training set is used as input to the classifier to perform supervised classification on all overlapping sliding windows, and the output includes three types of surface identification results: sea ice, interglacial waterways and thin ice. This results in the final output surface identification results having improved class confidence and consistency in polar sea ice observation scenarios. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0019] Figure 1This is a flowchart of a surface identification method based on photon point clouds in a polar sea ice observation scenario according to an embodiment of the present invention. Figure 2 This is a flowchart illustrating the process of removing two-ended outlier subsets according to an embodiment of the present invention. Figure 3 This is a flowchart of a surface identification method based on photon point clouds in another polar sea ice observation scenario according to an embodiment of the present invention. Figure 4 This is a flowchart illustrating the sea ice freeboard calculation using the identification results, as described in an embodiment of the present invention. Figure 5 This is a schematic diagram comparing the surface identification results of this invention with the ATL07 product. Figure 6 This is an internal structural diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] It should be noted that when a component is referred to as being "fixed to" or "set on" another component, it can be directly on the other component or there may be an intermediate component. When a component is considered to be "connected to" another component, it can be directly connected to the other component or there may be an intermediate component present. The terms "vertical," "horizontal," "upper," "lower," "left," "right," and similar expressions used in this specification are for illustrative purposes only and do not represent the only possible implementation.

[0022] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this invention, "a plurality of" means at least two, such as two, three, etc., unless otherwise explicitly specified.

[0023] In this invention, unless otherwise explicitly specified and limited, "above" or "below" the second feature can mean that the first feature is in direct contact with the second feature, or that the first feature and the second feature are in indirect contact through an intermediate medium. Furthermore, "above," "over," and "on top" the second feature can mean that the first feature is directly above or diagonally above the second feature, or simply indicates that the first feature is at a higher horizontal level than the second feature. "Below," "below," and "under" the second feature can mean that the first feature is directly below or diagonally below the second feature, or simply indicates that the first feature is at a lower horizontal level than the second feature.

[0024] Unless otherwise defined, all technical and scientific terms used in this specification have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used in this specification includes any and all combinations of one or more of the associated listed items.

[0025] The following is combined Figures 1-6 This invention describes a surface identification method based on photon point clouds in polar sea ice observation scenarios.

[0026] like Figure 1 As shown in one embodiment, a surface identification method based on photon point clouds in a polar sea ice observation scenario includes the following steps: Step S110: Obtain polar single-photon laser altimeter cloud data after geophysical correction, calculate the orbital distance and residual elevation of each photon, and construct a two-dimensional photon profile of orbital distance-residual elevation.

[0027] The polar single-photon laser altimeter cloud data originates from raw observation data of spaceborne lidar (such as the ATLAS instrument on the ICESat-2 satellite). After geophysical corrections (such as atmospheric delay correction and tidal correction), large-scale elevation trend terms are removed. In this embodiment, discrete photon points are projected onto a two-dimensional coordinate system consisting of orbital distance (as the horizontal axis) and residual elevation (as the vertical axis), forming a two-dimensional photon profile. The photon distribution pattern in this two-dimensional photon profile intuitively reflects the mixed characteristics of surface undulation and noise background. Signal photons typically exhibit a continuous stripe distribution along the surface contour, while background noise photons show random dispersion characteristics.

[0028] Step S120: Divide the two-dimensional photon profile into first-scale windows along the track direction. Within any first-scale window, calculate the density gradient field of the local density distribution of photons. Extract the minimum value of the boundary of the connected graph corresponding to the change from negative to positive density gradient as the global density segmentation threshold. Filter out photons with local density lower than the global density segmentation threshold to obtain the intermediate signal photon set.

[0029] The first-scale window is designed to adapt to surface features at different scales in the polar environment. The window length can be set according to the satellite's operating speed and sampling rate, for example, set to 1000 photon counts along the orbital direction or a fixed physical distance (e.g., 50 meters). Within the window, the local density in the neighborhood of each photon is calculated, thereby constructing a density gradient field. The density gradient field reflects the rate of change of photon density in space, and its zero-crossing point from negative to positive usually corresponds to the "trough" position of the density distribution, i.e., the boundary between the signal photon accumulation area and the background noise area. In this embodiment, the minimum value of the boundary of the connected graph is defined as the global density segmentation threshold, which can adaptively adjust the segmentation boundary dynamically according to the signal-to-noise ratio within the window. Compared with the traditional fixed threshold method, this method can avoid threshold setting deviations caused by differences in surface reflectivity in polar sea ice observation scenarios, initially filtering out large-scale low-density background noise and retaining potential signal photons.

[0030] Step S130: Define a second-scale window along the track direction for the intermediate signal photon set and perform local analysis. Construct an elevation window with positive and negative preset elevation thresholds based on the average elevation or modal elevation of the intermediate signal photon set. Remove double-ended outliers that correspond to high-density isolated tips of transient high-energy mutations, where the density fluctuations are at the bottom of the low-density tail within the second-scale window and the local density exceeds the preset high-density threshold outside the elevation window. This is done to simultaneously suppress diffuse reflection noise and transient artifacts, and obtain a high-confidence signal set.

[0031] After initial filtering, two types of interference may still remain in the intermediate signal photon set: one is the low-density noise trail formed by diffuse reflection from the ice and snow surface under polar day conditions, and the other is the high-density isolated artifact caused by cosmic rays hitting the detector or electronic transient failures. This embodiment performs local analysis within a more refined second-scale window to identify anomalies at both ends of the density distribution. Photons at the bottom of the low-density trail have a density value lower than the main signal region and belong to diffuse reflection noise; photons at the top of the high-density isolated tip have a density value abnormally higher than the main signal region and are spatially distributed in isolated patches, belonging to transient artifacts. By eliminating these outlier subsets at both ends, this embodiment achieves simultaneous suppression of diffuse reflection noise and transient artifacts, overcoming the limitation of traditional single-threshold methods in simultaneously preserving weak signals and eliminating strong artifacts, thus improving the purity of the signal set.

[0032] Step S140: Using the local photon density structure within the high-confidence signal set as a robust constraint, the resampled elevation is extracted along the track direction using an overlapping sliding window, and a sliding classification feature set with the overlapping sliding window as the statistical unit is generated.

[0033] The local photon density structure within the high-confidence signal set reflects the strength and distribution of surface signals. This embodiment introduces this density structure as a weight constraint to perform weighted averaging or robust estimation of photon elevations within overlapping sliding windows, extracting resampled elevations that represent the true location on the surface. Simultaneously, photon characteristics within each sliding window, such as elevation variance, density peak value, and track slope, are statistically analyzed to construct a multi-dimensional sliding classification feature set, providing input for subsequent clustering analysis.

[0034] Step S150: In the feature space formed by the sliding classification feature set, based on the constraints of the local curvature of the local density field of each sliding window and the time-series persistence along the track, the local density maxima window is selected as the candidate clustering physical centroid.

[0035] In the feature space, real-world surface types (such as sea ice and waterways) typically correspond to local maxima of photon density. This embodiment identifies boundary regions with drastic density changes by calculating the local curvature of the density field. Combined with the constraint of along-track temporal persistence (i.e., real-world surface types should exhibit some continuity along the track direction, rather than isolated jumps), the most representative windows are selected as candidate physical centroids for clustering. These centroids will serve as potential sample centers for the subsequent supervised classification training set.

[0036] Step S160: Register the candidate cluster physical centroid with the spatiotemporally matched optical high-resolution remote sensing image in the real pixel space. Registration includes: projecting the photon feature set of the sliding window and the spectral texture feature set of the optical high-resolution remote sensing image pixels into a common feature subspace, and performing gated filtering based on the statistical distance and observation geometric angle difference within the common feature subspace.

[0037] Because laser altimetry data and optical imagery differ in imaging mechanisms, resolution, and observation time, simple geographic coordinate matching often results in errors. This embodiment first projects the photon feature set (e.g., elevation, density, waveform features) and the spectral texture feature set (e.g., reflectivity, texture mean) of optical image pixels into a common feature subspace, eliminating dimensional and distributional differences between different modal data. In the common feature subspace, the statistical distance (e.g., Mahalanobis distance) between the sliding window sample and candidate image pixels is calculated and weighted in conjunction with observation geometric angle differences (e.g., differences in solar altitude angle and observation zenith angle) to form a comprehensive matching cost. Registration is considered successful only when the comprehensive matching cost is below a preset gating threshold. This gating filter eliminates erroneous registration pairs with excessively large feature space distances or inconsistent observation geometric conditions, resolving the spatial mismatch problem.

[0038] Step S170: For the registration results that have passed the gated filtering, the label consistency is checked using multiple time phases before and after the acquisition time of the polar single-photon laser altimeter cloud data based on the optical high-resolution remote sensing image, and conflict resolution is performed by combining the spatial label continuity of adjacent centroids along the track.

[0039] The spectral properties of optical images can be transiently altered by factors such as cloud cover and instantaneous changes in illumination, leading to unreliable single-phase labels. This embodiment introduces a label consistency check, acquiring spectral labels for registered pixels at multiple time phases before and after the altimeter measurement (e.g., cloudless images within a week before and after), constructing a label time series. If the number of consecutive time phases with the same label in this series is less than a preset threshold, the registration result is deemed inconsistent and discarded. Furthermore, utilizing the spatial continuity of land surface types, labels for adjacent centroids along the track are checked. If a centroid label differs from both before and after it, it is considered an isolated erroneous label and undergoes smoothing correction. This layer of verification mitigates the risk of mislabeling from both temporal and spatial dimensions.

[0040] Step S180: Input the multispectral reflectance of the high-resolution optical remote sensing image pixels into the ice and snow radiative transfer model, and perform physical consistency verification on the registration results verified above using theoretical photon return intensity and theoretical elevation broadening. Assign physical labels to the physical centroids of the clusters that pass the verification to construct a training set.

[0041] Even after passing the first two layers of verification, the registration result still needs to meet physical constraints. This embodiment inputs the multispectral reflectance of the optical image into an ice and snow radiative transfer model (such as the DISORT model) to simulate the theoretical photon return intensity and theoretical elevation broadening under current ice and snow parameters and observation geometry. The simulated theoretical values ​​are compared with the actual observation values ​​from the lidar. If the deviation is within the allowable range, the physical consistency verification is considered successful. This layer of verification ensures the matching of labels and physical attributes, prevents labeling errors caused by spectral confusion, and ultimately constructs a high-confidence training set.

[0042] Step S190: Use the training set as input to the classifier to perform supervised classification on all overlapping sliding windows and output the surface recognition result.

[0043] The process involves training a classifier (such as a support vector machine, random forest, or neural network) using a constructed high-confidence training set to learn the distribution patterns of different land surface types (such as snow-covered sea ice, interglacial channels, and thin ice) in the feature space. The trained classifier then predicts all overlapping sliding windows and outputs the physical label of the entity corresponding to each window, thereby achieving land surface identification in polar sea ice observation scenarios.

[0044] In one embodiment, after forming the two-dimensional photon profile of the photon and before dividing the first-scale window, the method of this embodiment further includes a step of performing a nonlinear adaptive compression transformation on the two-dimensional photon profile. Specifically, a nonlinear adaptive compression transformation is applied to the orbital coordinates in the two-dimensional photon profile, and this transformation is expressed as follows: , In the formula, The coordinates along the track before transformation. These are the transformed coordinates along the track. Based on the scaling ratio, The background photon rate is based on neighborhood statistics. The preset environmental adjustment coefficient is used; furthermore, different base scaling ratios are set for the strong and weak beams of the single-photon altimeter system. .

[0045] Each parameter in the formula has a clear physical meaning and a range of values. As a basic scaling factor, its value is usually set between 0.8 and 1.2, and it is used to normalize the basic scale of the coordinates along the track. The background photon rate is obtained based on neighborhood statistics and reflects the environmental noise level of the current observation area. For example, it can be obtained by calculating the number of photons in the neighborhood of the photon count per unit time. This is an environmental adjustment coefficient used to control the contribution weight of the background photon rate to the degree of coordinate compression. Its value is usually between 0.01 and 0.1, and the specific value can be adjusted according to the actual lighting conditions of the polar environment (such as polar day or polar night).

[0046] This embodiment sets different base scaling ratios for strong beams and weak beams. The reason is that spaceborne single-photon altimeter systems (such as ICESat-2) typically emit two laser beams, one strong and one weak. The strong beam has higher energy and a higher density of returned photon signals, while the weak beam has lower energy and a lower density of returned photon signals. If a uniform scaling ratio is used, the data from the strong beam will be over-compressed in dense regions, while the data from the weak beam will be under-compressed in sparse regions, thus affecting the accuracy of subsequent density calculations. By differentiating the settings, the data from the two laser beams can have relatively consistent density distribution characteristics in the transformed topological space.

[0047] This transformation addresses the problem of sparse and uneven photon distribution in polar sea ice observation scenarios. In these scenarios, surface reflection signals are weak, and the distribution of photon point clouds along the orbital direction is often extremely sparse and uneven. Directly calculating density based on the original orbital coordinates leads to overestimation of photon density in high-background-noise regions and underestimation of signal photon density in low-background-noise regions, causing subsequent density segmentation thresholds to fail. By introducing a nonlinear adaptive compression transformation, the background photon rate is utilized... Dynamically adjusting the coordinate axis scale allows for significant compression of the coordinates along the track in high-background-noise regions and smaller compression in low-background-noise regions, thereby balancing the photon density distribution in different areas. This lays the foundation for calculating the accurate density gradient field and global density segmentation threshold in subsequent step S120. Without this transformation, the density gradient field calculated in subsequent steps will be distorted, making it impossible to accurately extract the minimum values ​​at the boundaries of the connected graph, ultimately leading to a decrease in surface identification accuracy.

[0048] like Figure 2 As shown, in one embodiment, the specific implementation process of removing two-ended outliers includes the following steps: Step S132: Within the second scale window, the density value corresponding to the preset low percentile is used as the low-end cutoff threshold to remove photons with local densities lower than the low-end cutoff threshold, thereby removing diffuse reflection noise.

[0049] Within the second scale window, this embodiment uses a statistical quantile method to define the noise boundary. The preset low percentile is, for example, 5%, which is the 5th percentile of the local density of all photons within the calculation window. This density value is used as the low-end truncation threshold. Photons with a local density below this low-end truncation threshold are considered to be at the bottom of a low-density tail, belonging to diffuse reflection noise caused by the high reflectivity background of polar ice surfaces, and are therefore discarded.

[0050] Step S134: Outside the elevation window, using the density value corresponding to the preset high percentile as the preset high density threshold, photons with local densities higher than the preset high density threshold and spatially isolated patches are removed to eliminate transient high-energy mutation artifacts.

[0051] The preset high percentile is specifically set to, for example, 95%, which is the 95th percentile of the local density of all photons within the calculation window. This density value is used as the high-end truncation threshold (i.e., the preset high-density threshold). For photons with a local density higher than this high-end truncation threshold, it is further determined whether they are outside an elevation window constructed based on the average or modal elevation of the intermediate signal photon set, plus or minus a preset elevation threshold (e.g., plus or minus 10m). Only when a high-density photon is outside this elevation window and spatially appears as an isolated patch, for example, if the area-to-perimeter ratio of its connected domain is lower than a preset shape index threshold, or if the number of photons contained in the high-density cluster is much smaller than the number of photons in the main signal cluster within the window, is it determined to be an isolated tip artifact corresponding to a transient high-energy mutation and is therefore removed. This transient high-energy mutation is usually caused by cosmic rays hitting the detector or by electronic transient failures. By introducing the elevation window constraint, the erroneous removal of real surface signal photons in a high-density state by relying solely on the density threshold is avoided. By employing the aforementioned dual-end removal strategy, this embodiment can simultaneously suppress diffuse reflection noise and transient artifacts at both ends of the density distribution, overcoming the limitation of traditional single-threshold methods in simultaneously preserving weak signals and removing strong artifacts, thereby improving the purity of the signal set.

[0052] In one embodiment, the process of defining a second-scale window for local analysis also includes a step of triggering a pulse denoising mechanism based on the proportion of signal photons.

[0053] The process involves calculating the ratio of signal photons to the total number of photons within the second-scale window to obtain the signal photon proportion. When the signal photon proportion exceeds a preset saturation threshold (e.g., 90% or 95%), detector saturation is determined. Detector saturation often leads to secondary afterpulse peaks following the primary surface reflection peak, creating false elevation signals. In this case, a afterpulse denoising mechanism based on dynamic optimization of residual elevation is triggered: the second derivative is calculated from the histogram of the photon residual elevation distribution within the second-scale window. In the histogram, the primary surface reflection peak corresponds to a high-value region of photon count, and the afterpulse peak corresponds to a secondary high-value region, with a trough between them. The minimum or zero-crossing point of the second derivative can identify the location of this trough. The elevation value corresponding to the minimum point of this trough is established as the adaptive truncation threshold, and photons with residual elevations below this adaptive truncation threshold are identified as afterpulse photons generated by detector saturation and are discarded. This mechanism utilizes the physical characteristic that the post-pulse photon elevation is usually lower than the actual ground surface echo elevation, and dynamically optimizes to adaptively determine the cutoff position, eliminating systematic noise introduced by detector saturation and improving the generation logic of high-confidence signal sets.

[0054] It should be understood that the specific values ​​of the lower percentile, higher percentile, and saturation threshold mentioned above are merely examples, and those skilled in the art can adjust them according to the photon noise level of the actual observation environment. For example, in a polar night environment with low background noise, the lower percentile can be appropriately lowered; in a polar day environment with high reflectivity, the saturation threshold can be appropriately increased. By combining double-ended rejection with a post-pulse denoising mechanism, this embodiment constructs a complete denoising chain, ensuring the data quality for subsequent feature extraction and classification training.

[0055] In one embodiment, the process of extracting resampled elevations along the track direction using overlapping sliding windows with the local photon density structure within the high-confidence signal set as a robust constraint is specifically achieved by introducing a Tukey dual-weight robust estimator.

[0056] In each overlapping sliding window, to suppress the interference of residual noise photons on elevation extraction, this embodiment calculates the resampled elevation according to the following robust density-weighted formula. : , In the formula, For the first overlapping sliding window The residual elevation of a high-confidence signal photon. This represents the local density value of the photon. This represents the total number of high-confidence signal photons within the overlapping sliding window. The Tukey robust weighting factor, which is iteratively updated based on photon elevation residuals, is calculated using a double-squared function: , in, For the first The residual of each photon relative to the weighted median elevation within the window This is a robust scaling estimate of the residuals (such as the median absolute deviation, MAD). To adjust the constant, its value is typically between 4 and 6. The physical meaning of this weighting factor is: when the photon elevation residual is small, the weight is close to 1, and the photons are fully preserved; when the residual exceeds a preset threshold... At this point, the weight is truncated to 0, and the photon is completely removed. The weight factor is then updated iteratively. This embodiment adaptively reduces the contribution of outliers (such as incompletely filtered noise photons or anomalous echoes caused by surface roughness) to the resampled elevation, thereby achieving a robust estimate of the true location of the ground surface using local density structure. Compared to simple arithmetic mean or median filtering, the Tukey dual-weight estimator suppresses extreme value interference while preserving surface details, thus improving the accuracy of the resampled elevation.

[0057] In one embodiment, the construction of the sliding classification feature set also includes extracting track slope features and residual elevation skewness features.

[0058] The track slope feature is extracted by first-order forward difference of the resampled elevation of adjacent overlapping sliding windows, and its calculation formula is as follows: , In the formula, This is the index of the overlapping sliding window along the track direction. For the first Resampled elevation of overlapping sliding windows.

[0059] This feature reflects the degree of surface inclination along the track direction. Physically, sea ice surfaces typically exhibit dramatic slope changes due to snow cover and ice ridges; while interglacial channels, as linear open water bodies, tend to have flatter surfaces along the track direction, with smaller and more stable slope characteristics. Therefore, the track slope feature can enhance the morphological separability between sea ice and channels.

[0060] The residual elevation skewness feature is extracted from the third-order central moments of the residual elevation of high-confidence signal photons within the overlapping sliding window. The calculation formula is as follows: , In the formula, The residual elevation of photons within the window. The average elevation of photons within the window. The standard deviation is denoted as σ. Skewness measures the asymmetry of elevation distribution. Physically, the elevation distribution of sea ice surfaces is often positively skewed (long tails biased towards the side with higher elevation) due to the presence of ice ridges, while the elevation distribution of interglacial channels is approximately symmetrical or slightly negatively skewed. This embodiment introduces residual elevation skewness features to characterize the microstructural differences of different surface types from a statistical distribution perspective, enriching the dimensions of the sliding classification feature set and providing more discriminative input features for subsequent supervised classification. It should be understood that the above feature calculation method is only an example. Without departing from the concept of this invention, other statistical features such as second-order central moments (variance) or fourth-order central moments (kurtosis) can also be used to describe surface characteristics.

[0061] In one embodiment, the specific process of projecting the photon feature set of the sliding window and the spectral texture feature set of optical high-resolution remote sensing image pixels into a common feature subspace, and the specific algorithm for gating filtering based on statistical distance within the common feature subspace are detailed below.

[0062] Since the lidar photon feature set (including elevation, density, waveform features, etc.) and the optical image spectral texture feature set (including reflectivity, texture mean, variance, etc.) are two heterogeneous data modes, direct distance calculation results in differences in dimensions and distribution. This embodiment uses canonical correlation analysis to process the two feature sets. Canonical correlation analysis calculates the covariance matrix between the two sets of variables and finds the linear combination coefficient that maximizes the correlation between the two sets of variables, thereby extracting the canonical correlation variables. The space spanned by these canonical correlation variables is the common feature subspace, in which the lidar photon features and optical image features achieve maximum alignment and coupling.

[0063] Within the common feature subspace, this embodiment calculates the Mahalanobis distance between the sliding window samples and candidate image pixels. Mahalanobis distance considers the covariance structure between each feature dimension, eliminating distance distortion caused by feature correlation. Compared to Euclidean distance, it is more robust to outliers. The formula for calculating Mahalanobis distance is: , In the formula, Let be the feature vector of the sliding window samples in the common feature subspace. For the feature vectors of candidate image pixels, Let be the sample covariance matrix.

[0064] To comprehensively consider feature similarity and observation geometry, this embodiment weights and fuses Mahalanobis distance and observation geometric angle difference into a comprehensive matching cost. Observation geometric angle difference includes differences in solar altitude angle and observation zenith angle; excessively large geometric angle differences can lead to instability in spectral features. The formula for calculating the comprehensive matching cost can be expressed as: , In the formula, and The preset weighting coefficients, This is the normalized value of the observed geometric angle difference. Only when the overall matching cost... Only when the value falls below a preset gating threshold is the registration pair deemed to have passed the gating filter. This filtering mechanism eliminates erroneous registration pairs with excessively large feature space distances or severely mismatched observation geometry, thus resolving the spatial mismatch problem caused by relying solely on geographic coordinate matching.

[0065] like Figure 3 As shown, in one embodiment, the specific process of performing label consistency verification using multiple temporal phases of high-resolution optical remote sensing imagery and resolving conflicts by combining the spatial label continuity of adjacent centroids along the track includes the following steps: Step S172: For each image pixel corresponding to the physical centroid of a candidate cluster, obtain the spectral label of that pixel at least one time phase before and after the acquisition time of the polar single-photon laser altimeter point cloud data, forming a label time series. If the number of consecutive time phases with the same label in the label time series is less than a preset consistency threshold, the registration result is determined to be inconsistent with the label and is removed.

[0066] In polar environments, varying lighting conditions, cloud shadows, flares, or transient sea ice drift can lead to transient errors in the spectral labels of single-temporal optical images. This embodiment introduces a time-dimensional verification mechanism. For example, multiple cloud-free images are acquired within a week before and after the altimeter measurement, and a spectral classifier is used to extract labels for each temporal phase (e.g., snow-covered sea ice, waterways, thin ice). After constructing the label time series, the number of consecutive occurrences of the dominant label in the series is counted. If a registered pixel is labeled as a waterway at the altimeter measurement but is labeled as sea ice in both preceding and following phases, and the number of consecutive phases is less than a preset threshold (e.g., 3 phases), then the registration result is determined to be affected by transient interference, the labels are inconsistent, and it is discarded. This step utilizes the temporal stability of land cover types to prevent mislabeling caused by transient changes in spectral attributes.

[0067] Step S174: For centroids that pass the label consistency check, check their physical labels along the track direction and compare them with the labels of their adjacent centroids. If the label of a centroid is different from the labels of its adjacent centroids, adjust the label of the centroid using a neighborhood weighted smoothing strategy.

[0068] Land surface types typically exhibit spatial continuity, and isolated labels often indicate classification errors. This embodiment performs a spatial continuity check on the centroids that have passed the aforementioned verification. The centroid sequence is traversed along the track direction; if a certain centroid is found... The label is "waterway," and its adjacent centroids are... and If all the tags are sea ice, then it is determined that... For isolated erroneous labels, a neighborhood-weighted smoothing strategy is used for adjustment. The specific logic is as follows: calculate... and , In the feature space, the distance is weighted by the inverse of the distance, and the label features of adjacent centroids are weighted and assigned accordingly. The strategy leverages the spatial continuity constraint of sea ice and waterways to eliminate isolated erroneous labels and improve the spatial consistency of the training set samples.

[0069] like Figure 4 As shown, in one embodiment, after outputting the surface identification result, the method of this embodiment further includes a step of calculating the sea ice freeboard using the identification result. Sea ice freeboard refers to the height of sea ice above the water surface and is a key parameter for calculating sea ice thickness and estimating sea ice mass balance. This embodiment calculates the sea ice freeboard through the following steps: Step S410: Extract the cluster physical centroid corresponding to the overlapping sliding window identified as an inter-ice waterway type as the waterway anchor point.

[0070] After supervised classification, the surface identification results include category labels for each overlapping sliding window (e.g., sea ice, interglacial channels, thin ice, etc.). This embodiment filters out windows labeled "interglacial channels" and extracts the corresponding cluster physical centroids as channel anchors. Channel anchors represent the actual physical observation locations of sea surface height, and their resampled elevation is theoretically equal to the local sea surface height. It should be understood that the extraction criteria for channel anchors are not limited to a single category label; they can also be combined with the photon density distribution pattern within the window (e.g., low variance, low skewness) for secondary verification to ensure that the anchors are located in actual open water or at the center of newly formed interglacial channels, thus eliminating misidentified anchors caused by cloud shadow interference.

[0071] Step S420: Construct a Gaussian attenuation weight function with the distance along the track from each channel anchor point to the current sea ice type window to be solved as the independent variable, and perform a distance-weighted moving average on the resampled elevation corresponding to the channel anchor point to obtain the local reference sea surface height.

[0072] Because sea level varies slightly in space (e.g., tilted by tides and geostrophic currents), the average elevation of all channel anchor points cannot be simply taken as the reference sea level. This embodiment introduces a Gaussian attenuation weighting function, the expression of which is: , In the formula, This represents the distance along the track from the waterway anchor point to the current window to be solved. ξ This is the attenuation factor. ξThe setting determines the spatial smoothing scale of the reference sea surface, and its value is typically set between 500 meters and 2000 meters. Smaller... ξ The value can capture high-frequency sea surface fluctuations, but is easily affected by individual abnormal anchor points; larger values... ξ The value smoothing effect is stronger, but it may ignore local sea surface details. By using a distance-weighted moving average, waterway anchors that are closer to the current window are given greater weight, thereby estimating a local reference sea surface height that is closer to the actual physical state.

[0073] Step S430: For the overlapping sliding window identified as sea ice type, the resampled elevation is differentially compared with the local reference sea level at the corresponding location to obtain the sea ice freeboard along the track.

[0074] Among them, sea ice freeboard The calculation formula is: ,in For the resampling elevation of the sea ice window, This is the local reference sea level height estimated in step S420. This difference operation eliminates the influence of the sea level height background field and extracts the height information of sea ice relative to the sea level.

[0075] In one embodiment, when there is no cluster physical centroid identified as inter-ice channel type within a certain track interval, the resampled elevation of the overlapping sliding window identified as thin ice type within that interval is extracted as the physical upper limit constraint of the local reference sea level. Combined with the local reference sea level of spatially adjacent intervals, the local reference sea level of that interval is estimated by one-dimensional spline interpolation with physical upper limit constraint. Thin ice is shallow sea ice with a sea ice thickness less than a preset threshold.

[0076] In polar regions with thick ice cover, interglacial channels are sparse or even completely absent, preventing the extraction of usable channel anchor points in step S410 and resulting in a gap in the estimation of local reference sea level. To address this, this embodiment introduces "thin ice type" as a physical upper limit constraint. Thin ice (such as Nile ice and nascent ice) has a smaller thickness and, due to buoyancy, its surface elevation is slightly higher than the sea surface, but closer to sea level than thick ice. Therefore, the resampled elevation of thin ice can serve as the physical upper limit for the local reference sea level. This embodiment defines thin ice as shallow sea ice with a thickness less than a preset threshold (e.g., 0.3 meters).

[0077] In the interpolation process, the resampled elevations of all thin ice windows within the thick ice-covered area are first extracted as constraints. Then, combined with the known local reference sea level heights of spatially adjacent intervals (with channel anchor points), one-dimensional spline interpolation with physical upper limit constraints is performed. When fitting the spline curve, this interpolation algorithm enforces a constraint that the interpolation result must not exceed the physical upper limit set for the thin ice elevation. This mechanism utilizes the static equilibrium relationship between the thin ice and the sea surface, filling the gap in sea level height estimation through physical constraints in the absence of direct physical anchor points, and solving the systematic error problem caused by the lack of a reference surface in the calculation of sea ice freeboard in thick ice-covered areas.

[0078] In one specific embodiment, to verify the effectiveness of the surface identification method based on photon point clouds in the polar sea ice observation scenario provided by the present invention in practical applications, this embodiment selects the Beaufort Sea region in the Arctic as the study area and uses the raw photon point cloud data acquired by the ATLAS lidar carried by the ICESat-2 satellite during the summer polar day to conduct a sea ice freeboard inversion experiment.

[0079] The experimental data selected was GT1R strong beam data from the ATLAS Level 2A product (ATL03), with an orbital distance coverage range of 75° to 80° North latitude. First, step S110 was executed to acquire geophysically corrected polar single-photon laser altimeter cloud data, calculate the orbital distance and residual elevation of each photon, and construct a two-dimensional photon profile of orbital distance versus residual elevation. In this embodiment, the geophysical corrections included atmospheric delay correction, ocean tidal correction, and reverse pressure correction; the residual elevation was the relative elevation value after deducting the mean sea level height.

[0080] Subsequently, a nonlinear adaptive compression transformation is applied to the orbital coordinates in the two-dimensional photon profile. In this embodiment, the base scaling factor a0 is set to 1.0 for strong beams, and the environmental adjustment coefficient μ is set to 0.05. The background photon rate Rbg is obtained by statistically analyzing the photon count rate within every 100-meter interval along the orbital direction. This transformation balances the photon density distribution between high and low background noise regions.

[0081] Next, step S120 is executed, dividing the two-dimensional photon profile into a first-scale window along the track direction. In this embodiment, the length of the first-scale window is set to 1000 photon counts along the track direction. Within the window, the density gradient field of the local photon density distribution is calculated, and the minimum value of the connected graph boundary corresponding to the transition from negative to positive density gradient is extracted as the global density segmentation threshold. In mixed ice regions, this adaptive threshold, compared to the traditional fixed threshold method, can retain more weak signal photons.

[0082] For the intermediate signal photon set, step S130 is executed to define a second-scale window along the track direction for local analysis. In this embodiment, the physical length of the second-scale window is set to 50 meters. When removing outliers at both ends, the lower percentile is set to 5%, and the higher percentile is set to 95%. For windows where the signal photon proportion is higher than 90%, detector saturation is determined, triggering a post-pulse denoising mechanism. The second derivative is used to identify the minimum point of the trough as an adaptive truncation threshold. This step removes the low-density background noise introduced by strong diffuse reflection from ice and snow surfaces under polar day conditions, as well as the high-density isolated artifacts caused by cosmic rays.

[0083] Then, step S140 is executed, using the local photon density structure within the high-confidence signal set as a robust constraint, to extract resampled elevations along the track direction using an overlapping sliding window. In this embodiment, the step size of the overlapping sliding window is set to 10 meters, and the window width is set to 30 meters. A Tukey dual-weight robust estimator is introduced to calculate the resampled elevations, and track-side slope features and residual elevation skewness features are extracted to construct a sliding classification feature set.

[0084] Within the feature space, steps S150 to S180 are executed to screen candidate cluster physical centroids and perform three-layer verification. In this embodiment, the high-resolution optical remote sensing image used for registration is Sentinel-2 MSI data. In the gated filtering stage, canonical correlation analysis (CCA) is used to extract the common feature subspace, and the weighted fusion cost of Mahalanobis distance and observation geometric angle difference is calculated. In the label consistency verification stage, multi-temporal images within one week before and after the altimeter are acquired to construct a label time series, and the consistency threshold is set to three consecutive temporal phases. In the physical consistency verification stage, the multispectral reflectance of the optical image is input into the DISORT snow radiative transfer model to simulate the theoretical photon return intensity.

[0085] Finally, step S190 is executed, where the constructed training set is input into the random forest classifier to perform supervised classification on all overlapping sliding windows, outputting the surface identification results. The identification results include three types: snow-covered sea ice, interglacial waterways, and thin ice.

[0086] Based on the surface identification results, sea ice freeboard calculation is performed. The centroid of the window identified as an inter-ice channel type is extracted as the channel anchor point. A Gaussian attenuation weight function is constructed, with the attenuation factor σ set to 1000 meters. The local reference sea level height is obtained by distance-weighted moving average of the channel anchor point elevation. For the thick ice-covered area without channel anchor points, the elevation of the thin ice type window is extracted as a physical upper limit constraint, and one-dimensional spline interpolation with physical upper limit constraint is performed.

[0087] To quantitatively evaluate the technical effectiveness of the method of this invention, the sea ice freeboard results obtained by inversion can be compared and verified with the ATL07 product officially released by ICESat-2, using airborne high-precision laser altimetry data as the true reference. In thin ice areas (sea ice thickness less than 0.3 meters), due to the nonlinear adaptive compression transformation and Tukey's dual-weight robust estimator's preservation of weak signal photons and noise suppression, the root mean square error of the sea ice freeboard inverted by the method of this invention is lower than that of the ATL07 product. In thick ice-covered areas (sea ice thickness greater than 1.5 meters), the one-dimensional spline interpolation method with physical upper limit constraints compensates for the reference sea surface estimation bias caused by the lack of channel anchor points, and the sea ice freeboard inversion accuracy is also better than that of the ATL07 product. In inter-ice channel areas, the three-layer verification mechanism avoids label mistransfer, improves the purity of the training set, and the recognition recall rate is better than that of the traditional density threshold method. The above comparative verification shows that the present invention has good accuracy and robustness in surface identification and sea ice freeboard inversion in polar sea ice observation scenarios.

[0088] like Figure 5 As shown, to more intuitively demonstrate the technical advantages of the method of the present invention, this embodiment presents a visual comparison of the surface identification results of the method of the present invention (left column) and the ATL07 product officially released by ICESat-2 (right column) in two different polar sea ice observation scenarios. In the figure, blue represents sea ice, red represents interglacial channels, and yellow represents thin ice.

[0089] from Figure 5 The comparison results clearly show that: First, the ATL07 product (right column) can only distinguish between two types: sea ice and interglacial channels (see right column illustration), and cannot identify thin ice; while the method of the present invention (left column) can accurately output identification results for three types, including snow-covered sea ice, interglacial channels and thin ice, with finer classification granularity.

[0090] Second, in the uplink scenario, the method of the present invention accurately identified thin ice (yellow trajectory segment) at the edge of the inter-ice channel, while ATL07 missed or misjudged it as sea ice / inter-ice channel; in the downlink scenario, the method of the present invention identified a large area of ​​continuous thin ice, which is more consistent with the actual physical distribution characteristics of sea ice (such as the ice edge zone).

[0091] The comparative results fully demonstrate that the method of the present invention effectively overcomes the problems of single classification categories and easy misjudgment in polar sea ice observation scenarios by introducing elevation window constraints to eliminate double-ended outlier subsets and constructing training sets by combining multi-layer physical verification. It significantly improves the accuracy of surface identification, especially achieving effective identification of the key type of thin ice.

[0092] Figure 6This example illustrates a schematic diagram of the physical structure of an electronic device, which can be a smart terminal. Its internal structure diagram can be as follows: Figure 6 As shown, the electronic device includes a processor, memory, and a network interface connected via a system bus. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The network interface is used to communicate with external terminals via a network connection. When the computer program is executed by the processor, it implements the surface identification method based on photon point clouds in the polar sea ice observation scenario of any of the above embodiments.

[0093] Those skilled in the art will understand that Figure 6 The structure shown is merely a block diagram of a portion of the structure related to the present invention and does not constitute a limitation on the electronic device to which the present invention is applied. A specific electronic device may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0094] On the other hand, the present invention also provides a computer storage medium storing a computer program, which, when executed by a processor, implements the surface identification method based on photon point clouds in the polar sea ice observation scenario of any of the above embodiments.

[0095] In another aspect, a computer program product or computer program is provided, which includes computer instructions stored in a computer-readable storage medium. A processor of an electronic device reads the computer instructions from the computer-readable storage medium, and when the processor executes the computer instructions, it implements the surface identification method based on photon point clouds in polar sea ice observation scenarios according to any of the above embodiments.

[0096] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. This computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided by this invention can include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory may include random access memory (RAM) or external cache memory.

[0097] By way of illustration and not limitation, RAM is available in a variety of forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0098] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0099] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the appended claims.

Claims

1. A surface identification method based on photon point clouds in polar sea ice observation scenarios, characterized in that, include: Acquire geophysically corrected polar single-photon laser altimetry point cloud data, calculate the orbital distance and residual elevation of each photon, and construct a two-dimensional photon profile of orbital distance and residual elevation; The two-dimensional photon profile is divided into first-scale windows along the track direction. Within any first-scale window, the density gradient field of the local density distribution of photons is calculated. The minimum value of the boundary of the connected graph corresponding to the change from negative to positive density gradient is extracted as the global density segmentation threshold. Photons with local density lower than the global density segmentation threshold are filtered out to obtain the intermediate signal photon set. For the intermediate signal photon set, a second-scale window is defined along the track direction to perform local analysis. An elevation window with positive and negative preset elevation thresholds is constructed based on the average elevation or modal elevation of the intermediate signal photon set. Double-ended outliers corresponding to high-density isolated tips of transient high-energy mutations are removed from the second-scale window, where the density fluctuation is at the bottom of the low-density tail and the local density exceeds the preset high-density threshold outside the elevation window. This is done to simultaneously suppress diffuse reflection noise and transient artifacts, and obtain a high-confidence signal set. Using the local photon density structure within the high-confidence signal set as a robust constraint, resampled elevations are extracted along the track direction using overlapping sliding windows, and a sliding classification feature set with the overlapping sliding windows as statistical units is generated. Within the feature space formed by the sliding classification feature set, based on the constraints of the local curvature of the local density field of each sliding window and the temporal persistence along the track, the local density maxima window is selected as the candidate cluster physical centroid. The candidate cluster physical centroid is registered with the spatiotemporally matched optical high-resolution remote sensing image in the real pixel space. The registration includes: projecting the photon feature set of the sliding window and the spectral texture feature set of the optical high-resolution remote sensing image pixels into a common feature subspace, and performing gated filtering based on the statistical distance and observation geometric angle difference in the common feature subspace. For the registration results that pass through the gated filtering, the label consistency is checked using the optical high-resolution remote sensing image at multiple time phases before and after the acquisition time of the polar single-photon laser altimeter cloud data, and conflict resolution is performed by combining the spatial label continuity of adjacent centroids along the track. The multispectral reflectance of the high-resolution optical remote sensing image pixels is input into the ice and snow radiative transfer model. The physical consistency of the registration results verified above is checked by using theoretical photon return intensity and theoretical elevation broadening. The physical centroids of the clusters that pass the verification are assigned physical labels to the entities to construct a training set. The training set is used as input to the classifier to perform supervised classification on all the overlapping sliding windows, and the surface identification results are output, which include three types: sea ice, interglacial waterways, and thin ice.

2. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The two-dimensional photon profile is divided into a first-scale window, which includes: A nonlinear adaptive compression transformation is applied to the orbital coordinates in the two-dimensional photon profile, and this transformation is expressed as: ， In the formula, The coordinates along the track before transformation. These are the transformed coordinates along the track. Based on the scaling ratio, The background photon rate is based on neighborhood statistics. The preset environmental adjustment coefficient is used; and, the basic scaling ratio is set differently for the strong beam and weak beam of the single-photon altimeter system. .

3. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The exclusion of the double-ended outlier subset, which includes high-density isolated tips corresponding to transient high-energy mutations, where density fluctuations within the second-scale window are at the bottom of a low-density tail and outside the elevation window with local densities exceeding a preset high-density threshold, includes: Within the second scale window, the density value corresponding to the preset low percentile is used as the low-end cutoff threshold to remove photons with local densities lower than the low-end cutoff threshold, thereby eliminating diffuse reflection noise. Outside the elevation window, the density value corresponding to the preset high percentile is used as the preset high density threshold to remove photons with local densities higher than the preset high density threshold and spatially isolated patches, thus eliminating transient high-energy mutation artifacts.

4. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The step of extracting resampled elevations along the orbital direction using an overlapping sliding window, with the photon local density structure within the high-confidence signal set as a robust constraint, includes: Within each of the overlapping sliding windows, a Tukey dual-weight robust estimator is introduced to calculate the resampled elevation according to the following robust density-weighted formula. : ， In the formula, For the first overlapping sliding window The residual elevation of a high-confidence signal photon. This represents the local density value of the photon. This refers to the Tukey robust weighting factor, which is updated iteratively based on photon elevation residuals. This represents the total number of high-confidence signal photons within the overlapping sliding window.

5. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The sliding classification feature set also includes: The track slope features are extracted by the first-order forward difference of the resampled elevations of adjacent overlapping sliding windows, and the residual elevation skewness features are extracted by the third-order central moments of the residual elevations of high-confidence signal photons within the overlapping sliding windows.

6. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, During the local analysis performed within the defined second-scale window, the ratio of the number of signal photons to the total number of photons within the second-scale window is calculated as the signal photon proportion. When the signal photon proportion is higher than a preset saturation threshold, the detector is determined to be saturated, and a post-pulse denoising mechanism based on dynamic optimization of residual elevation is triggered within the second-scale window: the second derivative of the histogram of photon residual elevation distribution within the second-scale window is calculated, the minimum point of the valley between the primary surface reflection peak and the secondary post-pulse peak is identified, the elevation value corresponding to the minimum point of the valley is established as the adaptive truncation threshold, and photons with residual elevation lower than the adaptive truncation threshold are determined to be post-pulse photons generated by detector saturation and are removed.

7. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The step of projecting the photon feature set of the sliding window and the spectral texture feature set of the optical high-resolution remote sensing image pixels onto a common feature subspace includes: Canonical correlation analysis is performed on the photon feature set and the spectral texture feature set respectively to extract canonical correlation variables to form a common feature subspace. Gating filtering is performed based on the statistical distance in the common feature subspace, including calculating the Mahalanobis distance between the sliding window sample and the candidate image pixel in the common feature subspace, and weighting and fusing the Mahalanobis distance with the observation geometric angle difference to form a comprehensive matching cost. Only registration pairs with a comprehensive matching cost lower than a preset gating threshold are retained.

8. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, Label consistency verification is performed using multiple temporal phases of high-resolution optical remote sensing imagery, and conflict resolution is achieved by combining spatial label continuity of adjacent centroids along the track, including: For each candidate cluster physical centroid corresponding image pixel, the spectral label of the pixel is obtained at least one time phase before and after the acquisition time of the polar single-photon laser altimetry point cloud data, forming a label time series. If the number of consecutive time phases with the same label in the label time series is less than a preset consistency threshold, the registration result label is determined to be inconsistent and is removed. For centroids that pass the label consistency check, check their physical labels along the track direction and compare them with the labels of their immediate and adjacent centroids. If the label of a centroid is different from the labels of its immediate and adjacent centroids, then adjust the label of that centroid using a neighborhood-weighted smoothing strategy.

9. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 1, characterized in that, The output surface identification result then includes: Extract the cluster physical centroids corresponding to the overlapping sliding windows identified as interglacial waterways as waterway anchor points; A Gaussian decay weighting function is constructed with the distance along the track from each of the aforementioned channel anchor points to the current sea ice type window to be solved as the independent variable. The resampled elevation corresponding to the channel anchor points is then subjected to a distance-weighted moving average to obtain the local reference sea surface height. For overlapping sliding windows identified as sea ice types, their resampled elevations are differiated from the local reference sea level height at the corresponding location to obtain the sea ice freeboard along the track.

10. The surface identification method based on photon point clouds in polar sea ice observation scenarios according to claim 9, characterized in that, When there is no cluster physical centroid identified as inter-ice channel type within a certain track interval, the resampled elevation of the overlapping sliding window identified as thin ice type within that interval is extracted as the physical upper limit constraint of the local reference sea level. Combined with the local reference sea level of spatially adjacent intervals, the local reference sea level of that interval is estimated by one-dimensional spline interpolation with physical upper limit constraint. The thin ice is shallow sea ice with a thickness less than a preset threshold.

Images