CN114200449A - 基于qr分解的变量投影法在地表形变中的应用 - Google Patents

Abstract

本发明公开了一种基于QR分解的变量投影法在地表形变中的应用方法，将研究区在时间与空间上的地表形变拟合作为研究重点，在应用SBAS‑InSAR时序分析技术探究研究区地表形变的基础上，采用一种具有可分离结构的复杂指数模型拟合地表在空间位置上的沉降形势，利用基于QR分解的变量投影算法解算模型参数，在空间位置上以高精度逼近沉降剖面曲线，同时结合Logistic模型进行了地表在时间上的预测，实现了对地表形变状况在时间和空间上的全面监测，为地表形变监测提供了可靠的技术支持。

Description

基于QR分解的变量投影法在地表形变中的应用

技术领域

本发明属于卫星遥感技术领域，特别涉及基于QR分解的变量投影法的地表形变应用方法。

背景技术

传统监测地表形变的方法虽然精度较高，但只能监测有限个地面点的形变而无法大面积地监测地表整体形变情况，因此具有一定的局限性。合成孔径雷达差分干涉测量技术(Differential Interferometric Synthetic Aperture Radar,D-InSAR)监测地面形变具有全天时、大范围的优势，打破了传统方法的局限性(李珊珊等,2013；Gabrile A.K.etal.,1989；Massonnet D.et al.,1993)。然而D-InSAR技术采用的SAR图像数量较少，缺乏多余观测，且当时间跨度较大时容易产生失相干问题，难以获得连续的形变场。为克服这一问题，Berardino等人在D-InSAR的基础上提出了SBAS-InSAR时序分析技术，通过对大量数据中高相干点目标的时序分析，获取研究区的地表形变速率和时间形变序列(BerardinoP.et al.,2002)。该项技术是一种基于最小二乘意义上的数据后处理技术，可有效提高结果的可靠性，从而可以应用于监测地表长时间缓慢形变(王超等,2002)。

Golub和Pereyra于1973年针对可分离非线性最小二乘问题提出变量投影法(Variable Projection,VP)(Golub G H et al.,1973)，此后许多学者对这种将变量分离后再进行解算的算法进行了深入研究(Kaufman L.,1975；Ruhe A.et al.,1980；CorradiC.,1981；Dianne P.O’Leary et al.,2013)，此算法在数值分析、神经网络、电子通信与工程、环境科学和时间序列分析等领域得到了广泛的发展和应用。学者们在不同科研领域对变量投影法进行了多次尝试，并且取得了较为丰富的研究成果。根据实测资料预测地表沉降常用模型为Logistic曲线方程，通过大量工程实例该模型已被验证其科学性和合理性(朱志铎等,2009)。然而，Logistic模型多数用来预测时间序列上的沉降过程而无法拟合空间位置上的地表状态。

发明内容

本发明的目的在于克服上述缺陷，提供一种基于QR分解的变量投影法在地表形变中的应用方法，将研究区在时间与空间上的地表形变拟合作为研究重点，在应用SBAS-InSAR时序分析技术探究研究区地表形变的基础上，采用一种具有可分离结构的复杂指数模型拟合地表在空间位置上的沉降形势，利用基于QR分解的变量投影算法解算模型参数，在空间位置上以高精度逼近沉降剖面曲线，同时结合Logistic模型进行了地表在时间上的预测，实现了对地表形变状况在时间和空间上的全面监测，为地表形变监测提供了可靠的技术支持。

为实现上述发明目的，本发明提供如下技术方案：

一种基于QR分解的变量投影法在地表形变中的应用，包括以下步骤：

S1利用SBAS-InSAR时序分析方法，通过设定时间阈值生成匹配像对，进而获取原始地表沉降数据；

S2利用原始地表沉降数据，基于QR分解的变量投影算法进行空间位置上的地表沉降拟合，具体方法为：

S2.1将原始地表沉降数据代入可分离指数模型得到可分离指数模型残差函数；所述可分离指数模型残差函数为可分离非线性方程组形式；

S2.2在最小二乘原则下将可分离指数模型残差函数最小化，利用LM迭代求得可分离指数模型残差函数中的非线性参数估值

S2.3计算可分离指数模型残差函数中的线性参数估值

S2.4将非线性参数估值

和线性参数估值

代入可分离指数模型，得到解算完成的可分离指数模型，根据解算完成的可分离指数模型绘制地表沉降剖面图和沉降拟合图；

S3利用原始地表沉降数据，基于Logistic模型进行时间序列中的地表沉降拟合，具体方法为：

S3.1将原始地表沉降数据代入Logistic模型得到Logistic模型残差函数；所述Logistic模型残差函数为非线性方程组形式；

S3.2在最小二乘原则下将Logistic模型残差函数最小化，利用LM迭代求得Logistic模型残差函数中的参数估值；

S3.3将参数估值代入Logistic模型，得到解算完成的Logistic模型，根据解算完成的Logistic模型绘制地表沉降随时间的变化曲线和曲线拟合图。

进一步的，所述步骤S1中，利用SBAS-InSAR时序分析方法获取原始地表沉降数据的具体方法为：

S1.1通过设定时间阈值对遥感影像进行匹配，得到匹配像对；

S1.2对匹配像对进行差分干涉，生成干涉像对；

S1.3对由干涉像对组成的干涉图依次进行滤波增强处理，相位解缠后，根据干涉图估算地表形变速率和残余地形，并根据地表形变速率和残余地形计算时间序列上的地表形变量；

S1.4在时间序列上的地表形变量中去除大气相位，并转换为地理坐标系后，得到原始地表沉降数据。

进一步的，所述步骤S2.1中，可分离指数模型为：

可分离指数模型残差函数为：r＝y-η＝y-Φβ；

其中，y为地表沉降值，η为包含空间位置变量u，非线性参数α以及线性参数β的地表沉降估值，非线性函数矩阵

r为可分离指数模型残差，i表示模型中非线性函数的项数，i＝1，2，……n；非线性参数α包含了

λ和z，α＝(λ₁,λ₂,…,λ_n；z)，β＝(β₀,β₁,…,β_n)。

进一步的，所述步骤S2.2中，利用LM迭代求得可分离指数模型残差函数中的非线性参数估值

的方法为：

S2.2.1对可分离指数模型残差函数求解关于非线性参数α的一阶偏导，得到第一雅可比矩阵J；

S2.2.2给定非线性参数初值α₀和第一阻尼因子调整系数τ，设置收敛准则：第一迭代终止阈值ε和第一最大迭代次数k_max；

S2.2.3将α₀代入第一雅可比矩阵J和非线性函数矩阵Φ，对非线性函数矩阵Φ进行QR分解，得到基于QR分解的第一变量投影函数

令第一阻尼因子μ＝τ·max{diag(J^TJ)}，计算第一迭代步长h_LM＝-(J^TJ+μI)^-1J^Tr，α_k+1＝α_k+h_LM，k≥0；

S2.2.4如果h_LM使得

则接受α_k+1并减小第一阻尼因子μ，进入下一步；否则，增大第一阻尼因子μ并重新计算h_LM，直至

S2.2.5若

或k＞k_max，则终止迭代，输出此时的α_k+1，即为非线性参数估值

否则，重复执行步骤S2.2.3和S2.2.4。

进一步的，所述步骤S2.2.3中，对非线性函数矩阵Φ进行QR分解，得到基于QR分解的变量投影函数

的方法为：

将Φ进行QR分解得到分解式Φ＝QR，令Q＝[Q₁Q₂]，Q₁、Q₂分别为m×p和m×(m-o)阶矩阵，m为地表沉降数据个数，o为Φ的秩，R为上三角矩阵；

令线性参数β＝Φ⁺y，其中Φ⁺是Φ的广义逆矩阵；

结合分解式和β＝Φ⁺y得到非线性函数矩阵Φ基于QR分解的变量投影函数

进一步的，所述步骤S2.3中，线性参数估值

其中

为将非线性参数估值

代入后所得非线性函数矩阵的广义逆矩阵。

进一步的，所述步骤S3.1中，Logistic模型为

Logistic模型残差函数为r′＝y-S；其中，y为地表沉降值，S为包含时间变量t和参数A₁,A₂,A₃的地表沉降估值，r′为Logistic模型残差。

进一步的，所述步骤S3.2中，利用LM迭代求得Logistic模型残差函数中的参数估值的方法为：

S3.2.1对Logistic模型残差函数求解关于参数A的一阶偏导，得到第二雅可比矩阵J′，所述参数A为包含参数A₁,A₂和A₃的参数；

S3.2.2给定参数初值A₀和第二阻尼因子调整系数τ′，设置收敛准则：第二迭代终止阈值ε′和第二最大迭代次数k′_max；

S3.2.3将A₀代入第二雅可比矩阵J′和r′，进而得到残差平方和

令第二阻尼因子μ′＝τ′·max{diag(J′^TJ′)}，计算第二迭代步长h′_LM＝-(J′^TJ′+μ′I)^-1J′^Tr′，A_k+1＝A_k+h_LM，k≥0；

S3.2.4如果h′_LM使得

则接受A_k+1并减小第二阻尼因子μ′，进入下一步；否则，增大第二阻尼因子μ′并重新计算h′_LM，直至

S3.2.5若

或k＞k′_max，则终止迭代，输出此时的A_k+1，即为参数A的估值

否则，重复执行步骤S3.2.3和S3.2.4。

进一步的，所述步骤S1.2中，利用以下公式对对匹配像对进行差分干涉，生成干涉像对：

δφ_j(p,q)＝φ(t_B,p,q)-φ(t_A,p,q)；

其中，t_A,t_B为某两景SAR影像的获取时间，(p,q)表示像素点位置，φ表示对应像素点的相位，δφ表示对应像素点的相位差，j为在t_A和t_B时刻获得的两景SAR影像干涉后的干涉像对序列号。

本发明与现有技术相比具有如下有益效果：

(1)本发明在以SBAS-InSAR方法得到的地表形变数据基础上，采用QR分解的变量投影法和一种指数模型，实现了地表形变状况在空间上的监测，同时并配合Logistic模型对沉降区进行时间上的拟合，实现了对地表形变状况在时间上的监测，监测结果精确全面；

(2)本发明利用基于QR分解的变量投影算法，将一种具有可分离结构的指数模型中的线性与非线性参数分离开来依次求解，对研究区内较为典型的沉降区进行地理位置上的沉降剖面拟合，解算精度较高，且在顾及精度的同时提高了计算效率；

(3)本发明以青海祁连县为研究区对技术方案进行了验证，对指数模型和Logistic模型采用LM算法求解，具有较高的拟合精度。

附图说明

图1为本发明实施例1中祁连县范围及周边地理环境图；

图2为本发明实施例1中时间序列地表形变图与年平均形变速率图，其中，图2(1)～(6)为时间序列地表形变图，图2(7)为年平均形变速率图；

图3为本发明实施例1中终期地表形变图；

图4为本发明实施例1中可分离指数模型拟合沉降剖面图；

图5为本发明实施例1中沉降剖面预测图；

图6为本发明实施例1中Logistic模型拟合时间序列沉降值。

具体实施方式

下面通过对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。

本发明针对D-InSAR技术缺乏多余观测，和因时间跨度较大而产生的失相干问题，将SBAS-InSAR时序分析技术应用于地表形变监测中，基于不同时期的遥感影像数据，进行地表沉降变化监测，并将一种可分离非线性模型、Logistic模型和基于QR分解的变量投影法应用于沉降剖面拟合，为监测地表形变提供技术支撑。

本发明一种基于QR分解的变量投影法在地表形变中的应用，包括如下步骤：

(1)将SBAS-InSAR时序分析方法应用在地表形变监测中，掌握地表沉降状况与规律；本步骤基于不同时期的遥感数据选取多幅影像，根据时间和空间基线阈值配对连接产生干涉，并进行滤波和相位解缠等工作，最后计算地表形变量和形变速率；

本步骤具体方法为：

通过设定时间阈值对遥感影像进行匹配，匹配像对用线连接生成时空基线连接图。运用公式(6)对匹配像对做差分干涉，生成干涉像对：

δφ_j(p,q)＝φ(t_B,p,q)-φ(t_A,p,q) (6)。

式中，其中，t_A,t_B为某两景SAR影像的获取时间，(p,q)表示像素点位置，φ表示对应像素点的相位，δφ表示对应像素点的相位差，j为在t_A和t_B时刻获得的两景SAR影像干涉后的干涉像对序列号。对由干涉像对组成的干涉图进行滤波增强处理以提高干涉质量，并进行相位解缠以提取地表形变信息，估算形变速率和残余地形，计算时间序列上的形变量，估算并去除大气相位，并将结果转换为地理坐标系，得到原始地表沉降数据。

(2)利用具有可分离结构的指数函数模型，将基于QR分解的变量投影算法应用于模型的线性与非线性参数求解，对地表在空间位置上的沉降剖面进行拟合；

(2.1)将原始地表沉降数据代入可分离指数模型得到如式(1)所示的残差函数，形成可分离非线性方程组，可分离指数模型表达式为：

可分离指数模型残差函数为：r＝y-η＝y-Φβ； (1)；

其中，y为地表沉降值，η为包含空间位置变量u，非线性参数α以及线性参数β的地表沉降估值，非线性函数矩阵

r为可分离指数模型残差，i表示模型中非线性函数的项数，i＝1，2，……n；非线性参数α包含了非线性参数λ和z，α＝(λ₁,λ₂,…,λ_n；z)，β＝(β₀,β₁,…,β_n)。

(2.2)在最小二乘原则下将可分离指数模型残差函数最小化，利用LM迭代求得非线性参数估值

依据最小二乘原则，线性与非线性参数的求解即求使得可分离指数模型残差函数的平方和最小的

和

即：

将矩阵Φ通过QR分解法分解为：Φ＝QR，令Q＝[Q₁Q₂]，Q₁、Q₂分别为m×p和m×(m-o)阶矩阵，m为地表沉降数据个数，o为Φ的秩，R为上三角矩阵；令式(2)中的线性参数β＝Φ⁺y，结合Φ的QR分解式得到变量投影函数如下：

式中，Φ⁺是Φ的广义逆矩阵。

利用LM迭代求得非线性参数估值

的具体过程如下：

(2.2.1)对残差函数求关于非线性参数的一阶偏导，得到第一雅可比矩阵J；

(2.2.2)给定非线性参数初值α₀，第一阻尼因子调整系数τ，设置收敛准则——第一迭代终止阈值ε和第一最大迭代次数k_max；

(2.2.3)代入初值α₀得到数值型矩阵J和Φ，对Φ进行QR分解，进而得到如式(3)的基于QR分解的第一变量投影函数

令第一阻尼因子μ＝τ·max{diag(J^TJ)}，计算第一迭代步长h_LM＝-(J^TJ+uI)^-1J^Tr，α_k+1＝α_k+h_LM；

(2.2.4)若h_LM使得

则接受α_k+1并减小第一阻尼因子，进入下一步，否则，增大第一阻尼因子并重新计算h_LM，直至使得

(2.2.5)若

或k＞k_max，则终止迭代，否则，重复执行(2.2.3)和(2.2.4)。

(2.3)线性参数β按下式求解：β＝Φ⁺y (4)；

(2.4)将线性和非线性参数估值代入可分离指数模型，绘制地表沉降剖面图和沉降拟合图。

(3)利用Logistic模型对地表随时间推移的沉降状况进行拟合；

(3.1)将原始地表沉降数据代入Logistic模型得到Logistic模型残差函数，形成非线性方程组；

Logistic模型表达式为

Logistic模型残差函数为r′＝y-s (5)；

其中，y为地表沉降值，s为包含时间变量t和参数A₁,A₂,A₃的地表沉降估值，r′为Logistic模型残差。

(3.2)在最小二乘原则下将可Logistic模型残差函数最小化，利用LM迭代求得参数估值：

(3.2.1)对Logistic模型残差函数求解关于参数A的一阶偏导，得到第二雅可比矩阵J′；所述参数A为包含参数A₁,A₂和A₃的参数

(3.2.2)给定参数初值A₀和第二阻尼因子调整系数τ′，设置收敛准则：迭代第二终止阈值ε′和第二最大迭代次数k′_max；

(3.2.3)代入初值A₀得到数值型矩阵J′和r′，得到残差平方和

令第二阻尼因子μ′＝τ′·max{diag(J′^TJ′)}，计算第二迭代步长h_L′_M＝-(J′^TJ′+μ′I)^-1J′^Tr′，A_k+1＝A_k+h_L′_M，k≥0；

(3.2.4)如果h′_LM使得

则接受A_k+1并减小第二阻尼因子μ′，进入下一步；否则，增大第二阻尼因子μ′并重新计算h′_LM，直至

(3.2.5)若

或k＞k′_max，则终止迭代，输出此时的A_k+1，即为参数A的估值

否则，重复执行步骤(3.2.3)和(3.2.4)。

(3.3)将参数估值代入Logistic模型，绘制地表沉降随时间的变化曲线和曲线拟合图。

实施例1：

一、SBAS-InSAR技术获取形变信息

案例以青海祁连县为研究区，祁连县及周边地理环境如图1所示。所用SAR数据为7幅Sentinel-1升轨数据，C波段，VV极化方式，影像覆盖范围在37°55′N～38°25′N，99°40′E～100°30′E，时间跨度为2019年2月19日至2020年6月25日。影像数据时间信息如表1所示。基于SBAS-InSAR技术对数据进行了处理，获取了研究区在时间序列上的6幅地表形变图及1幅年平均形变速率图，如图2所示。

表1实验数据时间信息

在地表形变图和年平均形变速率图中，抬升区域用正值表示，沉降区域用负值表示。从图2中的年平均形变速率图中可知，城市及周边地区抬升速率较快，山区部分主要表现为沉降，且大多数区域地表形变速率集中在0—20mm/年，沉降速率较慢，表明地表相对比较稳定。观察时间序列形变图可知，自2019年12月起，山区地表开始产生明显形变，且随着时间的推移，形变量逐渐增大。由图2(6)地表形变图，研究区中的沉降区域与抬升区域划分明显，城市地区主要表现为抬升，而山区主要为沉降，这与年平均形变速率图所得结论较为一致。研究区最大沉降值超过了230mm，但多数沉降区沉降值集中在[-60,0]区间内。为更直观地显示研究区地表沉降情况，在研究区中选取一线段QL贯穿了祁连县城区部分，线段两端点地理坐标分别为(38.201°N，100.231°E)、(38.162°N，100.264°E)，如图3中所示，提取线段上各点的沉降值，利用可分离指数模型进行空间位置的沉降拟合与分析，绘制沉降剖面图，并利用Logistic模型拟合该线段上一点在时间序列上的沉降值。

二、基于QR分解的变量投影法在空间位置中的地表沉降拟合

实验利用如下可分离指数模型拟合沉降曲线：

其中，α＝(λ₁,λ₂,…,λ_n；z₁,z₂,…,z_n)^T，

此模型在曲线拟合中具有优良的性能，且具有线性参数与非线性参数可分离的结构，因此可以利用变量投影算法进行求解。令可分离指数模型中的n＝4，则可分离指数模型中有9个非线性参数及5个线性参数。采用变量投影算法拟合线段QL上提取的前300组沉降值数据并对地表形变进行预测，利用第301～400组数据对预测结果进行验证。将数据代入可分离指数模型，形成非线性方程组，给定相同的非线性参数迭代初值0.5，基于本发明的QR分解变量投影算法(VP_QR)拟合地表沉降方法、现有技术中的奇异值分解变量投影算法(VP_SVD)和满秩分解变量投影算法(VP_FRD)，使用LM对参数进行迭代求解。线性参数通过公式(4)求得。

可分离指数模型拟合各点沉降值结果如图4所示。表2列出了各算法对应的迭代次数、残差2-范数(r_2-norm)、均方根误差(RMSE)以及运算时间。由图4可知，线段穿过的祁连县区域地表沉降较小，最大沉降约为0.07mm且部分区域有抬升现象。从图中可以判断VP_FRD解算结果较差，而VP_QR与VP_SVD算法对应曲线几近重合，与原始数据拟合程度较高，均方根误差值数量级为10^-6，各点形变值与估值误差极小，表明参数解算精度较高，模型具有很好的拟合优度。在解算精度相当的前提下，VP_QR算法有最短的运算时间，计算效率更高。由此可见，本发明采用可分离指数模型与VP_QR算法在拟合空间位置上的地表沉降剖面具有良好的适用性。

表2迭代次数、残差2-范数(r_2-norm)、均方根误差(RMSE)及运算时间

可分离指数模型对各点沉降值的预测结果如图5所示。表3列出了各算法在预测过程中对应残差2-范数(r_2-norm)及均方根误差(RMSE)。由图5可知，VP_FRD预测效果较差，而VP_QR与VP_SVD算法对应的预测曲线与原始沉降数据具有一致的趋势，其预测精度为10^-6，形变值与其预测值误差极小。

表3模型预测的残差范数(r_2-norm)及均方根误差(RMSE)

三、Logistic模型在时间序列中的沉降拟合

选取线段QL上第417个像元，提取该像元每一时期的沉降值，得到6个沉降数据。Logistic模型如式(5)，将6个沉降值代入Logistic模型，给定Logistic模型参数的迭代初值，使用LM迭代算法进行求解，拟合该点各时期沉降值，如图6所示。从图中可以判断曲线具有较好的拟合效果，Logistic模型不仅适用于拟合时序列上的地表沉降，还可用于预测地面点随时间推移的沉降趋势。

以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (9)

1.一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，包括以下步骤：

S1利用SBAS-InSAR时序分析方法，通过设定时间阈值生成匹配像对，进而获取原始地表沉降数据；

S2利用原始地表沉降数据，基于QR分解的变量投影算法进行空间位置上的地表沉降拟合，具体方法为：

S2.1将原始地表沉降数据代入可分离指数模型得到可分离指数模型残差函数；所述可分离指数模型残差函数为可分离非线性方程组形式；

S2.2在最小二乘原则下将可分离指数模型残差函数最小化，利用LM迭代求得可分离指数模型残差函数中的非线性参数估值

S2.3计算可分离指数模型残差函数中的线性参数估值

S2.4将非线性参数估值

和线性参数估值

代入可分离指数模型，得到解算完成的可分离指数模型，根据解算完成的可分离指数模型绘制地表沉降剖面图和沉降拟合图；

S3利用原始地表沉降数据，基于Logistic模型进行时间序列中的地表沉降拟合，具体方法为：

S3.1将原始地表沉降数据代入Logistic模型得到Logistic模型残差函数；所述Logistic模型残差函数为非线性方程组形式；

S3.2在最小二乘原则下将Logistic模型残差函数最小化，利用LM迭代求得Logistic模型残差函数中的参数估值；

S3.3将参数估值代入Logistic模型，得到解算完成的Logistic模型，根据解算完成的Logistic模型绘制地表沉降随时间的变化曲线和曲线拟合图。

2.根据权利要求1所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S1中，利用SBAS-InSAR时序分析方法获取原始地表沉降数据的具体方法为：

S1.1通过设定时间阈值对遥感影像进行匹配，得到匹配像对；

S1.2对匹配像对进行差分干涉，生成干涉像对；

S1.3对由干涉像对组成的干涉图依次进行滤波增强处理，相位解缠后，根据干涉图估算地表形变速率和残余地形，并根据地表形变速率和残余地形计算时间序列上的地表形变量；

S1.4在时间序列上的地表形变量中去除大气相位，并转换为地理坐标系后，得到原始地表沉降数据。

3.根据权利要求1或2所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S2.1中，可分离指数模型为：

可分离指数模型残差函数为：r＝y-η＝y-Φβ；

其中，y为地表沉降值，η为包含空间位置变量u，非线性参数α以及线性参数β的地表沉降估值，非线性函数矩阵

r为可分离指数模型残差，i表示模型中非线性函数的项数，i＝1，2，……n；非线性参数α包含了λ和z，α＝(λ₁,λ₂,…,λ_n；z)，β＝(β₀,β₁,…,β_n)。

4.根据权利要求3所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S2.2中，利用LM迭代求得可分离指数模型残差函数中的非线性参数估值

的方法为：

S2.2.1对可分离指数模型残差函数求解关于非线性参数α的一阶偏导，得到第一雅可比矩阵J；

S2.2.2给定非线性参数初值α₀和第一阻尼因子调整系数τ，设置收敛准则：第一迭代终止阈值ε和第一最大迭代次数k_max；

S2.2.3将α₀代入第一雅可比矩阵J和非线性函数矩阵Φ，对非线性函数矩阵Φ进行QR分解，得到基于QR分解的第一变量投影函数

令第一阻尼因子μ＝τ·max{diag(J^TJ)}，计算第一迭代步长h_LM＝-(J^TJ+μI)^-1J^Tr，α_k+1＝α_k+h_LM，k≥0；

S2.2.4如果h_LM使得

则接受α_k+1并减小第一阻尼因子μ，进入下一步；否则，增大第一阻尼因子μ并重新计算h_LM，直至

S2.2.5若

或k＞k_max，则终止迭代，输出此时的α_k+1，即为非线性参数估值

否则，重复执行步骤S2.2.3和S2.2.4。

5.根据权利要求4所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S2.2.3中，对非线性函数矩阵Φ进行QR分解，得到基于QR分解的变量投影函数

的方法为：

将Φ进行QR分解得到分解式Φ＝QR，令Q＝[Q₁Q₂]，Q₁、Q₂分别为m×p和m×(m-o)阶矩阵，m为地表沉降数据个数，o为Φ的秩，R为上三角矩阵；

令线性参数β＝Φ⁺y，其中Φ⁺是Φ的广义逆矩阵；

结合分解式和β＝Φ⁺y得到非线性函数矩阵Φ基于QR分解的变量投影函数

6.根据权利要求5所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S2.3中，线性参数估值

其中

为将非线性参数估值

代入后所得非线性函数矩阵的广义逆矩阵。

7.根据权利要求1或2所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S3.1中，Logistic模型为

Logistic模型残差函数为r′＝y-s；其中，y为地表沉降值，s为包含时间变量t和参数A₁,A₂,A₃的地表沉降估值，r′为Logistic模型残差。

8.根据权利要求7所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S3.2中，利用LM迭代求得Logistic模型残差函数中的参数估值的方法为：

S3.2.1对Logistic模型残差函数求解关于参数A的一阶偏导，得到第二雅可比矩阵J′，所述参数A为包含参数A₁,A₂和A₃的参数；

S3.2.2给定参数初值A₀和第二阻尼因子调整系数τ′，设置收敛准则：第二迭代终止阈值ε′和第二最大迭代次数k′_max；

S3.2.3将A₀代入第二雅可比矩阵J′和r′，进而得到残差平方和

令第二阻尼因子μ′＝τ′·max{diag(J′^TJ′)}，计算第二迭代步长h′_LM＝-(J′^TJ′+μ′I)^{- 1}J′^Tr′，A_k+1＝A_k+h_LM，k≥0；

S3.2.4如果h′_LM使得

则接受A_k+1并减小第二阻尼因子μ′，进入下一步；否则，增大第二阻尼因子μ′并重新计算h′_LM，直至

S3.2.5若

或k＞k′_max，则终止迭代，输出此时的A_k+1，即为参数A的估值

否则，重复执行步骤S3.2.3和S3.2.4。

9.根据权利要求2所述的一种基于QR分解的变量投影法在地表形变中的应用，其特征在于，所述步骤S1.2中，利用以下公式对对匹配像对进行差分干涉，生成干涉像对：

δφ_j(p,q)＝φ(t_B,p,q)-φ(t_A,p,q)；

其中，t_A,t_B为某两景SAR影像的获取时间，(p,q)表示像素点位置，φ表示对应像素点的相位，δφ表示对应像素点的相位差，j为在t_A和t_B时刻获得的两景SAR影像干涉后的干涉像对序列号。

Images