Airborne glacier tomography processing method based on refraction correction and coherent stacking

Abstract

The application relates to the technical field of radar remote sensing and cryosphere detection, and provides an airborne glacier tomography processing method based on refraction correction and coherence stacking. The method introduces refraction path correction, improves the positioning accuracy of an under-ice scatterer in three-dimensional space, and avoids systematic errors caused by a traditional straight-line propagation assumption; coherence weighting tomographic focusing is used, the robustness of tomographic imaging under a low signal-to-noise ratio condition is improved, and the continuous interface under the ice is stably extracted; the method is not dependent on a specific orbit configuration or a specific implementation algorithm, has good universality and scalability, and is suitable for various glacier radar detection scenes.

Description

Technical Field

[0001] This application relates to the fields of radar remote sensing and cryosphere detection technology, and in particular to an airborne glacier tomography processing method based on refraction correction and coherent stacking. Background Technology

[0002] Glacier thickness and its internal structure are crucial parameters for studying glacier dynamics, assessing changes in glacier volume, and monitoring the cryosphere environment. High-altitude glaciers are widely distributed in high-altitude regions such as the Qinghai-Tibet Plateau, playing a vital role in regional water resource regulation and supply. In recent years, influenced by global warming, the rate of glacier ablation has accelerated significantly, drawing widespread attention to changes in glacier water volume and its potential environmental risks. However, due to limitations in methods for acquiring glacier thickness data, accurate estimations of glacier volume and its changes still face considerable uncertainty.

[0003] Synthetic Aperture Radar (SAR), as an active microwave remote sensing technology, possesses all-weather, all-day imaging capabilities and has been widely used in glacier remote sensing. Compared to indirect model estimation methods based on optical imagery, topographic data, or ice flow velocity, radar detection can directly acquire scattering information from within and beneath the ice, helping to improve the reliability of ice thickness and volume estimation. Traditional airborne or ground-based ice-penetrating radars typically only provide one-dimensional or two-dimensional reflection profiles along the survey line. In large glacier areas, interpolation methods are still needed to infer unmeasured areas, further introducing uncertainty.

[0004] In recent years, Tomographic Synthetic Aperture Radar (TomoSAR) technology has enabled three-dimensional imaging of the internal scattering structure of glaciers through multi-track or multi-view observations, providing a new technical means for glacier thickness inversion and internal structure analysis. However, in practical applications, the complex structure of the glacier medium causes radar waves to refract at the air-ice interface, altering their propagation path and velocity. If this is not adequately considered during imaging, it can easily introduce spatial positioning errors for subglacial targets. Furthermore, with increasing detection depth, the signal-to-noise ratio of the bedrock echo signal under the ice decreases significantly, making tomographic imaging results susceptible to noise and spurious scattering, affecting the stability and reliability of the three-dimensional focusing results. Summary of the Invention

[0005] In view of this, this application provides an airborne glacier tomography processing method based on refraction correction and coherent stacking to solve the problem that the recognition accuracy and stability of the existing glacier radar tomography method need to be improved.

[0006] A first aspect of this application provides an airborne glacier tomography processing method based on refraction correction and coherence stacking, comprising:

[0007] The multi-track SAR echo data acquired by tomographic synthetic aperture radar (SAR) is processed for imaging, and the imaging results are calibrated to obtain a calibrated complex radar image.

[0008] A three-dimensional imaging grid is constructed based on the digital elevation model of the ice surface for a predetermined glacier region; the predetermined glacier region includes an area covering a predetermined depth range below the ice surface.

[0009] Three-dimensional tomography focusing is performed using the calibrated complex radar images in the three-dimensional imaging grid to generate three-dimensional focused images of each SAR track in the preset area of ​​the glacier;

[0010] Phase coherence superposition of the three-dimensional focused images of each SAR track within the target imaging layer in the preset glacier region yields the three-dimensional scattering intensity distribution of the target in the preset glacier region; wherein, the target imaging layer is the layer in the preset glacier region with the preset target depth;

[0011] Based on the three-dimensional scattering intensity distribution of the target, spatial continuity analysis of the scatterer under the ice is performed to obtain imaging results.

[0012] A second aspect of this application provides an airborne glacier tomography processing device based on refraction correction and coherence stacking, comprising:

[0013] The pre-imaging module is configured to perform imaging processing on the multi-track SAR echo data acquired by tomographic synthetic aperture radar SAR, and to calibrate the imaging results to obtain a calibrated complex radar image.

[0014] The construction module is configured to build a 3D imaging mesh of a preset glacier region based on the digital elevation model of the ice surface; wherein the preset glacier region includes an area covering a predetermined depth range below the ice surface;

[0015] The focusing module is configured to perform three-dimensional tomographic focusing in the three-dimensional imaging grid using the calibrated complex radar images to generate three-dimensional focused images of each orbit SAR of the preset glacier region;

[0016] The overlay module is configured to perform phase coherence overlay on the three-dimensional focused images of each SAR track within the target imaging layer in the preset region of the glacier to obtain the three-dimensional scattering intensity distribution of the target in the preset region of the glacier; wherein, the target imaging layer is a layer in the preset region of the glacier with a preset target depth;

[0017] The imaging module is configured to perform spatial continuity analysis on the sub-ice scatterer based on the three-dimensional scattering intensity distribution of the target, and obtain imaging results.

[0018] The beneficial effects of this application embodiment compared with the prior art are as follows: By introducing refraction path correction, this application embodiment improves the positioning accuracy of subglacial scatterers in three-dimensional space and avoids the systematic errors caused by the traditional linear propagation assumption; by using coherence-weighted tomographic focusing, it improves the robustness of tomographic imaging under low signal-to-noise ratio conditions, which is conducive to the stable extraction of subglacial continuous interfaces; the method does not depend on specific orbital configurations or specific implementation algorithms, has good versatility and scalability, and is suitable for various glacier radar detection scenarios. Attached Figure Description

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

[0020] Figure 1 This is a schematic flowchart of an airborne glacier tomography processing method based on refraction correction and coherent stacking provided in an embodiment of this application.

[0021] Figure 2 This is a schematic flowchart of another airborne glacier tomography processing method based on refraction correction and coherent stacking provided in the embodiments of this application.

[0022] Figure 3 This is a schematic diagram of an airborne glacier tomography processing device based on refraction correction and coherent stacking, provided in an embodiment of this application.

[0023] Figure 4 This is a schematic diagram of the electronic device provided in the embodiments of this application. Detailed Implementation

[0024] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0025] The following will describe in detail, with reference to the accompanying drawings, an airborne glacier tomography processing method and apparatus based on refraction correction and coherent stacking according to embodiments of this application.

[0026] As mentioned above, in practical applications, radar waves are refracted at the air-ice interface, which may alter their propagation path and speed. Furthermore, the signal-to-noise ratio of the bedrock echo signal under ice decreases significantly with increasing depth. These factors all affect the stability and reliability of the 3D focusing results. Specifically, the glacier TomoSAR imaging method in related technologies typically assumes that the radar wave propagates in a homogeneous medium during focusing and then performs simplified corrections for refraction effects after imaging. While this approach is computationally efficient, it struggles to achieve precise phase history matching, limiting the accuracy of target localization under ice and the effective detection capability of deep-seated scatterers.

[0027] In other words, the glacier radar tomography method in related technologies has the following problems: radar waves are refracted at the air-ice medium interface. If the linear propagation assumption is still used for tomographic focusing, it is easy to cause systematic shifts in the depth and lateral positions of the scattering bodies under the ice. In the interior of the glacier and the subglacial region, the signal-to-noise ratio of the radar echo signal is low, and the multi-track tomography results are easily affected by noise and false scattering, making it difficult to stably identify the bedrock under the ice.

[0028] In view of this, this application provides an airborne glacier tomography processing method based on refraction correction and coherent stacking. This method can reasonably consider the influence of medium refraction during tomography and improve the three-dimensional focusing robustness of glacier radar imaging under low signal-to-noise ratio conditions, so as to meet the application requirements of accurate detection of ice thickness and internal structure in complex mountain glacier areas.

[0029] Figure 1 This is a schematic flowchart of an airborne glacier tomography processing method based on refraction correction and coherent stacking, provided in an embodiment of this application. Figure 1 As shown, the method includes the following steps:

[0030] In step S101, the multi-track SAR echo data acquired by tomographic SAR is processed for imaging, and the imaging results are calibrated to obtain a calibrated complex radar image.

[0031] In step S102, a three-dimensional imaging grid of the preset area of ​​the glacier is constructed based on the digital elevation model of the ice surface.

[0032] The glacier pre-defined area includes the region covering a predetermined depth below the ice surface.

[0033] In step S103, three-dimensional tomography focusing is performed on the calibrated complex radar images in the three-dimensional imaging grid to generate three-dimensional focused images of each SAR track in the preset area of ​​the glacier.

[0034] In step S104, the three-dimensional focused images of each SAR track are superimposed with phase coherence within the target imaging layer in the preset region of the glacier to obtain the three-dimensional scattering intensity distribution of the target in the preset region of the glacier.

[0035] The target imaging layer is a layer with a predetermined target depth within a predetermined region of the glacier.

[0036] In step S105, spatial continuity analysis of the sub-ice scatterer is performed based on the target's three-dimensional scattering intensity distribution to obtain imaging results.

[0037] In some embodiments of this application, the method may be executed by a server or by a terminal device with certain processing capabilities.

[0038] In some embodiments of this application, multi-track SAR echo data acquired by tomographic SAR can be processed for imaging, and the imaging results can be calibrated to obtain calibrated complex radar images. Furthermore, a three-dimensional imaging grid for a predetermined glacier region can be constructed based on a digital elevation model of the ice surface.

[0039] Finally, spatial continuity analysis of the subglacial scatterer can be performed based on the target's three-dimensional scattering intensity distribution to obtain imaging results. These imaging results may include, for example, the internal structure of the glacier or information about subglacial interfaces.

[0040] In some embodiments of this application, three-dimensional tomographic focusing can be performed using calibrated complex radar images within a three-dimensional imaging grid to generate three-dimensional focused images of each SAR orbit in a predetermined glacier region. Furthermore, phase coherence superposition of the three-dimensional focused images of each SAR orbit within the target imaging layer of the predetermined glacier region yields the three-dimensional scattering intensity distribution of the target in the predetermined glacier region.

[0041] According to the technical solution provided in the embodiments of this application, the positioning accuracy of subglacial scatterers in three-dimensional space is improved by introducing refraction path correction, avoiding systematic errors caused by the traditional straight-line propagation assumption; the robustness of tomographic imaging under low signal-to-noise ratio conditions is improved by coherence-weighted tomographic focusing, which is conducive to the stable extraction of subglacial continuous interfaces; the method does not depend on specific orbital configurations or specific implementation algorithms, has good versatility and scalability, and is suitable for various glacier radar detection scenarios.

[0042] In some embodiments of this application, the calibrated complex radar images include complex radar images calibrated for each track SAR.

[0043] The complex radar image after SAR calibration for each track is determined in the following way:

[0044] Analyze the SAR echo data of this orbit to obtain the SAR image of this orbit;

[0045] A temporal back-projection algorithm is used to focus the local SAR image onto the ice surface to obtain a two-dimensional imaging result.

[0046] Phase error correction is performed on the two-dimensional SAR imaging results of this orbit based on the phase center dual positioning algorithm to obtain the complex radar image after calibration of this orbit SAR.

[0047] In other words, the raw SAR echo data can be analyzed first, and then the Time Domain Backprojection (TDBP) algorithm can be used to focus the SAR image onto the ice surface for two-dimensional imaging. Finally, based on the Phase Center Double Localization (PCDL) algorithm, the phase error caused by inaccurate orbit information is corrected to obtain the calibrated complex radar image for each orbit.

[0048] In some embodiments of this application, the three-dimensional imaging grid of the predetermined glacier region is discretized into multiple imaging layers along the vertical direction. Each imaging layer corresponds to a two-dimensional plane used to represent the distribution of potential scatterers at different depths. The vertical spacing between adjacent imaging layers is less than the minimum resolution of tomographic SAR.

[0049] In other words, the three-dimensional imaging grid of the glacier's preset area, constructed based on the digital elevation model of the ice surface, is discretized into multiple imaging layers along the depth direction. Each imaging layer corresponds to a two-dimensional plane, which is used to represent the distribution of potential scatterers at different depths.

[0050] In some implementations, the vertical spacing between adjacent imaging layers It should be less than the minimum resolution of tomographic SAR, which can be the Rayleigh resolution achievable by tomographic SAR. ,in Indicates the radar operating wavelength. Indicates radar slant range, This represents the maximum spatial baseline length between different orbits. In one example, the interval between adjacent imaging layers can be set to 10 meters. The maximum depth of the 3D grid should be greater than the maximum depth of the glacier.

[0051] In some embodiments of this application, performing three-dimensional tomographic focusing using calibrated complex radar images within a three-dimensional imaging grid may include:

[0052] The complex radar images after SAR calibration for each track are defocused to obtain the corrected range compressed image. The defocusing process includes remapping the complex radar images after SAR calibration for each track back to the range compressed domain according to the inverse process of time-domain back-projection imaging.

[0053] For the target imaging unit in the three-dimensional imaging grid, the position of the corresponding radar signal ice surface incident point is calculated based on the position of the radar antenna on each track and in each direction at any time; where the target imaging unit is any spatial voxel unit in the three-dimensional imaging grid;

[0054] Based on the location of the incident point on the ice surface, the refraction path of the radar signal in the air and ice medium is corrected to obtain the refractive index-weighted distance of the radar signal from the antenna to the target imaging unit.

[0055] Based on the refractive index-weighted distance from the radar signal from the antenna to each imaging unit, the corrected range-compressed image is refocused onto the three-dimensional imaging grid to obtain the three-dimensional focused image of each SAR track. Here, one imaging unit corresponds to one spatial voxel unit in the three-dimensional imaging grid.

[0056] In some implementations, calculating the location of the radar signal incident point on the ice surface based on the position of the radar antenna at each time along each track and in each direction may include:

[0057] The initial incident point is determined by the intersection of the line connecting the target antenna and the target imaging unit with the ice surface; the target antenna is the radar antenna corresponding to any orbit and any azimuth at any time; the connecting line is the spatial straight line between the phase center of the target antenna and the center point of the target imaging unit.

[0058] The initial incident angle is determined by setting the angle between the target antenna and the target imaging unit and the vertical direction.

[0059] Iteratively update the incident point and incident angle starting from the initial incident point and initial incident angle;

[0060] In response to the determination that the change in the updated incident point position obtained from two adjacent iterations is less than a preset threshold, the iteration stops, and the updated incident point position at the time of stopping the iteration is determined to be the radar signal ice surface incident point position corresponding to the target antenna.

[0061] Iterative updates of the incident point and incident angle, starting from the initial incident point and initial incident angle, can include:

[0062] In the In the next iteration, by solving the equation Determine the updated angle of incidence ; For the refractive index of ice medium, , The line connecting the target antenna and the target imaging unit at the incident point The projected length on the tangent plane of the ice surface. For the target antenna to the incident point The perpendicular distance from the plane tangent to the ice surface; From the target imaging unit to the incident point The perpendicular distance from the plane tangent to the ice surface;

[0063] Using formula The updated spatial incident point is determined based on the updated incident angle; where... Point of incidence The unit vector perpendicular to the azimuth direction of the target antenna on the tangent plane of the ice surface;

[0064] Will Projecting the image vertically onto the ice surface corresponding to the digital elevation model of the ice surface, we obtain the first... The incident point after the next iteration update .

[0065] In some other implementations, the three-dimensional focused images of each SAR track are determined in the following manner:

[0066] The radar signal propagation optical path, considering refraction, is constructed based on the location of the radar signal incident point on the ice surface corresponding to the target antenna. ;in, For distance units, The position of the target antenna in the three-dimensional imaging grid space. The location of the radar signal incident on the ice surface corresponding to the target antenna. The location of the target imaging unit. The refractive index of the ice medium is used; the target antenna is the radar antenna corresponding to any orbit and any azimuth at any given time.

[0067] Construct an azimuth-direction matched filter function corresponding to the optical path length of radar signal propagation;

[0068] Based on the corrected range-compressed image, a three-dimensional time-domain back-projection focusing process is performed using an azimuth-matched filter function to obtain three-dimensional focused images of each SAR track in the preset glacier area; this three-dimensional focused image is different from the two-dimensional complex image, as it is a complex image generated oriented towards the three-dimensional imaging grid.

[0069] Among them, the target imaging unit voxel points The complex value in a three-dimensional focused image is ; For azimuth-time index, For target point The azimuth time corresponding to the direction pointed by the center of the radar antenna beam. This is the synthetic aperture length, measured in units of the number of azimuth sampling points. This is the weighting function for the azimuth antenna pattern. voxel point Corresponding index Distance unit, For index With distance unit The corresponding corrected distance compressed image, For radar operating wavelength, It is an exponential function. It is the symbol for imaginary numbers.

[0070] In other words, when performing three-dimensional tomographic focusing, the complex radar images after SAR calibration of each track can first be defocused to generate a corrected range compressed image.

[0071] In this context, defocusing refers to remapping the phase-calibrated complex radar image back to the range compression domain using the inverse process of back-projection imaging. This process involves reverse-expanding the coherent accumulation of the focused pixels along the azimuth direction, redistributing the complex values ​​of each pixel to the corresponding range cells, thereby obtaining range-compressed radar data consistent with the actual flight trajectory.

[0072] The mathematical expression for this process can be ;in, For distance-compressed images, This indicates the location index in terms of time. Represents a distance unit; This represents the complex value of the focused pixel after phase calibration; The azimuth antenna pattern weighting function; The azimuth time index is The position coordinates of the radar antenna in space. The position coordinates of the target imaging unit, denoted by [symbol]. Represents the magnitude of a vector. For radar operating wavelength, It is an exponential function. It is the symbol for imaginary numbers.

[0073] Next, for each imaging unit in the 3D mesh, based on the radar antenna's position on each orbit and in each direction at any given time, the corresponding radar signal ice surface incident point position can be calculated. Based on the ice surface incident point coordinates, the radar signal propagation path within the air and ice media is corrected for refraction, and the refractive index-weighted distance from the antenna to that imaging unit is calculated. Let... This indicates the coordinates of the radar antenna's position in space. Indicates the position coordinates of the target imaging unit. This indicates the location of the radar wave's incident point on the ice surface. The optical path of the radar signal between the antenna and the target element should be minimized. .in, For the refractive index of ice medium, Let be the set of all points on the ice surface, where The index function is used to find the minimum value.

[0074] In some implementations, the ice surface incident point can be determined iteratively. The intersection of the line connecting the antenna and the imaging unit with the ice surface can be used as the initial incident point. The angle between the connecting line and the vertical direction is recorded as the initial angle of incidence. In the first In this iteration, the update of the incident angle can be obtained by solving the following quartic equation. .

[0075] Seek Then, the angle of incidence It can be represented as The update of the incident point position can be calculated using the following formula. .

[0076] Will Projecting it onto the ice surface will yield the result. The incident point after +1 iterations The iteration is considered to have converged when the change in the incident point position between two consecutive iterations is less than a preset threshold. This preset threshold should be less than the horizontal resolution of the radar image; in one example, the preset threshold can be set to 0.1 meters.

[0077] Considering the optical path length of the refracted radar signal, it can be written as... .

[0078] Finally, the corrected distance-compressed image can be refocused onto the three-dimensional imaging grid based on the calculated corrected refraction path to generate a three-dimensional image of the glacier.

[0079] Specifically, for each imaging unit in the 3D imaging grid, the propagation optical path can be calculated using the method described above. A corresponding azimuth-matched filter function is then constructed based on this propagation optical path, and subsequently, a 3D time-domain back-projection focusing process is performed using this azimuth-matched filter function to obtain the 3D focused images of each SAR track in the preset glacier region. This 3D time-domain back-projection focusing process can be represented as follows: .

[0080] In some embodiments of this application, the three-dimensional scattering intensity distribution of the target in the predetermined region of the glacier can be determined in the following manner:

[0081] Within the target imaging layer, the three-dimensional focused images of each SAR track are paired up to form an interferometric image;

[0082] The average coherence tomographic image of the target imaging layer is obtained by weighting and superimposing the interferometric images based on phase coherence.

[0083] The three-dimensional scattering intensity distribution of the target in the preset region of the glacier is determined based on the average coherence tomography amplitude values ​​of each imaging layer.

[0084] Spatial continuity analysis of subglacial scatterers based on the target's three-dimensional scattering intensity distribution can include:

[0085] For each horizontal pixel coordinate, determine the vertical coordinate with the highest scattering intensity to generate a preliminary estimated image of the glacier depth distribution;

[0086] For each horizontal pixel position in the preliminary estimated image, the position of the imaging unit with the maximum average coherence tomography amplitude value is determined in the vertical direction, and the vertical coordinates corresponding to the imaging unit are used as the scattering depth of the horizontal pixel position.

[0087] Based on the assumption that the ice thickness varies continuously and smoothly in space, the spatial gradient between adjacent pixels is calculated by combining the scattering depths corresponding to all horizontal pixel positions.

[0088] Remove radar pixels whose gradient values ​​are greater than a preset gradient threshold;

[0089] Spatial interpolation is performed on the scattering depth of the remaining radar pixels to generate the ice layer thickness distribution results;

[0090] The imaging results are determined based on the ice thickness distribution.

[0091] In other words, after completing 3D refocusing, to enhance the detectability of deep weak scatterers and reduce the effects of amplitude imbalance and spatial decorrelation, the 3D focusing results can be weighted and superimposed based on interferometric phase coherence. Specifically, this includes:

[0092] First, within the same depth layer, the 3D focusing results acquired from different flight trajectories are paired to form an interferometric image. For example, for a fixed depth... , take the first Track and First The complex coherence of the orbit's three-dimensional focusing result is calculated, and it is defined as follows: ;in, For the first Track and First Complex coherence of orbit 3D focusing results and They represent the first Track and the first Orbit at voxel position Complex scattering value at that point, express The complex conjugate, This represents the multi-view averaging operator in the horizontal direction.

[0093] Then, the interferometric images can be weighted and superimposed based on phase coherence to obtain an average coherence three-dimensional tomographic image. For a given depth layer, all interferometric pairs that satisfy the coherence threshold condition are superimposed, and the expression is: ;in, For the complex coherence after superposition, Indicates depth The set of interference pairs that satisfy the average coherence threshold. ; This represents the set of horizontal coordinates of all points located inside the glacier.

[0094] Furthermore, after obtaining the coherence-weighted three-dimensional tomography results, spatial continuity analysis is performed on the subglacial scatterers to extract information about the glacier's internal structure or subglacial interfaces. Specifically, this includes:

[0095] For each horizontal pixel location, the location of the imaging unit with the maximum average coherence in the vertical direction is determined, and its expression is: ;in, Indicates the height of the ice surface. This represents the lower bound corresponding to the maximum search depth.

[0096] Then, based on the assumption that the ice thickness varies continuously and smoothly in space, the spatial gradient of the scattering depth is calculated, and radar pixels with gradient values ​​greater than a preset threshold are removed to eliminate the influence of noise and clutter. This gradient threshold should be greater than half of the Rayleigh resolution. In one example, 10 meters can be selected as the gradient threshold.

[0097] Finally, spatial interpolation is performed on the scattering depth of the remaining pixels after removing outlier pixels to generate the ice thickness distribution result. The glacier ice reserve is then estimated based on the ice thickness distribution result. In one example, the Kriging interpolation method can be used to calculate the ice thickness, or other interpolation methods can be used; no restrictions are placed here.

[0098] When using the Kriging interpolation method to calculate ice thickness, the radar pixels can first be divided into several spatially connected components, and the connected component with the largest area can be selected as a reliable sample set. Then, ordinary Kriging interpolation is performed on this set to obtain a preliminary estimate of the ice thickness. If necessary, the excluded areas can be compensated by using the block Kriging interpolation method to generate a spatially continuous ice thickness distribution result.

[0099] Figure 2 This is a schematic flowchart of another airborne glacier tomography processing method based on refraction correction and coherent stacking provided in an embodiment of this application. Figure 2 As shown, the method can be divided into three main processes: 3D image focusing, phase calibration, and ice thickness and ice volume estimation.

[0100] When performing a 3D image focusing process, an uncalibrated range-compressed single-view image can first be obtained based on SAR echo data. Then, the TDBP algorithm is used to focus the range-compressed single-view image to obtain a single-view image focused on the ice surface. This single-view image focused on the ice surface can be used as the phase calibration processing object.

[0101] Meanwhile, the three-dimensional image focusing process also includes receiving the calibrated range compressed image output from the phase calibration process, and using the propagation path calculation method to refocus the calibrated range compressed image to obtain the refocused SAR image, and then obtaining the InSAR coherence image.

[0102] When performing the phase calibration process, the single-view image focused on the ice surface output by the 3D image focusing process can be acquired first, and PCDL phase calibration can be performed on it to obtain the calibrated single-view image and the corrected flight trajectory. The calibrated single-view image and the corrected flight trajectory can be defocused to obtain the calibrated distance compressed image.

[0103] When performing the ice thickness and ice volume estimation process, the InSAR coherence image output from the 3D image focusing process can be obtained. The image can then be processed sequentially by image stacking, 3D coherence matrix construction, bedrock pixel extraction, ice thickness distribution map generation, and total ice volume estimation to obtain the estimation results.

[0104] All of the above-mentioned optional technical solutions can be combined in any way to form the optional embodiments of this application, and will not be described in detail here.

[0105] The following are embodiments of the apparatus described in this application, which can be used to execute the embodiments of the method described in this application. For details not disclosed in the apparatus embodiments of this application, please refer to the embodiments of the method described in this application.

[0106] Figure 3 This is a schematic diagram of an airborne glacier tomography processing device based on refraction correction and coherent stacking, provided in an embodiment of this application. Figure 3 As shown, the device includes:

[0107] The pre-imaging module 301 is configured to perform imaging processing on the multi-track SAR echo data acquired by tomographic synthetic aperture radar SAR, and to calibrate the imaging results to obtain a calibrated complex radar image.

[0108] The construction module 302 is configured to construct a three-dimensional imaging grid of a preset glacier region based on the digital elevation model of the ice surface; wherein the preset glacier region includes an area covering a predetermined depth range below the ice surface.

[0109] The focusing module 303 is configured to perform three-dimensional tomographic focusing in a three-dimensional imaging grid using calibrated complex radar images to generate three-dimensional focused images of each orbit of SAR for a preset area of ​​the glacier.

[0110] The overlay module 304 is configured to perform phase coherence overlay on the three-dimensional focused images of each SAR track within the target imaging layer in the preset region of the glacier to obtain the three-dimensional scattering intensity distribution of the target in the preset region of the glacier; wherein, the target imaging layer is a layer in the preset region of the glacier with a preset target depth.

[0111] Imaging module 305 is configured to perform spatial continuity analysis on the sub-ice scatterer based on the target's three-dimensional scattering intensity distribution to obtain imaging results.

[0112] According to the technical solution provided in the embodiments of this application, the positioning accuracy of subglacial scatterers in three-dimensional space is improved by introducing refraction path correction, avoiding systematic errors caused by the traditional straight-line propagation assumption; the robustness of tomographic imaging under low signal-to-noise ratio conditions is improved by coherence-weighted tomographic focusing, which is conducive to the stable extraction of subglacial continuous interfaces; the method does not depend on specific orbital configurations or specific implementation algorithms, has good versatility and scalability, and is suitable for various glacier radar detection scenarios.

[0113] In some implementations, the calibrated complex radar image includes the calibrated complex radar image of each SAR track; wherein, the calibrated complex radar image of each SAR track is determined as follows: the SAR echo data of the track is analyzed to obtain the SAR image of the track; the SAR image of the track is focused onto the ice surface plane using a time-domain back-projection algorithm to obtain a two-dimensional imaging result; the phase error of the two-dimensional imaging result of the track is corrected based on the phase center dual positioning algorithm to obtain the calibrated complex radar image of the track.

[0114] In some implementations, the three-dimensional imaging grid of the glacier's pre-defined region is discretized into multiple imaging layers along the vertical direction. Each imaging layer corresponds to a two-dimensional plane used to represent the distribution of potential scatterers at different depths. The vertical spacing between adjacent imaging layers is less than the minimum resolution of tomographic SAR.

[0115] In some implementations, three-dimensional tomographic focusing is performed on calibrated complex radar images within a three-dimensional imaging grid. This includes: defocusing the calibrated complex radar images of each SAR track to obtain a corrected range-compressed image; the defocusing process includes remapping the calibrated complex radar images of each SAR track back to the range-compressed domain according to the inverse process of time-domain back-projection imaging; calculating the corresponding radar signal ice surface incident point position for the target imaging unit in the three-dimensional imaging grid based on the position of the radar antenna in each track and in each direction; wherein, the target imaging unit is any spatial voxel unit in the three-dimensional imaging grid; correcting the refraction path of the radar signal propagation path within the air and ice medium based on the ice surface incident point position to obtain the refractive index-weighted distance of the radar signal from the antenna to the target imaging unit; and refocusing the corrected range-compressed image onto the three-dimensional imaging grid based on the refractive index-weighted distance of the radar signal from the antenna to each imaging unit to obtain a three-dimensional focused image of each SAR track.

[0116] In some implementations, the location of the radar signal incident point on the ice surface is calculated based on the position of the radar antenna at each time along each track and in each azimuth direction. This includes: determining the intersection of the line connecting the target antenna and the target imaging unit with the ice surface as the initial incident point; the target antenna being the radar antenna corresponding to any time along any track and in any azimuth direction; the connecting line being the spatial straight line between the phase center of the target antenna and the center point of the target imaging unit; determining the angle between the target antenna and the target imaging unit and the vertical direction as the initial incident angle; iteratively updating the incident point and incident angle starting from the initial incident point and initial incident angle; and stopping the iteration in response to determining that the change in the updated incident point location obtained from two adjacent iterations is less than a preset threshold, and determining that the updated incident point location at the time of stopping the iteration is the radar signal incident point location on the ice surface corresponding to the target antenna.

[0117] In some implementations, the incident point and incident angle are iteratively updated starting from the initial incident point and initial incident angle, including: in the first... In the next iteration, by solving the equation Determine the updated angle of incidence ; For the refractive index of ice medium, , The line connecting the target antenna and the target imaging unit at the incident point The projected length on the tangent plane of the ice surface. For the target antenna to the incident point The perpendicular distance from the plane tangent to the ice surface; From the target imaging unit to the incident point The perpendicular distance to the tangent plane of the ice surface; using the formula The updated spatial incident point is determined based on the updated incident angle; where... Point of incidence The unit vector perpendicular to the azimuth direction of the target antenna on the tangent plane of the ice surface; Projecting the image vertically onto the ice surface corresponding to the digital elevation model of the ice surface, we obtain the first... The incident point after the next iteration update .

[0118] In some implementations, the three-dimensional focused images of each SAR track are determined as follows: the radar signal propagation optical path, taking into account refraction, is constructed based on the location of the radar signal incident point on the ice surface corresponding to the target antenna. ;in, For distance units, The position of the target antenna in the three-dimensional imaging grid space. The location of the radar signal incident on the ice surface corresponding to the target antenna. The location of the target imaging unit. The refractive index of the ice medium is given; the target antenna is the radar antenna corresponding to any orbit and any azimuth time; an azimuth matched filter function corresponding to the propagation optical path of the radar signal is constructed; based on the corrected range compressed image, three-dimensional time-domain back-projection focusing processing is performed using the azimuth matched filter function to obtain the three-dimensional focused image of each orbit SAR of the preset glacier area; where the target imaging unit voxel points The complex value in a three-dimensional focused image is ; For azimuth-time index, For target point The azimuth time corresponding to the direction pointed by the center of the radar antenna beam. This is the synthetic aperture length, measured in units of the number of azimuth sampling points. This is the weighting function for the azimuth antenna pattern. voxel point Corresponding index Distance unit, For index With distance unit The corresponding corrected distance compressed image, For radar operating wavelength, It is an exponential function. It is the symbol for imaginary numbers.

[0119] In some implementations, the three-dimensional scattering intensity distribution of the target in the predetermined region of the glacier is determined as follows: within the target imaging layer, the three-dimensional focused images of each SAR orbit are paired to form an interferometric image; the interferometric images are weighted and superimposed according to phase coherence to obtain the average coherence tomography image of the target imaging layer; the three-dimensional scattering intensity distribution of the target in the predetermined region of the glacier is determined based on the amplitude value of the average coherence tomography image of each imaging layer.

[0120] In some implementations, spatial continuity analysis of the subglacial scatterer is performed based on the target's three-dimensional scattering intensity distribution, including: for each horizontal pixel coordinate, determining the vertical coordinate with the highest scattering intensity to generate a preliminary estimated image of the glacier depth distribution; for each horizontal pixel position in the preliminary estimated image, determining the position of the imaging unit with the largest average coherence tomography amplitude value in the vertical direction, and using the vertical coordinate corresponding to the imaging unit as the scattering depth of that horizontal pixel position; based on the assumption that the ice thickness changes continuously and smoothly in space, calculating the spatial gradient between adjacent pixels by combining the scattering depths corresponding to all horizontal pixel positions; removing radar pixels with gradient values ​​greater than a preset gradient threshold; spatially interpolating the scattering depths of the remaining radar pixels to generate the ice thickness distribution result; and determining the imaging result based on the ice thickness distribution result.

[0121] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0122] Figure 4 This is a schematic diagram of the electronic device provided in an embodiment of this application. For example... Figure 4 As shown, the electronic device 4 of this embodiment includes: a processor 401, a memory 402, and a computer program 403 stored in the memory 402 and executable on the processor 401. When the processor 401 executes the computer program 403, it implements the steps in the various method embodiments described above. Alternatively, when the processor 401 executes the computer program 403, it implements the functions of each module / unit in the various device embodiments described above.

[0123] Electronic device 4 can be a desktop computer, laptop, handheld computer, cloud server, or other electronic device. Electronic device 4 may include, but is not limited to, processor 401 and memory 402. Those skilled in the art will understand that... Figure 4 This is merely an example of electronic device 4 and does not constitute a limitation on electronic device 4. It may include more or fewer components than shown, or different components.

[0124] The processor 401 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc.

[0125] The memory 402 can be an internal storage unit of the electronic device 4, such as a hard disk or RAM of the electronic device 4. The memory 402 can also be an external storage device of the electronic device 4, such as a plug-in hard disk, Smart Media Card (SMC), Secure Digital (SD) card, Flash Card, etc., equipped on the electronic device 4. The memory 402 can also include both internal and external storage units of the electronic device 4. The memory 402 is used to store computer programs and other programs and data required by the electronic device.

[0126] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is merely an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiments can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0127] If integrated modules / units are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the methods of the above embodiments can also be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program may include computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. Computer-readable media may include: any entity or device capable of carrying computer program code, recording media, USB flash drives, portable hard drives, magnetic disks, optical disks, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc.

[0128] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. An airborne glacier tomography processing method based on refraction correction and coherent stacking, characterized in that, include: The multi-track SAR echo data acquired by tomographic synthetic aperture radar (SAR) is processed for imaging, and the imaging results are calibrated to obtain a calibrated complex radar image. A three-dimensional imaging grid is constructed based on the digital elevation model of the ice surface for a predetermined glacier region; the predetermined glacier region includes an area covering a predetermined depth range below the ice surface. Three-dimensional tomography focusing is performed using the calibrated complex radar images in the three-dimensional imaging grid to generate three-dimensional focused images of each SAR track in the preset area of ​​the glacier; Phase coherence superposition of the three-dimensional focused images of each SAR track within the target imaging layer in the preset glacier region yields the three-dimensional scattering intensity distribution of the target in the preset glacier region; wherein, the target imaging layer is the layer in the preset glacier region with the preset target depth; Based on the three-dimensional scattering intensity distribution of the target, spatial continuity analysis of the sub-ice scatterer is performed to obtain imaging results; The process of performing three-dimensional tomographic focusing using the calibrated complex radar image within the three-dimensional imaging grid includes: The complex radar images after SAR calibration for each track are defocused to obtain a corrected range compressed image; the defocusing process includes remapping the complex radar images after SAR calibration for each track back to the range compressed domain according to the inverse process of time-domain back-projection imaging. For the target imaging unit in the three-dimensional imaging grid, the position of the corresponding radar signal ice surface incident point is calculated based on the position of the radar antenna on each track and in each direction at the time; wherein, the target imaging unit is any spatial voxel unit in the three-dimensional imaging grid; Based on the location of the incident point on the ice surface, the refraction path of the radar signal in the air and ice medium is corrected to obtain the refractive index-weighted distance of the radar signal from the antenna to the target imaging unit. Based on the refractive index-weighted distance from the radar signal from the antenna to each imaging unit, the corrected range-compressed image is refocused onto the three-dimensional imaging grid to obtain a three-dimensional focused image of each SAR track.

2. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 1, characterized in that, The calibrated complex radar images include complex radar images calibrated for each SAR track; The complex radar image after SAR calibration for each track is determined in the following way: Analyze the SAR echo data of this orbit to obtain the SAR image of this orbit; A temporal back-projection algorithm is used to focus the local SAR image onto the ice surface to obtain a two-dimensional imaging result. Phase error correction is performed on the two-dimensional imaging results of the local SAR based on the phase center dual positioning algorithm to obtain the complex radar image after local SAR calibration.

3. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 1, characterized in that, The three-dimensional imaging grid of the glacier's preset area is discretized into multiple imaging layers along the vertical direction. Each imaging layer corresponds to a two-dimensional plane, which is used to represent the distribution of potential scatterers at different depths. The vertical spacing between adjacent imaging layers is less than the minimum resolution of tomographic SAR.

4. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 1, characterized in that, The location of the radar signal incident point on the ice surface is calculated based on the position of the radar antenna on each track and in each direction at that time, including: The initial incident point is determined by the intersection of the line connecting the target antenna and the target imaging unit with the ice surface; the target antenna is the radar antenna corresponding to any orbit and any azimuth at any time; the connecting line is the spatial straight line between the phase center of the target antenna and the center point of the target imaging unit. The initial incident angle is determined by setting the angle between the target antenna and the target imaging unit and the vertical direction. Iteratively update the incident point and incident angle starting from the initial incident point and initial incident angle; In response to the determination that the change in the updated incident point position obtained from two adjacent iterations is less than a preset threshold, the iteration stops, and the updated incident point position at the time of stopping the iteration is determined to be the radar signal ice surface incident point position corresponding to the target antenna.

5. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 4, characterized in that, Iteratively update the incident point and incident angle starting from the initial incident point and initial incident angle, including: In the In the next iteration, by solving the equation Determine the updated angle of incidence ; For the refractive index of ice medium, , The line connecting the target antenna and the target imaging unit at the incident point The projected length on the tangent plane of the ice surface. For the target antenna to the incident point The perpendicular distance from the plane tangent to the ice surface; From the target imaging unit to the incident point The perpendicular distance from the plane tangent to the ice surface; Using formula The updated spatial incident point is determined based on the updated incident angle; where... Point of incidence The unit vector perpendicular to the azimuth direction of the target antenna on the tangent plane of the ice surface; Will Projecting the image vertically onto the ice surface corresponding to the digital elevation model of the ice surface, we obtain the first... The incident point after the next iteration update .

6. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 1, characterized in that, The three-dimensional focused images of each SAR orbit were determined using the following method: The radar signal propagation optical path, considering refraction, is constructed based on the location of the radar signal incident point on the ice surface corresponding to the target antenna. ;in, For distance units, The position of the target antenna in the three-dimensional imaging grid space. The location of the radar signal incident on the ice surface corresponding to the target antenna. The location of the target imaging unit. The refractive index of the ice medium is used; the target antenna is the radar antenna corresponding to any orbit and any azimuth at any given time. Construct the azimuth matched filter function corresponding to the optical path of the radar signal propagation; Based on the corrected range compressed image, the three-dimensional time domain back projection focusing process is performed using the azimuth matched filtering function to obtain the three-dimensional focused image of each SAR track in the preset area of ​​the glacier. Among them, the target imaging unit voxel points The complex value in a three-dimensional focused image is ; For azimuth-time index, For target point The azimuth time corresponding to the direction pointed by the center of the radar antenna beam. This is the synthetic aperture length, measured in units of the number of azimuth sampling points. This is the weighting function for the azimuth antenna pattern. voxel point Corresponding index Distance unit, For index With distance unit The corresponding corrected distance compressed image, For radar operating wavelength, It is an exponential function. It is the symbol for imaginary numbers.

7. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 1, characterized in that, The three-dimensional scattering intensity distribution of the target in the pre-defined region of the glacier is determined using the following method: Within the target imaging layer, the three-dimensional focused images of each SAR track are paired up to form an interferometric image; The average coherence tomographic image of the target imaging layer is obtained by weighting and superimposing the interferometric images based on phase coherence. The three-dimensional scattering intensity distribution of the target in the preset region of the glacier is determined based on the average coherence tomography amplitude values ​​of each imaging layer.

8. The airborne glacier tomography processing method based on refraction correction and coherence stacking according to claim 7, characterized in that, Spatial continuity analysis of the subglacial scatterer based on the target's three-dimensional scattering intensity distribution includes: For each horizontal pixel coordinate, determine the vertical coordinate with the highest scattering intensity to generate a preliminary estimated image of the glacier depth distribution; For each horizontal pixel position in the preliminary estimated image, the position of the imaging unit with the maximum average coherence tomography amplitude value is determined in the vertical direction, and the vertical coordinates corresponding to the imaging unit are used as the scattering depth of the horizontal pixel position. Based on the assumption that the ice thickness varies continuously and smoothly in space, the spatial gradient between adjacent pixels is calculated by combining the scattering depths corresponding to all horizontal pixel positions. Remove radar pixels whose gradient values ​​are greater than a preset gradient threshold; Spatial interpolation is performed on the scattering depth of the remaining radar pixels to generate the ice layer thickness distribution results; The imaging results are determined based on the ice thickness distribution.

9. An airborne glacier tomography processing device based on refraction correction and coherent stacking, characterized in that, include: The pre-imaging module is configured to perform imaging processing on the multi-track SAR echo data acquired by tomographic synthetic aperture radar SAR, and to calibrate the imaging results to obtain a calibrated complex radar image. The construction module is configured to build a 3D imaging mesh of a preset glacier region based on the digital elevation model of the ice surface; wherein the preset glacier region includes an area covering a predetermined depth range below the ice surface; The focusing module is configured to perform three-dimensional tomographic focusing in the three-dimensional imaging grid using the calibrated complex radar images to generate three-dimensional focused images of each orbit SAR of the preset glacier region; The overlay module is configured to perform phase coherence overlay on the three-dimensional focused images of each SAR track within the target imaging layer in the preset region of the glacier to obtain the three-dimensional scattering intensity distribution of the target in the preset region of the glacier; wherein, the target imaging layer is a layer in the preset region of the glacier with a preset target depth; The imaging module is configured to perform spatial continuity analysis on the sub-ice scatterer based on the three-dimensional scattering intensity distribution of the target, and obtain imaging results; The process of performing three-dimensional tomographic focusing using the calibrated complex radar image within the three-dimensional imaging grid includes: The complex radar images after SAR calibration for each track are defocused to obtain a corrected range compressed image; the defocusing process includes remapping the complex radar images after SAR calibration for each track back to the range compressed domain according to the inverse process of time-domain back-projection imaging. For the target imaging unit in the three-dimensional imaging grid, the position of the corresponding radar signal ice surface incident point is calculated based on the position of the radar antenna on each track and in each direction at the time; wherein, the target imaging unit is any spatial voxel unit in the three-dimensional imaging grid; Based on the location of the incident point on the ice surface, the refraction path of the radar signal in the air and ice medium is corrected to obtain the refractive index-weighted distance of the radar signal from the antenna to the target imaging unit. Based on the refractive index-weighted distance from the radar signal from the antenna to each imaging unit, the corrected range-compressed image is refocused onto the three-dimensional imaging grid to obtain a three-dimensional focused image of each SAR track.

Images