[0001] 本发明涉及地表监测技术领域，具体涉及一种地表形变监测方法。

[0002] 目前，可用于地表形变监测的技术包括传统的GPS测量、水准测量、地面位移计测量、DInSAR和TS‑DInSAR技术。地表形变是一种多发的地质灾害，主要包括地面沉降(垂直向形变)、山体滑坡、火山活动位移、地震位移、矿区塌陷、冰川运动、建构筑物变形等，会对人类生命财产和经济发展造成严重的危害。对地表形变进行监测，为灾害的防治提供重要的参考和决策依据具有重要的实际意义。

[0003] 传统的GPS测量、水准测量、地面位移计测量技术工作效率低，空间分辨率低，成本较高，易受天气和地理环境限制，且难以满足对大范围区域开展高效动态监测的需求。差分合成孔径雷达干涉(Differential Interferometric Synthetic Aperture Radar,DInSAR)技术和时序差分合成孔径雷达干涉(Time  Series  Differential Interferometric Synthetic Aperture Radar,TS‑DInSAR)技术均属于对地表形变进行监测的计算。DInSAR技术具有全天时和全天候监测、工作效率高、监测范围广、空间分辨率高、成本低，不易受天气和地理环境限制的技术优势，但易受时间失相干、空间失相干、轨道误差、地形误差及大气延迟的影响，观测数据的质量受到制约，导致形变监测精度不高。TS‑DInSAR是在DInSAR技术基础上发展而来，利用覆盖同一地区的多幅合成孔径雷达(Synthetic Aperture Radar,SAR)影像进行相干点目标探测，并针对相干点上的差分干涉相位进行地表形变和高程误差建模，并利用时空滤波技术去除噪声、轨道误差和大气延迟的影响，形变监测精度最高可达毫米级。目前，TS‑DInSAR在地面沉降、山体滑坡、火山位移、矿区塌陷、冰川运动等监测中得到了较多的应用。

[0004] 但是，常规的TS‑DInSAR技术仅能监测地表沿雷达视线(Line Of Sight,LOS)方向的一维形变信息，而实际地表形变是同时包含了垂直向和水平向(水平向分为东西向和南北向)的三维形变，这一缺陷导致该技术无法获取全面的地表形变信息，严重制约了TS‑DInSAR技术的形变监测效果和进一步应用，尤其是在实际工程中的应用更是难以普及。例如，在进行地面沉降监测时，由于无法排除水平向形变而降低了沉降监测精度；在监滑坡形变时，由于滑坡位移是空间三维形变研坡向和垂直坡向的分了，TS‑DInSAR仅能获取空间三维形变在LOS向的一维分量，难以获取滑坡的准确形变信息。需要指出的是，大量研究表明，差分合成孔径雷达干涉技术对南北向形变不敏感，其所观测到的形变中南北向的贡献非常小，在监测沉降等信息时可忽略不计。

[0005] 针对这一问题，国内外学者提出了多维度小基线干涉(Multidimensional SBAS,MSBAS)技术和时序差分雷达干涉多约束建模与地表形变信息提取(Multi‑Constrained Time Series Interferometric Synthetic Aperture Radar,MCTS‑InSAR)技术，结合双源(如升轨SAR和降轨SAR)数据解算垂直向和东西向地表形变。虽然由于差分合成孔径雷达干涉技术对南北向形变不敏感，无法解算南北向地表形变，但能够区分垂直向和东西向形变，实现二维形变监测，可以进一步提高地表形变监测的精度和可靠性。

[0006] MSBAS技术是采用正则化约束的方式对二维形变解算模型附加约束条件，解决方程组秩亏的问题，让二维形变解算方程组可解，本质上市叠加单位矩阵方式减轻观测方程的病态性并使其满秩。但在该技术中正则化因子的选取难以把握，因子过大时会对时序二维形变结果造成过度平滑，使得解算结果与真实形变相差较大，因子过小时无法解决观测方程的病态问题。再者，正则化方法未考虑形变在时间维的演变趋势和物理机制，在形变较为复杂的区域监测精度受到影响。此外，当SAR影像时间采样率较低，尤其是某个时间段缺少SAR影像时，正则化方法的解算结果误差会增大，且无法恢复“数据空白”阶段的形变趋势，造成解算结果不可靠。

[0007] 为克服MSBAS技术的不足，有研究提出了采用单一的多项式形变模型对时序二维形变解算方程进行约束的方法，即上述的MCTS‑InSAR技术，在正则化的基础上增加形变约束条件，在一定程度上减轻了观测方程的秩亏和病态问题，同时可克服不同数据源联合解算过程中形变在时间维的传递问题。然而，该技术采用垂直向和东西向联合构建LOS向形变模型作为约束，使得观测方程具有较大的线性相关性，且当数据采样率不足时，观测矩阵依然存在秩亏和病态问题，因此，只能采用SVD分解法或者叠加正则化处理进行时序二维形变解算，得到的SVD最小范数解或正则化平滑后的结果与真实形变存在一定的差异，无法准确地进行时序二维形变信息获取，影响应用效果。并且，这种技术在垂直向和水平向采用统一的形变约束模型，未考虑实际中垂直向和水平向形变趋势及机制的差异性，导致最终的求解结果失真，尤其是当垂直向和水平向形变均较为显著时，该问题更为严重。

[0051] 以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

[0052] 实施例

[0053] 本发明解决上述技术问题的技术方案如下：

[0054] 本发明提供一种地表形变监测方法，如图1所示，所述地表形变监测方法包括：

[0055] S1：分别构建地表的垂直向综合形变模型和东西向综合形变模型；

[0056] S2：根据所述垂直向综合形变模型和东西向综合形变模型，构建二维形变解算模型；

[0057] S3：根据所述二维形变解算模型，得到垂直向形变时间序列和东西向形变时间序列；

[0058] S4：根据所述垂直向形变时间序列和东西向形变时间序列，得到垂直向年均形变速率和东西向年均形变速率。

[0059] 本发明具有以下有益效果：

[0060] 1、通过在垂直向和东西向分别构建约束模型，解决了观测方程秩亏的问题，直接采用最小二乘法即可进行时序二维形变模型解算，解算结果更加接近真实形变；

[0061] 2、将垂直向和东西向独立建模进行约束，这促使二维形变模型求解时的线性相关性降低，克服了原有方法存在的线性相关而导致的秩亏问题；

[0062] 3、在垂直向和东西向分别构建不同的、独立的综合性约束模型，更加符合垂直向和东西向形变机制与趋势存在差异的实际情况，可提高时序二维形变解算精度与可靠性，并且，可以分别实现垂直向和东西向形变在时间维的正确传递。

[0063] 由于不同地面形变中具体的垂直向与东西向形变组合类型以及各形变占比不同，不能确定具体的多项式形变组合模型，因此，在构建垂直向综合形变模型和东西向综合形变模型时，综合考虑地面形变的物理机制，即垂直向和东西向形变时间序列中同时包含线性形变、多阶非线性形变、周期性形变，构建适应性更强的综合性模型。由于垂直向和东西向形变模型在结构上相同，但二者之间相互独立互不影响，具有两组独立的模型参数，因此在时序二维形变解算过程中，根据真实的地面形变情况可解算出不同的形变模型参数a0～a5和b0～b5。此时，a0～a5和b0～b5分别表征了垂直向和东西向形变的不同机制和趋势。但如果在垂直向和东西向采用同一个模型(同一组模型参数)，则会强制将两个方向的形变同化，造成形变结果的失真。

[0064] 可选择地，所述步骤S1中，所述垂直向综合形变模型为：

[0065] δdi,U＝a0+a1(ti‑t0)+a2(ti‑t0)2+a3sin(2πti/T)+a4cos(2πti/T)+a5Bi[0066] 所述东西向综合形变模型为：

[0067] δdi,E＝b0+b1(ti‑t0)+b2(ti‑t0)2+b3sin(2πti/T)+b4cos(2πti/T)+b5Bi[0068] 其中，i表示第i个成像时刻；δdi,U表示第i个时刻的垂直向形变量；δdi,E表示第i个时刻的东西向形变量；a0至a5表示垂直向多项式形变模型参数；b0至b5表示东西向多项式形变模型参数。ti和t0表示形变时间序列第i个时刻和起始时刻；T表示正弦和余弦周期函数的周期；Bi表示第i个时刻SAR影像与起始时刻SAR影像间的空间基线。

[0069] 具体地，在具体实施中，假设一共有N副升、降轨SAR影像，一共形成S组干涉对以及K个相邻影像时间间隔(第一幅影像为起始时刻)，将升降轨形变时间序列根据成像时刻的先后顺序进行时序融合，以相邻两时刻间的垂直向与东西向形变增量作为待估参数建立时序二维形变解算方程，并将上述多项式模型作为约束引入解算方程中进行时序二维形变解算。

[0070] 因此，可选择地，所述步骤S2中，所述二维形变解算模型为：

[0071]

[0072] 其中，S表示S组干涉对，K表示K个相邻影像时间间隔，XU为垂直向多项式模型参数T列向量，且XU＝[a0,a1,a2,a3,a4,a5] ；XE为东西向多项式模型参数列向量，且XE＝[b0,b1,b2,T
b3,b4,b5]；GU表示Dmat和RU矩阵对应位置元素相乘的S×K阶矩阵；GE表示Dmat和RE矩阵对应位置元素相乘的S×K阶矩阵；Dmat是由0和1组成的S×K阶设计矩阵，1对应于干涉对时间跨度内包含的相邻影像间的时间间隔，其余元素均为0；RU和RE分别表示垂直向和东西向在雷达视线方向上投影系数组成的S×K阶矩阵，RU＝cosθk， k为升轨或降轨，θk为升降轨的入射角， 为升降轨的轨道方位角；ΔdU和ΔdE分别表示相邻两时刻间的垂直向和东西向形变量组成的K阶列向量，也是需要求解的参数；ΔdLOS表示S阶雷达视线方向上地表升轨和降轨形变时间序列组成的列向量，是已知的观测值；LU为垂直向形变系数矩阵；LE为东西向形变系数矩阵，且

[0073]

[0074] fU为垂直向(K+1)×6阶多项式形变模型矩阵；fE为东西向(K+1)×6阶多项式形变模型矩阵，并且

[0075]

[0076] t0～tK表示形变时间序列的起始时刻至第K时刻，T表示正弦和余弦周期函数的周期；B1～BK表示第1个时刻SAR影像至第K个时刻SAR影像与起始时刻SAR影像间的空间基线。

[0077] 在本发明所提供的具体实施例中，由于在形变监测过程中需要用到一些数据，例如升轨影像数据和降轨影像数据，因此在所述步骤S1之前，所述地表形变监测方法还包括：

[0078] A1：获取地表升轨影像数据和地表降轨影像数据；

[0079] A2：根据所述地表升轨影像数据和所述地表降轨影像数据，分别得到雷达视线方向上地表升轨形变时间序列初始向量和地表降轨形变时间序列初始向量；

[0080] A3：对所述雷达视线方向上地表升轨形变时间序列初始向量和地表降轨形变时间序列初始向量进行基准统一处理，得到雷达视线方向上地表升轨形变时间序列最终向量和地表降轨形变时间序列最终向量。

[0081] 可选择地，所述地表升轨影像数据和所述地表降轨影像数据均为SAR影像数据。

[0082] 可选择地，所述步骤S3中，将所述地表升轨形变时间序列最终向量和地表降轨形变时间序列最终向量输入所述二维形变模型中，结合最小二乘解算模型计算，分别得到垂直向形变时间序列和东西向形变时间序列。

[0083] 可选择地，所述最小二乘解算模型为：

[0084] X＝(GTG)‑1GTΔd

[0085] 其中，X为要求解的垂直向和东西向形变时间序列以及垂直向和东西向模型参数，且 G为所述二维形变解算模型的系数矩阵，且 Δd为所述二维形变解算模型的观测值矩阵，且 GU表示Dmat和RU矩阵对应位置元素相乘的S×K阶矩阵；GE表示Dmat和RE矩阵对应位置元素相乘的S×K阶矩阵；Dmat是由0和1组成的S×K阶设计矩阵，1对应于干涉对时间跨度内包含的相邻影像间的时间间隔，其余元素均为0；RU和RE分别表示垂直向和东西向在雷达视线方向上投影系数组成的S×K阶矩阵，RU＝cosθk， k为升轨或降轨，θk为升降轨的入射角， 为升降轨的轨道方位角；ΔdU和ΔdE分别表示相邻两时刻间的垂直向和东西向形变量组成的K阶列向量，也是需要求解的参数；ΔdLOS表示S阶雷达视线方向上地表升轨和降轨形变时间序列组成的列向量，是已知的观测值；LU为垂直向形变系数矩阵；LE为东西向形变系数矩阵，且[0086]

[0087] fU为垂直向(K+1)×6阶多项式形变模型矩阵；fE为东西向(K+1)×6阶多项式形变模型矩阵，并且

[0088]

[0089] t0～tK表示形变时间序列的起始时刻至第K时刻，T表示正弦和余弦周期函数的周期；B1～BK表示第1个时刻SAR影像至第K个时刻SAR影像与起始时刻SAR影像间的空间基线。

[0090] 具体地，为了便于阐述，将二维形变模型简化为：

[0091] G·X＝Δd

[0092] 其中，

[0093]

[0094] 最小二乘法的解算为：

[0095] X＝(GTG)‑1GTΔd

[0096] X为要求解的垂直向和东西向形变时间序列以及垂直向和东西向模型参数，且[0097] 由此便可解算出垂直向和东西向形变时间序列。

[0098] 可选择地，所述步骤S4中，将所述垂直向形变时间序列和所述东西向形变时间序列输入新的二维形变解算模型中，分别得到垂直向年均形变速率和东西向年均形变速率。

[0099] 可选择地，所述新的最小二乘解算模型为：

[0100] v＝(BTB)‑1·BTDT

[0101] 其中，v为年均形变速率，B为影像的时间间隔集，且 B为B的转置，D为形变时间序列集，且 d1至dK为形变时间序列，dK表示第K景影像所对应时间的垂直向形变量，ΔtK表示第K景影像到第一景影像的时间间隔。

[0102] 具体的，得到垂直向和东西向二维形变时间序列后，可分别以二者为观测值，通过最小二乘方法解算得到垂直向和东西向年形变速率。要解算形变速率，首先要建立形变速率和形变时间序列之间的关系式，以垂直向形变速率的解算为例，其观测模型为：

[0103]

[0104] 其中ΔtK表示第K景影像到第一景影像的时间间隔(单位为“天”)，v表示垂直向形变速率，是待求解的量，dK表示第K景影像所对应时间的垂直向形变量，d1至dK就是通过二维形变解算模型的简化模型和最小二乘法解算模型求解得到的垂直向形变时间序列， 为对应的改正数。将上述观测模型转换成矩阵形式为：

[0105]

[0106] 进一步地，简化表示为：

[0107] Vd＝B·v‑D

[0108] 其中，

[0109]

[0110] 那么新的最小二乘解算模型为：

[0111] v＝(BTB)‑1·BTD

[0112] 其中，v为年均形变速率，B为影像的时间间隔集，且 BT为B的转置，D为形变时间序列集，且 d1至dK为形变时间序列，dK表示第K景影像所对应时间的垂直向形变量，ΔtK表示第K景影像到第一景影像的时间间隔。

[0113] 由此便可解算出垂直向年均形变速率和东西向年均形变速率。

[0114] 可选择地，所述步骤A2中，采用常规时序干涉叠加方法，即常规TS‑DInSAR方法，得到雷达视线方向上地表升轨形变时间序列初始向量和地表降轨形变时间序列初始向量。

[0115] 在本发明所提供的地表形变监测方法中，首先需要获取覆盖同一地表区域的升轨SAR影像和降轨SAR影像。然后，分别对升轨SAR影像和降轨SAR影像进行常规TS‑DInSAR处理。以升轨SAR影像的处理为例，常规TS‑DInSAR以升轨SAR影像为基础，首先进行时序差分干涉处理，通过引用外部的数字表面模型(DSM)数据去除地形相位，得到雷达视线方向上的差分干涉相位，形成差分干涉图。然后利用SAR影像的反射强度信息探测监测区域内的相干目标点(稳定的地物点)，并提取相干目标点上的差分干涉相位。

[0116] 在得到相干目标点及其差分干涉相位后，TS‑DInSAR方法首先对相邻相干目标点之间的差分干涉相位进行求差，得到邻域差分相位，并基于此建立相邻相干目标点之间的雷达视线方向上形变速率增量(相邻两个相干目标点之间的形变差)解算模型，该模型为线性形变速率模型。基于此模型可得到相邻相干点目标之间的雷达视线方向上线性形变速率增量。

[0117] 在得到所有相邻相干点目标之间的雷达视线方向上线性形变速率增量后，可选择其中一个稳定的相干点目标为参考点，以参考点为起始点，逐次把相邻相干点目标间的线性形变速率增量进行累加，即可得到每个相干点目标雷达视线方向上的线性形变速率。用线性形变速率乘以每个TS‑DInSAR差分干涉图的时间间隔可得到与每个差分干涉图相对应的线性形变量，从差分干涉图中扣除线性形变量，然后再从剩余的相位中提取非线性形变量，最后将非线性形变量与线性形变量相加即可得到雷达视线方向上的地表形变时间序列。

[0118] 分别对升轨和降轨SAR影像进行上述TS‑DInSAR数据处理，可分别得到雷达视线方向上升轨地表形变时间序列和降轨地表形变时间序列。然后，以升轨结果中的参考点为准，在降轨结果中找到同一个点，用降轨结果中其它所有点的形变减去这个点的形变，完成升轨和降轨形变的基准统一。至此，就完成了时序二维形变建模和解算的前期准备工作，获取了时序二维形变解算模型方程式所需的观测值(已知量)，也就是经过统一基准后的雷达视线方向上升轨地表形变时间序列和降轨地表形变时间序列。

[0119] 除此之外，在本发明所提供的具体实施例中，以天津市区境内的GPS站点观测数据为例，通过将采用本发明所计算的形变监测结果与GPS形变监测结果进行对比分析，以展示本发明在时序二维形变监测中的效果和精度。

[0120] 本案例中获取了天津境内TJWQ和TJBD两个GPS站点2010至2020年、2012至2020年垂直向与东西向的精密形变时间序列数据。TJWQ点垂直向和东西向形变时间序列主要是线性形变，包含少量的周期形变和多阶非线性形变；TJBD点垂直向以周期性形变和多阶非线性形变为主，东西向以线性形变为主，包含少量其它非线性形变。以TJBD站点形变时间序列为基础模拟得到的两组形变时间序列组合(分别命名为SP1和SP2)。其中，SP1垂直向主要是线性形变，并有少量其它非线性形变，东西向主要是较小的周期性形变和多阶非线性形变；SP2垂直向和东西向均为较小的周期性形变和多阶非线性形变。

[0121] 将GPS垂直向和东西向形变时间序列分别根据Sentinel‑1A/B卫星升降轨轨道参数和周期合成两组LOS向形变时间序列。

[0122] 对于东西向形变时间序列解算结果，特别是对TJWQ和TJBD东西向较大线性形变时间序列的解算，原有方法解算结果与GPS真实形变时间序列出现了较大偏差，而本发明解算结果与GPS形变一致，明显优于原方法的效果。由此，可以确定的是，本发明能克服单一LOS向多项式形变模型约束时序二维形变解算方法的不足，同时，本发明对垂直向和东西向具有相同或不同形变时间序列规律的时序二维形变监测均具有良好的效果。

[0123] 进一步对本发明的可靠性进行定量分析，以相邻时间间隔内形变增量为基础，计算形变增量与GPS真实形变增量之间差异的均值和标准差，结果如表1所示。分析不同形变组合的解算精度，本发明技术对不同的垂直向和东西向形变时间序列解算精度基本相当，具有较好的解算可靠性与稳定性，与原有技术相比，改善效果明显。

[0124] 表1

[0125]

[0126] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。