技术领域
[0001] 本发明属于InSAR地表形变监测领域，具体涉及一种InSAR地表形变监测方法。
相关背景技术
[0002] 受SAR卫星在轨运行寿命、飞行计划以及重访周期等因素的影响，利用单一SAR卫星的SAR影像难以满足长时序InSAR形变监测对各SAR影像成像时间分布密度以及时间长度的要求，不同卫星SAR数据联合解算进行地表形变时序监测的方法应运而生。不同SAR卫星数据联合解算进行地表形变时序监测技术是指：以已配准至统一地理坐标系下的不同SAR卫星数据组合生成的干涉对的解缠相位作为观测值，以相邻SAR影像获取时刻间的形变速率作为参数，建立观测方程，然后进行参数估计获取每个SAR获取时刻的形变值。由于不同SAR 卫星成像几何参数的差异，各类SAR数据生成的干涉图属于不同的数据子集，子集间缺少形变基准，导致观测方程的系数矩阵是秩亏的，无法采用传统的最小二乘方法进行参数估计，特别是利用多源SAR数据进行多维形变时序分析必须考虑该问题。
[0003] 对上述问题的解决方法目前主要有奇异值分解和基于线性形变假设的内插或外推法。奇异值分解方法仅能获取具有数学意义的时间序列，整体的形变时间序列不正确，多类SAR数据的形变序列间存在系统偏移；基于线性形变的内插或外推方法过于依赖给定的线性形变假设，使用范围有限，若实际形变表现出明显的非线性特征，其结果也是有偏的。随着研究的深入，一些学者针对干涉组合不连续问题，提出了相应的解决办法，如利用约束条件连接间断数据子集的方法，但是由于其给定的约束条件过于简单，未考虑研究区的非线性形变特征，仅在短时间序列或简单形变地区时能够获取较好的时间序列结果；采用Tikhonov正则化约束来解决方程的病态问题，通过平滑因子来解决方程的秩亏问题。由于Tikhonov正则化理论是从优化角度和积分方程角度解决方程欠定问题，在时间序列形变求解过程中并未考虑研究区地表形变的实际变化特征，求得的结果也仅是数学意义上的最优解，不具有物理意义；在独立数据集间添加不满足短基线干涉对选取条件的额外干涉对参与建模，从而解决系统方程奇异问题，并对各个干涉对进行加权，控制各干涉对在时序形变解算中的贡献，但是解算结果依然容易被附加干涉对的质量影响，无法保证可靠性。
[0004] 综上，现有技术在解决子集间缺少形变基准问题时，存在获得的形变时间序列不具有物理意义，或依赖与附加的线性形变模型或附加的干涉对质量，不能获取符合地表形变规律的最优解的问题。
具体实施方式
[0036] 如图1所示，本发明的InSAR地表形变监测方法，具体包括以下步骤：
[0037] 步骤1：获得形变时间序列：对覆盖研究区域内所有SAR卫星采集到的SAR 影像进行处理，得到每个SAR数据子集对应的形变时间序列；
[0038] 步骤2：获取形变周期:根据步骤1中每个SAR数据子集对应的形变时间序列，计算得到对应的残余的形变时间序列；根据残余形变时间序列获得每个SAR数据子集残余时间序列的周期；并根据覆盖研究区域内所有SAR影像的时间顺序排列的时间轴，选取其中一个周期作为约束条件；
[0039] 步骤3：计算附加约束条件的形变时间序列：根据步骤2得到的约束条件，采用最小二乘法解算观测方程，获得非线性地表形变时间序列，并将其加上对应的线性形变，得到最终地表形变时间序列。
[0040] 上述技术方案通过选取具有物理意义的形变周期作为约束条件从而获取更具物理意义的形变时间序列；同时，该方法不依赖于附加的线性形变模型或附加的干涉对的质量，能够获取符合地表形变规律的最优解，从而实现研究区域内地表形变的高精度、长时间跨度的时序监测。
[0041] 具体的，步骤1中可采用PS-InSAR方法或者SBAS-InSAR方法对SAR影像进行处理，采用哪种方法主要依据研究区域内SAR数据影像数量的多少，能够准确得到每个SAR数据子集对应的形变时间序列。
[0042] 优选的，步骤2具体包括如下操作：
[0043] 2.1将步骤1获取的不同SAR数据子集的形变时间序列配准至统一的地理坐标系下，并以每个SAR数据子集下的解缠相位作为观测值，计算研究区域的线性形变速率；
[0044] 2.2从步骤1得到的每个SAR数据子集对应的形变时间序列中扣除该形变时间序列对应的线性形变，得到残余的形变时间序列；对残余形变时间序列采用快速傅里叶变换(FFT)进行频谱分析，获取每个SAR数据子集的残余时间序列的周期Tk,j(j＝1,2,…,L)；L为每一类SAR卫星的残余时间序列的周期的个数；
[0045] 2.3将覆盖研究区域内的各类SAR卫星下的所有SAR影像按时间顺序排列在同一时间轴上；从步骤2.2中得到的每个SAR数据子集的残余时间序列的周期中，选取周期T作为约束条件，所述周期T不小于所述时间轴上来自不同SAR 卫星的SAR影像之间的时间间隔。
[0046] 步骤2的操作方式统一了子集间的形变基准，获得了具有物理意义的形变周期作为约束条件。
[0047] 优选的，所述步骤2.1中，采用公式1计算研究区域的线性形变速率：
[0048]
[0049] 式中， 表示第k类SAR卫星由时刻tj与ti构成的干涉组合对应的解缠相位，k为研究区域内一类SAR卫星的序号；λk表示第k类SAR卫星的雷达波长； v表示线性形变速率；表示第k类SAR卫星由时刻tj与ti构成的干涉组合对应的残余相位。
[0050] 优选的，所述步骤2.2中，所述形变时间序列对应的线性形变的计算公式为：
[0051] v(ti-t0)；
[0052] 其中t0为步骤1中采集SAR影像的起始时刻，ti为步骤1中采集SAR影像的某一时刻；v表示线性形变速率。
[0053] 优选的，步骤3具体包括如下操作：
[0054] 3.1对公式2所示的观测方程和公式3所示的约束条件方程采用最小二乘法解算，获取非线性地表形变时间序列；
[0055]
[0056]
[0057] 式中，vs和vr表示相邻两景SAR影像间的形变速率，ts+1、ts以及tr+1和tr均表示相邻两景SAR影像的获取时刻；tm+n-tm＝aT,a为正整数； 表示第k类 SAR卫星由时刻tj与ti构成的干涉组合对应的残余相位。
[0058] 显然，当约束条件方程个数大于独立子集的个数时，公式(2)、(3)组成的方程组存在唯一解。因此，加入约束条件方程后的观测方程系数矩阵可达到满秩条件。
[0059] 3.2将所述非线性地表形变时间序列加上步骤2.3中的时间轴上时刻tc对应的线性形变v(tc-t0')，其中，t0'为所述时间轴上的起始时刻，tc为时间轴上任一时刻，得到最终的地表形变时间序列。
[0060] 步骤3的操作方式克服了观测方程的系数矩阵是秩亏的，无法采用最小二乘方法进行参数估计的问题，最终得到地表形变时间序列的最优解。
[0061] 实施例：
[0062] 由于采用多个SAR卫星得到多个SAR数据子集与采用一个SAR卫星得到多个数据在理论上的处理效果是类似的，考虑到多个卫星的数据完整性不能兼顾，本发明具体实施方案以美国南加州地区且仅考虑单一SAR卫星的多个SAR 数据子集为例来详细描述本发明的效果。
[0063] 首先对覆盖研究区域内采集到的SAR影像进行处理，得到每个SAR数据子集对应的形变时间序列，将所述每个SAR数据子集对应的形变时间序列配准至统一的地理坐标下，得到研究区域的线性形变速率；从得到的每个SAR数据子集对应的形变时间序列中扣除该线性形变时间序列对应的线性形变，得到残余形变时间序列；然后对每个SAR数据子集的残余形变时间序列进行FFT频谱分析获取对应的时间序列周期(如2所示)；最后选取周期作为约束条件，结合观测方程与约束条件方程，采用最小二乘法进行解算，得到非线性地表形变时间序列，加以处理最终获得地表形变时间序列。
[0064] 为了验证本发明的应用效果，选取BLSA、LBC1、PMHS、SACY 4个典型GPS (全球定位系统)点，采用三种方案进行时间序列形变的恢复，并与GPS监测结果进行了对比分析。方案一为常规的SBAS(短基线集技术)，方案二为基于线性内插或外推的最小二乘算法，方案三为附加约束条件的秩亏InSAR时序解算方法。如图3所示，相比于方案一和方案二，方案三获取的时间序列与 GPS时间序列之间的吻合程度更高，方案三获取的时间序列与GPS时间序列之间差值的标准差分别为5.0mm/a(BLSA)、3.6mm/a(LBC1)、6.3mm/a(PMHS)、 5.9mm/a(SACY)，精度均为毫米级，有效证明了本方法的可行性与有效性。
