Long-time-sequence glacier flow velocity monitoring method

Abstract

The invention relates to the technical field of image data processing, and discloses a long-time-sequence glacier flow velocity monitoring method, which comprises the following steps of: acquiring a glacier SAR (Synthetic Aperture Radar) image sequence and preprocessing to generate a plurality of multi-view SAR images; setting a time baseline threshold value, performing image pair combination on all the multi-view SAR images, performing offset estimation, and respectively generating a plurality of offset sequences in azimuth and distance directions; constructing a nonlinear glacier flow velocity model based on seasonal components, optimizing all offset sequences, generating an optimized offset sequence, calculating a glacier surface average flow velocity field, and generating a glacier flow velocity direction matrix in combination with second glacier cataloguing data; taking the optimized offset sequence as an observed quantity, setting a glacier flow direction as a constraint condition, and carrying out iterative solution on displacement distribution corresponding to the imaging moment of each multi-view SAR image by adopting a least square method to obtain optimal displacement distribution corresponding to the imaging moment of each multi-view SAR image; according to the method, the reliability and precision of glacier flow velocity calculation are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of image data processing, and particularly relates to a long-time-series glacier flow velocity monitoring method. Background Art

[0002] With the application of spaceborne earth observation technology in glacier resource surveys, glacier remote sensing has entered a rapid development stage and has currently become an important support for glacier inventory, flow velocity monitoring, and ablation estimation. Compared with traditional ground observation technologies, glacier remote sensing monitoring has advantages such as wide coverage, high spatio-temporal resolution, and low cost. As a supplement and extension of optical remote sensing for earth observation, Synthetic Aperture Radar (SAR) provides an opportunity for high-time-frequency flow velocity dynamic monitoring in the southeastern hinterland of the Qinghai-Tibet Plateau where marine glaciers are intensively developed due to its advantage of being unaffected by clouds, fog, rain, snow, and light and shadow. Taking advantage of the SAR remote sensing data acquisition, it can provide data guarantee for revealing the evolution law of the glacier change process under the background of climate warming on a long time scale.

[0003] As one of the regions where modern glaciers are intensively developed in the Hengduan Mountains, the glaciers in the Gongga Mountain area show an obvious trend of retreat, especially in the past ten years. As a typical monsoon marine glacier development area, the investigation of the current situation of the Gongga Mountain glaciers is crucial for comprehensively monitoring the response of glaciers under the background of climate change on the Qinghai-Tibet Plateau; at the same time, the evaluation of its dynamic evolution has important reference and significance for revealing the mechanism of glacier retreat in areas with similar hydrothermal conditions in southeastern Tibet. In addition, as a tourist resource, the Hailuogou Scenic Area on the east slope of Gongga Mountain attracts domestic and foreign tourists to stop and visit due to its beautiful glacier landscape. Based on this, carrying out the monitoring and analysis of the dynamic evolution of marine glaciers in the Gongga Mountain area is extremely crucial for systematically investigating the periglacial disaster-prone environment and dynamically evaluating the glacier disaster risk.

[0004] Existing research has confirmed that glacier movement is significantly correlated with meteorological factors such as temperature and rainfall, manifested as the periodic fluctuation of ice flow velocity in the cold season and warm season with a natural year as the cycle. Therefore, when reconstructing the glacier flow velocity field for a multi-year long time series, existing methods such as linear or high-order polynomial models of the observation equation are difficult to adapt to the significant periodic flow characteristics of glaciers. In addition, the assumption that the repeated observation values of the flow velocity in the glacier area follow or approximately follow a normal distribution in the long time series no longer holds. In short, when modeling the glacier flow velocity for a long time series, whether it is the refinement process of the offset stack sequence or the observation equation of glacier flow velocity modeling, if the periodic activity law facing seasonal meteorological conditions is not considered, it will inevitably reduce the reliability and accuracy of glacier flow velocity calculation. Summary of the Invention

[0005] In view of the above deficiencies in the prior art, the present invention provides a long-term glacier flow velocity monitoring method to solve the problems of poor reliability and low accuracy in glacier flow velocity calculation existing in the existing glacier flow velocity monitoring methods.

[0006] In order to achieve the above invention object, the technical solution adopted by the present invention is as follows: A long-term glacier flow velocity monitoring method includes the following steps: S1. Obtain a glacier SAR image sequence and perform preprocessing to generate a number of multi-looked SAR images; S2. Set a time baseline threshold, combine all multi-looked SAR images into image pairs, and then perform offset estimation to generate a number of offset sequences in the azimuth direction and the range direction respectively; S3. Construct a non-linear glacier flow velocity model based on seasonal components, optimize all offset sequences, and generate optimized offset sequences; S4. Based on the optimized offset sequences, calculate the average flow velocity field on the glacier surface, and combine with the second glacier inventory data to generate a glacier flow velocity direction matrix; S5. Use the optimized offset sequences as observables, set the glacier flow direction as a constraint condition, and use the least squares method to iteratively solve the displacement distribution corresponding to the imaging time of each multi-looked SAR image to obtain the optimal displacement distribution corresponding to the imaging time of each multi-looked SAR image.

[0007] The present invention has the following beneficial effects: The long-term glacier flow velocity monitoring method proposed by the present invention optimizes the offset sequences by constructing a non-linear glacier basin model based on seasonal components, improves the accuracy of outlier filtering in the offset sequences, and solves the problem that it is difficult for existing models such as linear or high-order polynomial models of glacier surface flow velocity observation equations to adapt to the significant periodic flow characteristics of glaciers during the reconstruction of the glacier flow velocity field over long time series across years. Thus, the reliability and accuracy of long-term glacier flow velocity calculation are improved. Description of the Drawings

[0008] Figure 1 It is a flow schematic diagram of a long-term glacier flow velocity monitoring method proposed by the present invention. Detailed Embodiments

[0009] The following describes the detailed embodiments of the present invention to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the detailed embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0010] As Figure 1As shown in the figure, a long-time-series glacier flow velocity monitoring method includes the following steps S1 - S5: S1. Obtain a glacier SAR image sequence and perform preprocessing to generate a number of multi-looked SAR images.

[0011] Specifically, step S1 specifically includes S11 - S12: S11. Obtain a glacier SAR image sequence, which includes a number of glacier SAR images.

[0012] S12. After mosaicking and sub-pixel registration of all glacier SAR images, perform multi-looking operations on the azimuth and range directions of the registered glacier SAR images according to the proportionality coefficient to generate a number of multi-looked SAR images.

[0013] In this embodiment, by registering and multi-looking all glacier SAR images, the offset estimation error caused by registration error and noise is reduced.

[0014] S2. Set a time baseline threshold, combine all multi-looked SAR images into image pairs, and then perform offset estimation to generate a number of offset sequences in the azimuth and range directions respectively.

[0015] Specifically, step S2 specifically includes S21 - S22: S21. Set a time baseline threshold and combine all multi-looked SAR images into image pairs to generate a number of image pairs.

[0016] In this embodiment, considering the difference in backscattering intensity between the warm and cold seasons on the glacier surface, the Hailuogou Glacier on the east slope of Gongga Mountain in the Hengduan Mountains on the southeastern edge of the Qinghai-Tibet Plateau is selected as a typical experimental site, and 90 archived Sentinel-1 satellite ascending-track SAR images (glacier SAR images) between January 2019 and January 2022 are selected as experimental data. Referring to the ablation period law of the Hailuogou Glacier from April to November every year, the time periods from April to June, May to November, and October to May of the following year are defined as independent intervals for this dataset respectively, and a time baseline threshold of 150 days is given within each imaging interval to freely combine image pairs. Therefore, through the selection of the time baseline of the image pairs in this step, the subsequent glacier flow velocity calculation is made more robust.

[0017] S22. Use the optical flow motion estimation method to calculate the offset of each image pair to generate offset sequences in the azimuth and range directions respectively, and finally obtain a number of offset sequences in the azimuth and range directions.

[0018] In this embodiment, the optical flow motion estimation method is used to calculate the offset of each image pair. Each image pair will respectively generate a set of offset results in the azimuth and range directions to generate offset sequences in the azimuth and range directions respectively, and finally generate a number of offset sequences in the azimuth and range directions.

[0019] S3. Construct a non - linear glacier flow velocity model based on seasonal components, optimize all offset sequences, and generate optimized offset sequences.

[0020] Specifically, step S3 includes S31 - S36: S31. Construct a non - linear glacier flow velocity model based on seasonal components, which includes a periodic function and an aperiodic function; among them, the periodic function is established by the periodic main control factors including temperature and rainfall through sine and cosine functions; the aperiodic function is established by the aperiodic main control factors including terrain, mass, and stress through a second - order polynomial function.

[0021] Specifically, the formula of the non - linear glacier flow velocity model based on seasonal components in step S31 is:

[0022] Wherein, represents the glacier flow velocity at the imaging time , represents the sine function, represents the cosine function, , , , all represent the coefficients of the periodic function, , , all represent the coefficients of the aperiodic function.

[0023] In this embodiment, by assuming the glacier flow velocity and decomposing based on the periodic main control factors and aperiodic main control factors, a non - linear glacier flow velocity model based on seasonal components composed of a periodic function and an aperiodic function is constructed. The periodic main control factors composed of temperature and rainfall in this model can be represented by the periodic function composed of the sine function and the cosine function , that is . Then, the aperiodic main control factors composed of terrain, mass, stress, etc. in this model can be represented by a second - order polynomial function, that is .

[0024] S32. Based on all offset sequences, select an observation sample subset by setting a time threshold, specifically: S321. Set the upper limit and the lower limit of the time threshold.

[0025] S322: Based on all offset sequences, select image pairs that satisfy the relationship to generate an observation sample subset; among them, represents the Time interval of the group of image pairs.

[0026] In this embodiment, the time interval of the group of image pairs is obtained by subtracting the imaging times of the two multi-look SAR images in the image pair.

[0027] S33. Based on the subset of observation samples, the coefficients of the periodic function and the non-periodic function are fitted pixel by pixel using the non-linear least squares method.

[0028] S34. Based on the coefficients of the fitted periodic function and non-periodic function, calculate the velocity confidence interval of the offset sequence in the subset of observation samples to filter the subset of observation samples.

[0029] Specifically, the specific process of filtering the subset of observation samples in step S34 is as follows: First, define the set of median times of all image pairs in the subset of observation samples as: , , , , , , respectively represent the imaging times of the 1st, 2nd, 3rd, 4th, , th multi-look SAR images, represents the median time of the image pair composed of the 1st and 2nd multi-look SAR images, represents the median time of the image pair composed of the 1st and 3rd multi-look SAR images, represents the median time of the image pair composed of the 1st and 4th multi-look SAR images, represents the median time of the image pair composed of the 2nd and 3rd multi-look SAR images, represents the th and the th multi-look SAR images.

[0030] Define the set of all image pairs as: , , , , , , respectively represent the th and the Image pairs composed of multiple-look SAR images, image pairs composed of the first and the second multiple-look SAR images, image pairs composed of the first and the third multiple-look SAR images, image pairs composed of the first and the fourth multiple-look SAR images, image pairs composed of the second and the third multiple-look SAR images, the image pairs composed of the th and the

[0031] multiple-look SAR images. , that is:

[0032] wherein, represents the instantaneous at the median time of the image pair composed of the th and the multiple-look SAR images, represents the instantaneous velocity at the median time of the image pair composed of the first and the second multiple-look SAR images, represents the instantaneous velocity at the median time of the image pair composed of the first and the third multiple-look SAR images, represents the instantaneous velocity at the median time of the image pair composed of the second and the third multiple-look SAR images, represents the instantaneous at the median time of the image pair composed of the

[0033] multiple-look SAR images. Then, substitute the median time of each image pair and its corresponding instantaneous velocity into the non-linear glacier flow velocity model based on seasonal components. After fitting the model coefficients using non-linear least squares method, calculate the upper limit and the lower limit

[0034] of the velocity at the median time of the observed sample subset within the 95% confidence interval. In this embodiment, the periodic function and non-periodic function coefficients of the model are estimated using the least squares method. The relationship between the median time of each image pair and its corresponding instantaneous velocity is quantified with these coefficients, and the uncertainty range of the estimated model coefficients is measured through the 95% confidence interval, that is, the upper limit and the lower limit

[0035] of the instantaneous velocity at the median time, providing a basis for filtering out outliers in the subsequent process. Finally, the upper limit Outlier rejection is performed to generate a filtered subset of the observed samples.

[0036] S35. Determine whether there are still outliers in the filtered subset of the observed samples. If so, execute step S34; otherwise, obtain the initial optimized coefficients of the periodic function and the non-periodic function of the non-linear glacier flow velocity model based on the seasonal component.

[0037] In this embodiment, through the cyclic refinement process, until no outliers are filtered from the subset of the observed samples, that is, when the number of the subset of the observed samples remains constant, the initial optimized coefficients of the model are obtained.

[0038] S36. Based on the initial optimized coefficients, calculate the flow velocity confidence intervals of all offset sequences to filter the offset sequences, and determine whether all pixels of each offset sequence have been optimized. If so, generate an optimized offset sequence and execute step S4; otherwise, execute step S32.

[0039] In summary, this step optimizes the offset sequence by constructing a non-linear glacier flow velocity model based on the seasonal component, thereby improving the reliability of the glacier flow velocity.

[0040] S4. Based on the optimized offset sequence, calculate the average flow velocity field on the glacier surface, and combine it with the second glacier inventory data to generate a glacier flow velocity direction matrix.

[0041] Specifically, step S4 specifically includes S41 - S42: S41. According to the optimized offset sequence, calculate the average flow velocity field on the glacier surface, that is:

[0042] Among them, represents the average flow velocity field on the glacier surface, represents the total number of optimized offset sequences, represents the th optimized offset sequence, represents the th time interval of the

[0043] In this embodiment, for the optimized offset sequence , its corresponding time interval is , so the calculation expression for obtaining the average flow velocity field on the glacier surface is . Among them, , , respectively represent the 1st, 2nd, and th optimized offset sequences; , , respectively represent the 1st, 2nd, and The time interval of an optimized offset sequence.

[0044] S42. Based on the second glacier inventory data, the glacier area and non-glacier area are obtained. Based on the average flow velocity field on the glacier surface, the glacier area is divided into a negative area and a positive area, and finally a glacier flow velocity direction matrix is generated.

[0045] In this embodiment, the glacier area and non-glacier area are divided based on the second glacier inventory data, and for coordinate pixels, if the flow velocity direction is in the positive area, it is marked as 1, otherwise it is marked as -1. If it is in the non-glacier area, it is marked as 0, so as to generate a glacier flow direction matrix as one of the constraint conditions for iterative solution in subsequent steps.

[0046] S5. Taking the optimized offset sequence as the observable quantity, setting the glacier flow direction as the constraint condition, and using the least squares method to iteratively solve the displacement distribution corresponding to the imaging time of each multi-look SAR image, so as to obtain the optimal displacement distribution corresponding to the imaging time of each multi-look SAR image.

[0047] Specifically, step S5 specifically includes S501 - S510: S501. Taking the optimized offset sequence as the observable quantity, that is:

[0048] Among them, represents the observable quantity, , , respectively represent the 1st, 2nd, th optimized offset sequences.

[0049] S502. Taking the displacement distribution corresponding to the imaging time of each multi-look SAR image as the quantity to be solved, that is:

[0050] Among them, represents the quantity to be solved, , , respectively represent the displacement distributions corresponding to the imaging times of the 1st, 2nd, th multi-look SAR images.

[0051] S503. According to the offset sequences of each image pair, construct an observation equation for the displacement on the glacier surface, that is:

[0052] Among them, represents order optimized offset sequence matrix, represents The coefficient matrix of the order represents the matrix of quantities to be solved

[0053] Among them, the offset sequence of each image pair is as follows: , represents the th and the th multi-view SAR image pair's offset sequence, represents the displacement distribution corresponding to the imaging time of the th multi-view SAR image, represents the displacement distribution corresponding to the imaging time of the th multi-view SAR image. In addition, the imaging time of the th multi-view SAR image is , the imaging time of the th multi-view SAR image is , and .

[0054] In this embodiment, the purposes of steps S501 - S503 are: by taking the displacement distributions at the imaging times of all multi-view SAR images as the quantities to be calculated, and the optimized offset sequence as the observed quantity, an observation equation is constructed from observed quantities and unknown quantities to be solved, so as to restore the displacement distributions at the imaging times of all multi-view SAR images. Among them, since each row vector of the matrix corresponds to an image pair, the th and the th column coefficients of each row of the coefficient matrix are 1 and -1 respectively, and the rest are all 0, that is: .

[0055] S504. Perform singular value decomposition on the coefficient matrix of the observation equation of the glacier surface displacement to generate a coefficient matrix equation of the singular value decomposition, that is:

[0056] Among them, represents the coefficient matrix of the singular value decomposition, represents order unitary matrix, represents order non-negative diagonal matrix, represents order unitary matrix, represents the transpose.

[0057] In this embodiment, due to the limitations of the time baseline and spatial baseline thresholds, First-order optimized offset sequence matrix They are often discontinuous, manifested as the fracture of a certain image pair in the time domain, causing the coefficient matrix to be singular, and then leading to an underdetermined system of equations, resulting in infinitely many sets of solutions. Therefore, it is necessary to perform singular value decomposition on the coefficient matrix to estimate the approximate solution in the sense of the minimum norm.

[0058] S505. Substitute the coefficient matrix equation obtained by singular value decomposition into the observation equation of the glacier surface displacement, and at the same time take the glacier flow direction as the first constraint condition to generate an observation equation that constrains the flow velocity direction, that is:

[0059] where, represents taking the minimum value, represents the optimized offset sequence matrix, represents the optimized offset sequence, represents the matrix of quantities to be settled, represents the glacier flow velocity direction matrix, represents the th optimized offset sequence.

[0060] S506. Use the least squares method to perform the first-round solution on the observation equation that constrains the flow velocity direction to obtain the initial solution of the displacement distribution corresponding to the imaging time of each multi-look SAR image.

[0061] S507. Convert the initial solutions of the displacement distributions corresponding to the imaging times of each multi-look SAR image to generate the initial flow velocities of the displacement distributions corresponding to the imaging times of each multi-look SAR image, that is:

[0062] where, represents the set of initial flow velocities of the displacement distributions corresponding to the imaging times of each multi-look SAR image, represents the initial flow velocity of the displacement distribution corresponding to the imaging time of the first multi-look SAR image, represents the initial flow velocity of the displacement distribution corresponding to the imaging time of the second multi-look SAR image, represents the th initial flow velocity of the displacement distribution corresponding to the imaging time of the multi-look SAR image.

[0063] In this embodiment, it is assumed that the ground surface is linearly deformed within the imaging times of two adjacent multi-look SAR images. The initial solutions of the displacement distributions corresponding to the imaging times of each multi-look SAR image can be converted into velocities to estimate the model coefficients again, so as to construct the second constraint condition for solving the observation equation of the glacier surface displacement, that is, the glacier displacement gradient constraint condition.

[0064] S508. Substitute the initial flow velocity of the displacement distribution corresponding to the imaging time of each multi-view SAR image and the imaging time into the non-linear glacier flow velocity model based on seasonal components for fitting estimation to generate the optimal coefficients of the periodic function and the non-periodic function.

[0065] S509. Based on the optimal coefficients, use the non-linear glacier flow velocity model based on seasonal components again to calculate the upper limit of the 95% confidence interval corresponding to the initial flow velocity corresponding to the imaging time of each multi-view SAR image and the lower limit , and generate the second constraint condition.

[0066] Specifically, the second constraint condition in step S510 is:

[0067] where , respectively represent the lower limits of the 95% confidence intervals corresponding to the initial flow velocities corresponding to the imaging times of the -th and -th multi-view SAR images, , respectively represent the upper limits of the 95% confidence intervals corresponding to the initial flow velocities corresponding to the imaging times of the -th and -th multi-view SAR images, represents the imaging time of the -th multi-view SAR image, represents the imaging time of the -th multi-view SAR image.

[0068] In this embodiment, the relationship between the imaging time of each image and its corresponding initial flow velocity is quantified based on the optimal coefficients, and the uncertainty range of the estimated model coefficients is measured by the 95% confidence interval, further ensuring the robustness of the non-linear glacier flow velocity model based on seasonal components.

[0069] S510. Based on the second constraint condition, use the least squares method to perform a second round of solution on the observation equation for constraining the flow velocity direction to obtain the optimal displacement distribution corresponding to the imaging time of each multi-view SAR image.

[0070] In this embodiment, the initial solution of the displacement distribution corresponding to the imaging time of each multi-view SAR image obtained in the first round is converted into an initial flow velocity, and then substituted into the non-linear glacier flow velocity model based on seasonal components together with the imaging time for calculation to calculate the upper and lower limits of the velocity corresponding to the imaging time, and use this as the displacement gradient constraint to perform a second round of solution to obtain the optimal solution, providing a better solution for the wide-area and normal glacier dynamic update.

[0071] In summary, a long-time series glacier flow velocity monitoring method proposed by the present invention effectively improves the problem that the significant intensity difference in SAR images under complex meteorological conditions of marine glaciers restricts the accurate monitoring of glacier flow velocity sequences, and provides an accurate data calculation method for glacier dynamic evolution investigation and trend assessment.

[0072] Specific embodiments are applied in the present invention to elaborate on the principle and implementation manner of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

[0073] Those of ordinary skill in the art will realize that the embodiments described herein are for helping readers understand the principle of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed in the present invention, and these deformations and combinations are still within the protection scope of the present invention.

Claims

1. A long-time series glacier flow velocity monitoring method, characterized in that, It includes the following steps: S1. Obtain a glacier SAR image sequence and perform preprocessing to generate a number of multi-looked SAR images; S2. Set a time baseline threshold. After combining image pairs of all multi-looked SAR images, perform offset estimation to respectively generate a number of offset sequences in the azimuth direction and the range direction; S3. Construct a non-linear glacier flow velocity model based on seasonal components, optimize all offset sequences, and generate optimized offset sequences; S4. Based on the optimized offset sequences, calculate the average flow velocity field on the glacier surface, and combine with the second glacier inventory data to generate a glacier flow velocity direction matrix; S5. Take the optimized offset sequences as observables, set the glacier flow direction as a constraint condition, and use the least squares method to iteratively solve the displacement distribution corresponding to the imaging time of each multi-looked SAR image to obtain the optimal displacement distribution corresponding to the imaging time of each multi-looked SAR image.

2. The long-time glacier flow velocity monitoring method according to claim 1, characterized in that Step S1 specifically includes: S11. Obtain a glacier SAR image sequence, which includes a number of glacier SAR images; S12. After mosaicking and sub-pixel registration of all glacier SAR images, perform multi-looking operations on the azimuth direction and the range direction of the registered glacier SAR images according to a proportionality coefficient to generate a number of multi-looked SAR images.

3. The long-time glacier flow velocity monitoring method according to claim 2, characterized in that, Step S2 specifically includes: S21. Set a time baseline threshold, combine image pairs of all multi-looked SAR images to generate a number of image pairs; S22. Use the optical flow motion estimation method to calculate the offset of each image pair to respectively generate offset sequences in the azimuth direction and the range direction, and finally obtain a number of offset sequences in the azimuth direction and the range direction.

4. The long-time glacier flow velocity monitoring method according to claim 3, wherein Step S3 specifically includes: S31. Construct a non-linear glacier flow velocity model based on seasonal components, which includes a periodic function and an aperiodic function; Among them, the periodic function is established by periodic main control factors including temperature and rainfall through sine and cosine functions; the aperiodic function is established by aperiodic main control factors including terrain, mass, and stress through a second-order polynomial function; S32. Based on all offset sequences, select an observation sample subset by setting a time threshold; S33. Based on the observation sample subset, fit the coefficients of the periodic function and the aperiodic function pixel by pixel using the non-linear least squares method; S34. Based on the fitted coefficients of the periodic function and the aperiodic function, calculate the flow velocity confidence interval of the offset sequences in the observation sample subset to filter the observation sample subset; S35. Determine whether there are still outliers in the filtered observation sample subset. If so, execute step S34. Otherwise, obtain the initial optimized coefficients of the periodic function and the aperiodic function of the non-linear glacier flow velocity model based on seasonal components; S36. Based on the initial optimized coefficients, calculate the flow velocity confidence interval of all offset sequences to filter the offset sequences, and determine whether all pixels of each offset sequence have been optimized. If so, generate optimized offset sequences and execute step S4. Otherwise, execute step S32.

5. The long-time glacier flow velocity monitoring method according to claim 4, characterized in that The formula of the non-linear glacier flow velocity model based on seasonal components in step S31 is: Among them, represents the glacier flow velocity at the imaging time , represents the sine function, represents the cosine function, , , , all represent the coefficients of periodic functions, , , all represent the coefficients of non-periodic functions.

6. The long-time series glacier flow velocity monitoring method according to claim 5, characterized in that Step S32 specifically includes: S321. Set the upper limit of the time threshold and the lower limit ; S322: Based on all offset sequences, select image pairs that satisfy the relationship to generate a subset of observation samples; Among them, represents the time interval of the -th group of image pairs.

7. The long-time glacier flow velocity monitoring method according to claim 6, characterized in that, The specific process of filtering the observation sample subset in step S34 is: First, define the set of median times for all image pairs in the observed sample subset as: , , , , , , represent the imaging times of the 1st, 2nd, 3rd, 4th, th, th multi-look SAR images respectively, represents the median time of the image pair composed of the 1st and 2nd multi-look SAR images, represents the median time of the image pair composed of the 1st and 3rd multi-look SAR images, represents the median time of the image pair composed of the 1st and 4th multi-look SAR images, represents the median time of the image pair composed of the 2nd and 3rd multi-look SAR images, represents the th and the th multi-look SAR images; Define the set of all image pairs as follows: , , , , , , respectively represent the image pair composed of the th and the th multi-look SAR images, the image pair composed of the 1st and the 2nd multi-look SAR images, the image pair composed of the 1st and the 3rd multi-look SAR images, the image pair composed of the 1st and the 4th multi-look SAR images, the image pair composed of the 2nd and the 3rd multi-look SAR images, the image pair composed of the th and the th multi-look SAR images; Secondly, obtain the instantaneous velocity corresponding to the median time of each image pair, and use it as the average flow velocity of each image pair to obtain the set of average flow velocities of all image pairs , that is: Among them, represents the instant corresponding to the median time of the image pair composed of the th and the th multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the 1st and the 2nd multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the 1st and the 3rd multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the 2nd and the 3rd multi-looked SAR images, represents the th and the th multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the Then, substitute the median time of each image pair and its corresponding instantaneous velocity into the non-linear glacier flow velocity model based on seasonal components. After fitting the model coefficients using the non-linear least squares method, calculate the upper limit of the velocity within the 95% confidence interval corresponding to the median time of the observed sample subset and the lower limit ; Finally, the upper speed limit corresponding to the median moment exceeding the 95% confidence interval of the observed sample subset and the lower limit of the outliers are removed to generate a filtered observed sample subset.

8. The long-time glacier flow velocity monitoring method according to claim 7, characterized in that Step S4 specifically includes: S41. Calculate the average flow velocity field on the glacier surface according to the optimized offset sequence, i.e.: Among them, represents the average flow velocity field on the glacier surface, represents the total number of optimized offset sequences, represents the th optimized offset sequence, represents the th time interval of the optimized offset sequence; S42. Obtain the glacier area and non - glacier area based on the second glacier inventory data, and divide the glacier area into a negative - direction area and a positive - direction area based on the average flow velocity field on the glacier surface, and finally generate a glacier flow velocity direction matrix.

9. The long-time glacier flow velocity monitoring method according to claim 8, characterized in that, Step S5 specifically includes: S501. Take the optimized offset sequence as the observable quantity, i.e.: Among them, represents the observable quantity, , , respectively represent the 1st, 2nd, th optimized offset sequences; S502. Take the displacement distribution corresponding to the imaging time of each multi - look SAR image as the quantity to be solved, i.e.: Among them, represents the quantity to be solved, , , respectively represent the displacement distributions corresponding to the imaging times of the first, second, and the -th multi-view SAR images; S503. Construct an observation equation for the displacement on the glacier surface according to the offset sequence of each image pair, i.e.: Among them, represents the offset sequence matrix optimized to the order, represents the coefficient matrix of the order, represents the matrix of quantities to be solved for Among them, the offset sequence of each image pair is as follows: , indicating the offset sequence of the image pair composed of the -th and the -th multi-look SAR images, indicating the displacement distribution corresponding to the imaging time of the -th multi-look SAR image, indicating the displacement distribution corresponding to the imaging time of the -th multi-look SAR image; S504. Perform singular value decomposition on the coefficient matrix of the observation equation for the displacement on the glacier surface to generate a coefficient matrix equation for singular value decomposition, i.e.: Among them, represents the coefficient matrix of singular value decomposition, represents a unitary matrix of order represents a non - negative diagonal matrix of order represents a unitary matrix of order represents the transpose; S505. Substitute the coefficient matrix equation for singular value decomposition into the observation equation for the displacement on the glacier surface, and at the same time take the glacier flow direction as the first constraint condition to generate an observation equation for the constrained flow velocity direction, i.e.: Among them, represents taking the minimum value, represents the optimized offset sequence matrix, represents the optimized offset sequence, represents the matrix of quantities to be settled, represents the matrix of glacier flow velocity directions, represents the th optimized offset sequence; S506. Use the least - squares method to perform the first - round solution on the observation equation for the constrained flow velocity direction to obtain the initial solution of the displacement distribution corresponding to the imaging time of each multi - look SAR image; S507. Convert the initial solution of the displacement distribution corresponding to the imaging time of each multi - look SAR image to generate the initial flow velocity of the displacement distribution corresponding to the imaging time of each multi - look SAR image, i.e.: Among them, represents the initial velocity set of the displacement distribution corresponding to the imaging time of each multi-look SAR image, represents the initial velocity of the displacement distribution corresponding to the imaging time of the first multi-look SAR image, represents the initial velocity of the displacement distribution corresponding to the imaging time of the second multi-look SAR image, represents the initial velocity of the displacement distribution corresponding to the imaging time of the S508. Substitute the initial flow velocity and imaging time of the displacement distribution corresponding to the imaging time of each multi - look SAR image into the non - linear glacier flow velocity model based on seasonal components for fitting estimation to generate the optimal coefficients of the periodic function and non - periodic function; S509. Based on the optimal coefficient, the upper limit of the 95% confidence interval corresponding to the initial flow velocity at the imaging time of each multi-look SAR image is calculated again using the non-linear glacier flow velocity model based on the seasonal component and the lower limit , to generate the second constraint condition; S510. Based on the second constraint condition, use the least - squares method to perform the second - round solution on the observation equation for the constrained flow velocity direction to obtain the optimal displacement distribution corresponding to the imaging time of each multi - look SAR image.

10. The long-time glacier flow velocity monitoring method according to claim 9, characterized in that The second constraint condition in step S510 is: Among them, , respectively represent the lower limit of the 95% confidence interval of the initial flow velocity corresponding to the imaging time of the th and th multi-look SAR images, , respectively represent the upper limit of the 95% confidence interval of the initial flow velocity corresponding to the imaging time of the th and th multi-look SAR images, represents the imaging time of the th multi-look SAR image, represents the imaging time of the th multi-look SAR image.

Images