[0001] 本发明属于矿区测量技术领域，具体涉及一种基于开采沉陷先验模型的矿区地表高时空分辨率三维形变估计方法。

[0002] 高时空分辨率矿区地表三维形变监测对于理解矿区开采沉陷动态机理和矿区灾害预防和精确评估具有重要作用。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)技术由于其高空间分辨率、高精度、不受云雨天气干扰等优势在矿区地表形变监测和预计中发挥着巨大的作用。然而，InSAR技术监测得到的形变是地表真实三维形变在雷达视线向的一维投影，而一维雷达视线向形变很难反应矿区真实地表形变。针对此问题，Li等于2015年提出一种名为SIP(single InSAR pair)的方法(参见文献1Li Z,Yang Z,Zhu J,Hu J,et al.Retrieving three-dimensional displacement fields of mining areas from a single InSAR pair[J].Journal of Geodesy,2015,89(1):17–
32.)，该方法基于矿区地表形变的物理特性，利用单个InSAR干涉对反演得到矿区地表三维形变。然而，该方法只能获得组成干涉对的两景合成孔径雷达(Synthetic Aperture Radar，SAR)影像期间的矿区地表形变。为了克服该局限，Yang等于2017年提出一种SGI(singlegeometry interferometric)方法将SIP方法发展到时序上(参见文献2Yang Z,Li Z,Zhu J,Feng G,et al.Deriving time-series three-dimensional displacements of mining areas from a single-geometry InSAR dataset[J].Journal of Geodesy,2017,
92(5):529-544.)，该方法基于一个单几何的InSAR观测数据集，利用加权最小二乘的方法对多时相干涉对进行联合解算，实现了矿区地表三维时序形变获取方法。

[0003] 然而，Yang的方法存在如下明显的局限：1)对于SGI计算得到的矿区地表时序形变结果，其时间分辨率完全由所选单几何InSAR数据的影像获取时间间隔决定。目前为止，SAR卫星的重访周期通常是固定不变的，但是星载/机载SAR卫星的数量却越来越多。而对于SGI方法而言，只要所选SAR数据类型固定，则所得到的时序形变结果的时间分辨率将不会高于这一SAR传感器的重访周期。对于高度非线性的矿区地表形变而言，这将对矿区地表形变的精确预计和建模，以及对开采导致的地质灾害的精确预计和评估造成严重影响。2)由于矿区大多处于农村以及郊区，其地表覆盖物通常变化较大，加上地表形变速率较快，使得基于单轨InSAR数据很难获得矿区地表时序形变结果。尤其是当由单轨数据集中存在时间断缺时，利用传统的基于最小二乘的算法将很难获得高度非线性的矿区地表时序形变结果。

[0061] 下面将结合实施例对本发明做进一步的说明。

[0062] 如图1所示，本发明实施例提供一种对矿区地表形变能更加精确预计的方法，包括如下步骤：

[0063] S1：基于每个轨道获取的SAR影像分别组成可干涉的InSAR干涉对集合。

[0064] 其中，L个轨道的SAR影像分别组成L个SAR影像集来覆盖待监测矿区，L为大于或等于2的正整数。每个SAR影像集的数量分别表示为：M1,M2，…，ML，SAR影像的总数量为M+1，即，M+1＝M1+…Mj…+ML。

[0065] 基于每个轨道的SAR影像集分别生成每个轨道对应的InSAR干涉对集合，每个InSAR干涉对集合中InSAR干涉对的数量分别为：G1,G2,…,GL，InSAR干涉对的总数来那个为G，即，G＝G1+…Gj…+GL,1≤j≤L。

[0066] 需要说明的是，本实施例步骤S1中每个轨道对应的InSAR干涉对数据集是根据不同InSAR数据源特征和待监测矿区形变量级构建的，具体是根据不同InSAR数据源的影像质量和雷达波长等特征以及待监测矿区开采引起的地表形变速率共同确定的干涉对时空基线阈值生成的，但阈值的选取和确定具有一定的经验性，即以上述内容为参考设定的经验值。而根据干涉对时空基线阈值如何构建InSAR干涉对数据集为现有技术内容，因此本发明对此不进行具体的说明。

[0067] S2：依靠外部数字高程模型(Digital Elevation Model，DEM)去除地形相位，利用“二轨法”差分InSAR(Differential InSAR，DInSAR)技术分别处理各生成的InSAR干涉对，从而获得多轨InSAR数据联合监测得到的待监测矿区地表一维雷达视线向多时相形变监测值。

[0068] 其中，得到的一维雷达视线向多时相形变监测值的总数量为G，其中，DInSAR技术是现有成熟技术，使用DInSAR技术处理InSAR干涉对也是本学科经典数据处理技术手段之一，因此，本发明对此不进行具体说明。

[0069] 在实际应用中，由于不同轨道获取的SAR影像的成像几何参数不同，而同一轨道的SAR影像集的成像几何参数大致相同，故本发明所有SAR平台的成像几何参数可分为L组。基于得到的不同SAR成像几何的一维雷达视线向多时相形变场进行解算前，是需要对其进行坐标转换的。本实施例选择用于去除InSAR干涉对中地形相位的外部数字高程模型(Digital Elevation Model,DEM)的坐标框架几何参数为参考，并选用GAMMA软件的DIFF&GEO模块进行处理，将监测得到的各成像几何不同的InSAR一维雷达视线向多时相形变场进行地理编码和数据重采样，得到统一地理坐标框架后的一维雷达视线向多时相形变场，再执行步骤S3。其中，外部DEM的坐标系统通常为WGS 1984标准大地坐标系。

[0070] S3：基于矿区地表水平移动与下沉梯度呈线性比例关系的先验模型分别构建矿区地表沿东西、南北方向的水平移动与下沉梯度之间的函数关系，并基于所述函数关系以及一维雷达视线向多时相形变监测值与矿区地表沿着垂直、东西、南北方向的三维多时相形变观测场的关系方程分别计算出矿区地表沿着垂直、东西、南北方向的三维多时相形变观测场。

[0071] 其中，矿区地表东西、南北方向的水平移动与下沉梯度之间的函数关系如下所示：

[0072]

[0073]

[0074] κE＝b·r(i,j),κN＝b·r(i,j)(i＝1:Rmax，j＝1:Cmax，h＝1:G)

[0075]

[0076] 式中， 和 分别表示第h个多时相形变观测场中(i,j)像元位置地表沿东西、南北和垂直方向的位移，κE和κN分别为东西、南北两个方向上的水平移动与下沉梯度之间的比例系数，b和r(i,j)分别表示(i,j)位置处地表形变的水平移动系数和主要影响半径，东西、南北方向上的b和r(i,j)分别相同； 和 表示第h个多时相形变观测场中(i,j)像元位置处的垂直下沉在东西和南北方向上的沉降梯度，分别表示每个像元沿东西和南北方向的地面分辨率，Rmax和Cmax分别表示影像的最大行列数；

[0077] 一维雷达视线向多时相形变监测值与矿区地表沿着垂直、东西、南北方向的三维多时相形变观测场的关系方程如下所示：

[0078]

[0079] 且：

[0080]

[0081]

[0082]

[0083] 式中，DLOS为一维雷达视线向多时相形变观测场，θ和α分别表示各多时相InSAR干涉对对应的SAR传感器的入射角和飞行方位角；UT、ET、NT分别表示矿区地表沿着垂直、东西、南北方向的三维多时相形变观测场，每个维度的多时相形变观测场的数量均为G，分别为G个多时相形变观测场中的第h个多时相形变观测场。

[0084] 将上述矿区地表水平移动与下沉梯度的比例关系函数带入上述一维雷达视线向形变与三维形变之间的关系方程中，并在InSAR监测得到的矿区全盆地多时域形变影像中分别进行逐行/列迭代求解，即可得到矿区地表全盆地三维多时相形变监测结果。

[0085] S4：构建矿区地表三维形变速率与三维多时相形变观测场之间的观测方程组；

[0086] 其中，将每景SAR影像按照时间先后排序，排序的SAR影像的时间t为：t＝[t0,t1,…ti…,tM]，其中，ti为排序的SAR影像中第i+1景SAR影像时间；

[0087] 再基于排序的SAR影像构建干涉对的时间差系数矩阵T。干涉对的时间差系数矩阵T由每个干涉对的时间差系数向量组成，每个干涉对的时间差系数向量由M个元素组成，按序对应所述M+1景SAR影像的时间间隔。因此，本实施例中，时间差系数矩阵T的维度为G×M，每行对应为一个干涉对的时间差系数向量。其中，第k个InSAR干涉对的时间差系数向量(本实施例中第k行)中第MFk个元素前以及第MSk-1个元素以后的所有元素均为0，从第MFk到第MSk-1个元素之间的值依次为：

[0088] 其中，T(k,MFk)、T(k,MSk-1)分别表示第MFk和第MSk-1元素的值，MFk、MSk分别对应第k个InSAR干涉对的主、辅SAR影像的时间索引， 分别为生成第k个InSAR干涉对的主、辅SAR影像的时间； 分别为基于排序的SAR影像的时间t中MFk的下一时间，MSk的上一时间。

[0089] 此外，假设排序后的相邻SAR影像期间，矿区地表沿着垂直、东西和南北方向的三维形变速率为恒定值，即相邻SAR影像期间在垂直、东西和南北方向的三维形变速率用三维形变速率矩阵表示，其具体为如下：

[0090] VU＝[VU1,…VUi…,VUM]T，VE＝[VE1,…VEi…,VEM]T和VN＝[VN1,…VNi…,VNM]T；

[0091] 其中，VUi、VEi、VNi分别表示基于排序的SAR影像中第i与第i+1景SAR影像在垂直、东西、南北方向的形变速率。

[0092] 基于上述干涉对的时间差系数矩阵T以及三维形变速率矩阵的形式，构建的观测方程组如下所示：

[0093]

[0094] S5：基于步骤S3获取的矿区地表沿着垂直、东西、南北方向的三维多时相形变观测场以及基于每个SAR影像时间构建的干涉对的时间差系数矩阵采用广义最小二乘法求解步骤S4中的所述关系方程组得到三维形变速率矩阵。

[0095] 其中，计算公式如下：

[0096]

[0097] 其中， 分别为求解出的三维时序形变速率矩阵，以及垂直、东西、南北方向的形变速率矩阵，(TT·P·T)+表示(TT·P·T)的Moore-Penrose伪逆；P为参与解算的各个干涉对所占权重。

[0098] S6：利用如下公式以及步骤S5计算出的三维形变速率矩阵计算出待监测矿区在垂直、东西、南北方向的高时空分辨率三维形变；

[0099]

[0100]

[0101]

[0102] DU(t0)＝DE(t0)＝DN(t0)≡0

[0103] 其中，DU(tk)、DE(tk)、DN(tk)分别表示基于排序的SAR影像中第k+1景SAR影像获取时刻待监测矿区地表在垂直、东西和南北方向的高时空分辨率三维形变累积量；DU(t0)、DE(t0)、DN(t0)分别表示基于排序的SAR影像中第1景SAR影像获取时刻待监测矿区地表在垂直、东西和南北方向的三维形变值；

[0104] 为了更加清楚的说明本发明及其优势，下述将结合部分图来进一步解释本发明提供的所述方法。在基于矿区开采沉陷先验模型解算得到三维多时相形变观测量后，三维时序形变是基于相同的逐点解算系统完成的。鉴于本发明的特征，本节以矿区单点沉降为例，对本发明进行具体实施例演示与分析。本实施例利用logistic模型(一种矿区地表单点沉降预计数学模型)模拟了地下开采导致的地表单点沉降，如图2a中灰色曲线所示，作为本实施例的原始观测数据模型。设置ALOS PALSAR卫星的重访周期：46天作为相邻观测的时间间隔，分别模拟得到两个观测数据集(用T1和T2表示)，如图2a中实心圆和十字星号所示。以某矿区的实际监测数据的时空基线为参考，通过给定的时空基线阈值模拟生成由T1和T2组成的InSAR干涉对集合。由这些InSAR干涉对组成以及时空基线情况绘制出如图2b所示的基线集网络，其中黑色网络由T1数据集数据组成，灰色网络由T2数据集数据组成，浅灰色阴影部分表示了在T1数据集网络中存在的时间断缺部位。利用logistic模型将生成的干涉对集合内的各干涉对涵盖时间间隔内包含的沉降量依次计算出来，并分别加入均值为0，均方根误差为3cm的高斯白噪声，从而模拟得到矿区多轨道多时相沉降监测数据集。

[0105] 对于单轨数据集T1，我们利用经典最小二乘算法进行解算，并和利用本方案得到的多轨道数据融合解算结果一并与相应的模拟真值进行对比，如图3所示。我们可以发现：当单轨数据中存在时间断缺时(如图3a)，由其解算得到的时序沉降结果在时间断缺附近将存在较大误差，并在此之后的时序形变解算结果中进行累积，从而导致整个时序形变解算结果出现较大错误；与模拟真值相比差异非常大，均方根误差达到了25cm，占形变量最大值的25％，如此大的误差在矿区地表形变监测中是不可以被接受的，从而使得时序形变解算失败。然而，当加入多轨数据进行融合解算时(如图3b)，所得到的结果能有效地重建由于某单轨数据中存在的时间断缺部分，消除由于时间断缺导致的解算错误，与模拟真值相比，均方根误差也得到大幅度减小，降低至1.5cm，只占最大形变值的1.5％；且相对于单轨数据，融合更对轨道的InSAR数据将能够有效提高监测结果的时间分辨率，在本实施例中，相对于单轨数据解算结果，基于两个轨道数据集融合解算的时序形变解算结果的时间分辨率提高了大约43.7％；得到了高精度的时序形变解算结果并增强了形变结果抵抗粗差的能力。

[0106] 需要强调的是，本发明所述的实例是说明性的，而不是限定性的，因此本发明不限于具体实施方式中所述的实例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，不脱离本发明宗旨和范围的，不论是修改还是替换，同样属于本发明的保护范围。