A long-time series glacier flow velocity monitoring method

Abstract

The present invention relates to the technical field of image data processing, and discloses a long-time-series glacier flow velocity monitoring method, including: acquiring a glacier SAR image sequence and preprocessing it to generate a plurality of multi-looked SAR images; setting a time baseline threshold, after combining image pairs of all multi-looked SAR images, performing offset estimation to respectively generate a plurality of offset sequences in the azimuth direction and the range direction; constructing a non-linear glacier flow velocity model based on seasonal components, optimizing all offset sequences to generate optimized offset sequences, calculating the average flow velocity field on the glacier surface, and combining with the second glacier inventory data to generate a glacier flow velocity direction matrix; using the optimized offset sequences as observables, setting the glacier flow direction as a constraint condition, and using 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; this method improves the reliability and accuracy of glacier flow velocity calculation.

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 now 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-frequency flow velocity dynamic monitoring in the southeastern hinterland of the Qinghai-Tibet Plateau where marine glaciers are concentrated due to its advantage of being unaffected by clouds, rain, snow, and light and shadow. Utilizing the advantages of SAR remote sensing data acquisition can provide data support for revealing the evolution law of glacier change processes under the background of climate warming on a long time scale.

[0003] As one of the regions where modern glaciers are concentrated in the Hengduan Mountains, the glaciers in the Gongga Mountain area show an obvious trend of retreat, especially in the past ten-odd 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 significance for revealing the mechanism of glacier retreat in glaciers with similar hydrothermal conditions in southeastern Tibet. In addition, as a tourist resource, the Hailuogou Scenic Area on the eastern slope of Gongga Mountain attracts domestic and foreign tourists to stop and visit because of its beautiful glacier landscapes. 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 ice margin disaster-prone environment and dynamically evaluating glacier disaster risks.

[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 purpose, the technical solution adopted by the present invention is as follows:

[0007] A long-term glacier flow velocity monitoring method includes the following steps:

[0008] S1. Obtain a sequence of glacier SAR images and perform preprocessing to generate a number of multi-looked SAR images;

[0009] 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 range direction respectively;

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

[0011] 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;

[0012] S5. Use the optimized offset sequences as observables, set the glacier flow direction as a constraint condition, and adopt 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.

[0013] The present invention has the following beneficial effects:

[0014] A 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 filtering outliers 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 long-term inter-annual glacier flow velocity field, thereby improving the reliability and accuracy of long-term glacier flow velocity calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a schematic flow chart of a long-term glacier flow velocity monitoring method proposed by the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0016] The specific embodiments of the present invention will be described below 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 specific embodiments. For those ordinary skilled 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.

[0017] As Figure 1 shown, a long-time-series glacier flow velocity monitoring method includes the following steps S1 - S5:

[0018] S1. Obtain a glacier SAR image sequence and perform preprocessing to generate a number of multi-looked SAR images.

[0019] Specifically, step S1 specifically includes S11 - S12:

[0020] S11. Obtain a glacier SAR image sequence, which includes a number of glacier SAR images.

[0021] 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 a proportionality coefficient to generate a number of multi-looked SAR images.

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

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

[0024] Specifically, step S2 specifically includes S21 - S22:

[0025] S21. Set a time baseline threshold and combine image pairs of all multi-looked SAR images to generate a number of image pairs.

[0026] 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 rule 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.

[0027] S22. Use the optical flow motion estimation method to calculate the offsets for each pair of images, respectively generating offset sequences in the azimuth direction and the range direction, and finally obtaining a number of offset sequences in the azimuth direction and the range direction.

[0028] In this embodiment, the optical flow motion estimation method is used to calculate the offsets for each pair of images. Each pair of images will respectively generate a set of offset results in the azimuth direction and the range direction, so as to generate offset sequences in the azimuth direction and the range direction respectively, and finally generate a number of offset sequences in the azimuth direction and the range direction.

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

[0030] Specifically, step S3 specifically includes S31 - S36:

[0031] 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.

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

[0033]

[0034] Where, 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.

[0035] In this embodiment, by assuming the glacier flow velocity and decomposing it based on the periodic main control factors and the 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 , the non-periodic main control factors such as terrain, mass, and stress in the model can be represented by a second-order polynomial function, that is .

[0036] S32. Based on all offset sequences, by setting a time threshold, select a subset of observation samples, specifically:

[0037] S321. Set the upper limit of the time threshold and the lower limit .

[0038] S322: Based on all offset sequences, select the image pairs that satisfy relationship to generate a subset of observation samples; where represents the time interval of the th group of image pairs.

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

[0040] S33. Based on the subset of observation samples, use the non-linear least squares method to fit the coefficients of the periodic function and the non-periodic function for each pixel.

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

[0042] Specifically, the specific process of filtering the subset of observation samples in step S34 is as follows:

[0043] First, define the set of the median times of all image pairs in the subset of observation samples as: , , [[ID=I47]] , , , , represent the imaging times of the 1st, 2nd, 3rd, 4th, , 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 Zhang and the median time of the image pair composed of Zhang multi-look SAR images.

[0044] Define the set of all image pairs as: , , , , , , respectively represent the Zhang and the Zhang multi-look SAR image pair, the first and second Zhang multi-look SAR image pair, the first and third Zhang multi-look SAR image pair, the first and fourth Zhang multi-look SAR image pair, the second and third Zhang multi-look SAR image pair, the Zhang and the Zhang multi-look SAR image pair.

[0045] 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:

[0046]

[0047] Among them, represents the Zhang and the instantaneous corresponding to the median time of the image pair composed of Zhang multi-look SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the first and second Zhang multi-look SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the first and third Zhang multi-look SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the second and third Zhang multi-look SAR images, represents the Zhang and the instantaneous velocity corresponding to the median time of the image pair composed of Zhang multi-look SAR images.

[0048] Then, substitute the median time of each image pair and its corresponding instantaneous velocity into the non-linear glacier flow velocity model based on the seasonal component, and after fitting the model coefficients using the non-linear least squares method, calculate the upper limit and lower limit of the velocity at the 95% confidence interval corresponding to the median time of the observed sample subset.

[0049] In this embodiment, the periodic function and non-periodic function coefficients of the model are estimated by the least squares method. The relationship between the median time of each image pair and its corresponding instantaneous velocity is quantified using these coefficients, and the uncertainty range of the estimated model coefficients, i.e., the upper limit of the instantaneous velocity corresponding to the median time, is measured by the 95% confidence interval. and the lower limit , which provides a basis for filtering outliers in the subsequent steps.

[0050] Finally, the velocity upper limit and the lower limit of the 95% confidence interval corresponding to the median time that exceed the observed sample subset are removed, and a filtered observed sample subset is generated.

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

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

[0053] 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.

[0054] In summary, this step constructs a nonlinear glacier flow velocity model based on the seasonal component to optimize the offset sequences, thereby improving the reliability of the glacier flow velocity.

[0055] 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.

[0056] Specifically, step S4 includes S41 - S42:

[0057] S41. Calculate the average flow velocity field on the glacier surface according to the optimized offset sequences, that is:

[0058]

[0059] where, 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 The time intervals of an optimized offset sequence.

[0060] 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, th optimized offset sequence; , , respectively represent the time intervals of the 1st, 2nd, th optimized offset sequence.

[0061] S42. According to the second glacier inventory data, the glacier area and non-glacier area are obtained, and 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.

[0062] 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, and 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 the iterative solution in the subsequent steps.

[0063] S5. Taking the optimized offset sequence as the observation 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.

[0064] Specifically, step S5 specifically includes S501 - S510:

[0065] S501. Taking the optimized offset sequence as the observation quantity, that is:

[0066]

[0067] Among them, represents the observation quantity, , , respectively represent the 1st, 2nd, th optimized offset sequence.

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

[0069]

[0070] Among them, represents the quantity to be solved, , , respectively represent the displacement distributions corresponding to the imaging times of the first, second, and th multi-view SAR images.

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

[0072]

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

[0074] Among them, the offset sequences of each image pair are: , represents the th and th multi-view SAR images, 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 .

[0075] 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 sequences as the observed quantities, and constructing an observation equation composed of observed quantities and unknown quantities to be solved, 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 th column coefficients of each row of the coefficient matrix are 1 and -1 respectively, and the rest are all 0, that is:

[0076] .

[0077] S504. Perform singular value decomposition on the coefficient matrix of the observation equation for the glacier surface displacement to generate a coefficient matrix equation for singular value decomposition, i.e.:

[0078]

[0079] where, represents the coefficient matrix of singular value decomposition, represents an order unitary matrix, represents a non - negative diagonal matrix of order, represents an order unitary matrix, represents the transpose.

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

[0081] S505. Substitute the coefficient matrix equation of singular value decomposition into the observation equation for the glacier surface displacement, and at the same time use the glacier flow direction as the first constraint condition to generate an observation equation for constraining the flow velocity direction, i.e.:

[0082]

[0083] where, represents taking the minimum value, 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.

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

[0085] 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.:

[0086]

[0087] where, An initial velocity set representing the displacement distribution corresponding to the imaging time of each multi-look SAR image The initial velocity representing the displacement distribution corresponding to the imaging time of the first multi-look SAR image The initial velocity representing 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 multi-look SAR image

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

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

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

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

[0092]

[0093] Among them, , Respectively represent the lower limits of the 95% confidence intervals corresponding to the initial velocities corresponding to the imaging times of the th and th multi-look SAR images , Respectively represent the upper limits of the 95% confidence intervals corresponding to the initial velocities corresponding to the imaging times 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

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

[0095] S510. Based on the second constraint condition, the least squares method is used to perform a second-round solution for the observation equation that constrains the flow velocity direction, and the optimal displacement distribution corresponding to the imaging time of each multi-look SAR image is obtained.

[0096] In this embodiment, the initial solution of the displacement distribution corresponding to the imaging time of each multi-look 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 the seasonal component together with the imaging time for calculation to calculate the upper and lower limits of the velocity corresponding to the imaging time, and this is used as the displacement gradient constraint for the second-round solution to obtain the optimal solution, providing a better solution for wide-area and normalized glacier dynamic update.

[0097] 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 of SAR images under complex meteorological conditions of marine glaciers restricts the accurate monitoring of the glacier flow velocity sequence, and provides an accurate data calculation method for glacier dynamic evolution investigation and trend assessment.

[0098] In the present invention, specific embodiments are used 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.

[0099] Those of ordinary skill in the art will realize that the embodiments described herein are for helping the reader 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 specific deformations and combinations that do not deviate from the essence of the present invention according to these technical revelations disclosed by 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 all the multi-looked SAR images into image pairs, 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 the 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. Specifically: S41. Calculate the average flow velocity field on the glacier surface according to the optimized offset sequences, that is: 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 according to the second glacier inventory data, and divide the glacier area into a negative area and a positive area based on the average flow velocity field on the glacier surface, and finally 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. Specifically: S501. Take the optimized offset sequences as observables, that is: Among them, represents the observed quantity, , , respectively represent the 1st, 2nd, th optimized offset sequences; S502. Take the displacement distribution corresponding to the imaging time of each multi-looked SAR image as the quantity to be solved, that is: Among them, represents the quantity to be solved, , , respectively represent the displacement distributions corresponding to the imaging times of the first, second, th multi-view SAR images; S503. According to the offset sequences of each image pair, construct an observation equation for the displacement on the glacier surface, that is: Among them, represents the offset sequence matrix optimized to the th order, represents the coefficient matrix of the th order, represents the matrix of the quantity to be solved of the Among them, the offset sequence of each image pair is as follows: , denotes the offset sequence of the image pair composed of the th and the th multi-look SAR images, denotes the displacement distribution corresponding to the imaging time of the th multi-look SAR image, denotes 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 of the singular value decomposition, that is: Among them, represents the coefficient matrix of singular value decomposition, represents an $n$-order unitary matrix, represents an $m$-order non-negative diagonal matrix, represents a $p$-order unitary matrix, represents the transpose; S505. Substitute the coefficient matrix equation of the 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, that is: Among them, represents taking the minimum value, represents the optimized offset sequence, represents the matrix of quantities to be solved, 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 of 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-looked SAR image; S50�. Convert the initial solution of the displacement distribution corresponding to the imaging time of each multi-looked SAR image to generate the initial flow velocity of the displacement distribution corresponding to the imaging time of each multi-looked SAR image, that is: 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 the imaging time of the displacement distribution corresponding to the imaging time of each multi-looked 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 the 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-view 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 of 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-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 the glacier SAR images, perform multi-looking operations on the azimuth direction and the range direction of the registered glacier SAR images according to the scale factor to generate a number of multi-looked SAR images.

3. The long-time glacier flow velocity monitoring method according to claim 2, wherein Step S2 specifically includes: S21. Set the time baseline threshold, combine all multi-look SAR images to generate several image pairs; S22. Use the optical flow motion estimation method to calculate the offsets for each image pair, respectively generating offset sequences in the azimuth direction and range direction, and finally obtaining several offset sequences in the azimuth direction and 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 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; S32. Based on all the offset sequences, select a subset of observation samples by setting a time threshold; S33. Based on the subset of observation samples, 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 coefficients of the fitted periodic function and aperiodic function, calculate the velocity confidence intervals of the offset sequences in the subset of observation samples to filter the subset of observation samples; S35. Determine whether there are still outliers in the filtered subset of observation samples. If so, execute step S34. Otherwise, obtain the initial optimized coefficients of the periodic function and aperiodic function of the non-linear glacier flow velocity model based on seasonal components; S36. Based on the initial optimized coefficients, calculate the velocity confidence intervals of all the offset sequences to filter the offset sequences, and determine whether all the pixels of each offset sequence have been optimized. If so, generate the 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 moment , 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, wherein 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 an observation sample subset; Among them, represents the time interval between the 7. The long-time glacier flow velocity monitoring method according to claim 6, wherein The specific process of filtering the subset of observation samples in step S34 is: First, define the set of median times for all image pairs in the observed sample subset as follows: , , , , , , respectively represent the imaging times of the 1st, 2nd, 3rd, 4th, th, 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 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 moment 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 2nd multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the 1st and 3rd multi-looked SAR images, represents the instantaneous velocity corresponding to the median time of the image pair composed of the 2nd and 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 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 of the 95% confidence interval corresponding to the median time beyond the observed sample subset and the lower limit of the outliers are removed to generate a filtered subset of the observed samples.

8. The long-time glacier flow velocity monitoring method according to claim 1, wherein 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 the -th multi-looked 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 the -th multi-looked SAR images, represents the imaging time of the -th multi-looked SAR image, represents the imaging time of the -th multi-looked SAR image.

Images