Sentinel-1A数据矿区地表形变监测适用性分析

白泽朝， 汪宝存， 靳国旺， 徐青， 张红敏， 刘辉

(1.信息工程大学地理空间信息学院，郑州 450001； 2.河南省地质矿产勘查开发局测绘地理信息院，郑州 450006； 3.华北水利水电大学测绘与地理信息学院，郑州 450046)

0 引言

煤矿开采导致的地面沉降，严重影响人类的生命财产安全。合成孔径雷达干涉测量(interferometric synthetic aperture Radar，InSAR)技术作为一项重要的监测手段[1-2]，在矿区地表形变监测中的应用越来越得到广泛认可，国内学者也对此进行了大量的实验研究。吴立新等[3]利用合成孔径雷达差分干涉测量(differential InSAR，DInSAR)技术获取了唐山市及开滦矿区形变图，分析了试验区沉陷的扩展及演变过程，并对DInSAR技术应用中存在的时间去相干、空间去相干等误差因素进行了分析和讨论； 董玉森等[4]采用JERS-1卫星数据通过差分处理获取6景差分干涉图，发现研究区有4处沉降区域； 朱建军等[5]介绍了InSAR技术在矿区地表形变监测中的应用现状及进展； 汪宝存等[6]以永城市为例，实验表明联合利用DInSAR和小基线集2种技术适合在矿区开展地表形变监测； 刘广等[7]利用重轨差分InSAR技术获得了峰峰矿区地表ENVISAT和JERS1的雷达形变结果，对C波段和L波段数据相干特性、相位特性以及干涉测量技术在矿区地表沉降监测中应用的可行性和局限性进行了研究。相关成果的取得主要集中在冬春季节，这是因为在这个季节植被干扰小，SAR数据保持着较好的相干性。但是在植被生长旺盛季节，由于失相干严重，导致利用InSAR技术开展矿区地表形变识别与监测难以顺利开展，因此在相当程度上阻碍了矿区地表形变InSAR监测技术的推广和应用。

此外，相关研究表明[8-10]，时间基线是影响干涉对相干性的重要因素之一。这一因素在不同地物类型表现不同，在植被区，时间去相关现象更为复杂，尤其是在夏季植被覆盖茂密时间段内相干性衰减非常迅速。由于卫星技术的迅速发展，2014年发射的Sentinel-1A卫星重访周期为12 d，如果与2016年发射的Sentinel-1B卫星进行组网观测，重访周期只有6 d。重访周期的缩短将有效提高SAR数据的相干性，削弱时间失相干在夏季植被覆盖区的影响。

针对InSAR监测矿区地表形变中形变位置确定、形变梯度估计和相干性之间的关系，为客观评价Sentinel-1A数据在矿区形变应用中的监测能力，需要尽可能保证1 a的监测周期，对于不同地表类型分别进行量化分析。本文选择矿区3个典型形变区域，引入形变梯度模型，通过对研究区不同地表类型进行相干性的定量化分析，以及采用目视解译的方法检测形变条纹的可见度，研究在1 a的监测周期内，矿区形变监测中Sentinel-1A数据相干性和形变位置确定、形变梯度估计之间的关系，旨在说明Sentinel-1A数据在矿区形变应用中的监测能力，分析Sentinel-1A数据能否实现全年的形变监测。

1 形变模型和数据获取

1.1 形变模型

引入InSAR可检测的最大最小形变梯度函数模型[11]，验证针对Sentinel-1A数据模型的适用性，以及定量化分析不同季节、地表类型和时间基线影响下相干性与可检测形变之间的关系。首先需要考虑InSAR技术本身的约束，当相干性γ=1，即没有失相干时，InSAR所能检测到的最大形变梯度为相邻像元不超过一个条纹[11]，即

(1)

式中：D为形变梯度；μmin为不同视数下像元最短边长，本文采用5视，μmin的取值为20 m，则最大的形变梯度为1.4×10-3。此外，研究表明[12-13]，相干性γ=1时，可监测到的最小形变梯度约为10-7，相当于250 km幅宽的Sentinel-1A影像内发生了25 cm的形变。这2个约束条件表明形变梯度模型必须经过(1，D)和(1，10-7)，因此线性函数模型可以定义为

Dmax(γ)=D+Kmax(γ-1)，

(2)

Dmin(γ)=10-7+Kmin(γ-1)，

(3)

式中：Dmax和Dmin分别为最大和最小形变梯度；Kmax和Kmin分别为最大和最小形变梯度斜率。本文采用Jiang等[12]建立的5视下形变梯度与相干性的经验函数模型

Dmax(γ)=0.002 0(γ-1)+0.001 4，

(4)

Dmin(γ)=-0.000 081(γ-1)+10-7。

(5)

对上述模型求取相干性阈值，令Dmax=Dmin，可以得到γ=0.32，即当相干性γ<0.32时，差分干涉图无法反映出任何地表形变信息，整体相干性优于0.32，微小的形变就有可能被恢复。

1.2 数据获取

实验数据采用2015年6月—2016年5月1 a间的Sentinel-1A卫星降轨获取的21景干涉宽带模式SAR影像，空间分辨率为5 m×20 m。数字高程模型为SRTM的3″数据，格网间隔为90 m×90 m。精密轨道数据采用成像21 d之后发布的精密轨道数据，定位精度可达5 cm。需要对相邻影像进行差分干涉处理，影像数据干涉对分布如图1所示。

图1 相邻影像干涉对分布

图1中150526表示2015年5月26日。为进行形变梯度函数模型验证分析，获取了SAR图像覆盖范围内某一矿区沉降中心位置和沉降区域外围2个水准点的观测成果，水准观测数据采用二等水准测量作业方式获取，作业时间为2015年7月31日—2016年7月29日，期间共监测9次，水准监测时间同SAR影像获取时间吻合较好。通过2个水准点之间做差获取矿区实际沉降量值。图2为沉降量随时间变化情况，采用二次多项式拟合，公式为

y=-0.009 3-0.005 6x+0.000 009x2。

(6)

图2 实际沉降量随时间变化情况

2 InSAR矿区形变监测定性分析

选取河南省焦作市某矿区为研究区，在实地调查的基础上，根据矿区地表覆盖情况，选择3个研究区域，分别为①伴有杂草的裸地(A)； ②农田(B)，春季种植小麦，夏季种植玉米； ③村落(C)，周边种植玉米。底图为WorldView-3光学影像。

(a) 研究区范围 (b) 区域A

(c) 区域B (d) 区域C

获取裸地、农田和村落3个区域相邻影像20组时间序列差分干涉图如图4—6所示，每组干涉对的成像时间、时间基线Bt以及平均相干性γ已分别在图上标出。

图4-1 区域A裸地时间序列差分干涉图

图4-2 区域A裸地时间序列差分干涉图

图5 区域B农田覆盖区时间序列差分干涉图

图6-1 区域C村落时间序列差分干涉图

图6-2 区域C村落时间序列差分干涉图

从图中可以较好地说明Sentinel-1A数据在1 a间相干性的变化情况，其中差分干涉对中最大时间基线为60 d，最小时间基线为12 d，大部分时间基线为12 d和24 d。此外空间基线对相干性存在一定的影响，干涉处理中视数对形变梯度模型的影响需要单独说明。对于空间基线的影响，由于Sentinel-1A卫星轨道精度较高且采用先进的轨道控制技术，干涉对空间基线最大为162 m，多数都在100 m范围内，故忽略干涉对空间基线的差异。对于视数的影响，由于Sentinel-1A数据常常采用视数为5进行处理，故不考虑其他多视情况。

图4—6中，裸地、农田和村落区域的相干性数值分别在0.52～0.94，0.42～0.93和0.51～0.95之间。可以看出，成像时间在夏季(150617—150816)时间段内有4组干涉对，其中150723—150816干涉对相干性为全年最低，农田区域最低为0.42，裸地区无明显形变，农田区和村落区均因相干性过低，无法目视辨认其形变条纹。成像时间在秋季和冬季(150816—160107)时间段内有11组干涉对，秋季裸地区从150828—150909干涉对出现较为清晰明显的形变条纹，而农田区和村落区形变条纹模糊，3类地表均可有效地目视判别出形变的位置和范围。在秋末干涉对相干性明显提高。冬季为1 a中相干性最好的时期，完全可以进行形变量的精确反演。成像时间在春季(160107—160530)时间段内有5组干涉对，60 d时间基线的160107—160307干涉对也能保持较高的相干性，轮廓清晰，裸地区和农田区形变条纹数可以有效目视辨认。在春末夏初，时间基线为24 d的160412—160506干涉对和160506—160530干涉对裸地区仍保持较好的相干性，形变条纹的轮廓也可清晰辨认； 而农田区相干性较低，形变条纹缺失； 村落区相干性较高，但由于周边是农田，形变条纹缺失。

图7为3种地表类型1 a间相干性变化情况，夏季相干性均较低，存在1 a中相干性最低的干涉对，在秋末干涉对相干性明显提高，冬季为1 a中相干性最好的时期，春季相干性同样较低。3种地表类型对比可以发现，农田区在夏季和春季相干性明显低于裸地区和村落区，冬季和秋季相干性基本一致。裸地区和村落区全年相干性相似，γ均低于0.5的有3组干涉对，时间为夏季和春季，且干涉对时间基线均为24 d。

图7 相干性变化折线图

3 InSAR矿区形变监测定量分析

对于上述目视判读分析结果可以通过形变梯度函数模型进行定量化分析，分别选取2种地表类型干涉对，由于裸地区域在夏季无明显形变，选取春季160506—160530干涉对，农田区域则选取夏季150617—150629干涉对。裸地区目标a和b相距约为160 m，形变量约为28 mm(1个条纹)，相干均值γ为0.58，计算形变梯度D为0.175×10-3。农田区目标c和d相距约800 m，形变量约为28 mm(1个条纹)，相干均值γ为0.57，计算形变梯度D为0.033×10-3。通过目视解译可知，验证区中裸地区域(图8(a))条纹轮廓清晰，条纹数较易辨认； 农田区域(图8(b))条纹轮廓清晰，但受强噪声影响条纹数不易辨认，形变信息丢失。

(a) 裸地区域 (b) 农田区域

经过目视判读和定量分析，将形变梯度模型(0.58，0.175×10-3)和真实数据(0.57，0.033×10-3)一并绘制在图9中，2条线之间为可预测范围，观察模型预测结果和上述分析吻合程度。可以看出，裸地区域(三角形点)的形变梯度可以监测到，与分析的结果吻合； 农田区域(正方形点)位于可监测的最小形变梯度线上，通过目视识别分析，可以确定形变轮廓，但形变条纹模糊无法辨认。模型估算与实际观测结果较吻合，表明该模型同样适用Sentinel-1A数据的形变估算。

图9 模型验证结果

计算水准点实际形变梯度及对应干涉对平均相干性，并将相干性和形变梯度作为模型输入值，根据目视判读每组干涉对形变梯度的可检测性与水准数据带入模型计算得到的可检测性进行对比。图10为对比结果，红色为目视判读与模型计算形变可检测性两者结果一致，黑色为目视判读与模型计算两者结果不一致。16组验证数据中，两者结果一致的有12组，一致率为75%； 剩余4组为实际水准数据带入模型，表现为可检测形变梯度，而实际情况通过目视判读形变条纹缺失，属于不可检测情况，两者结果存在矛盾。

图10 模型一致性验证结果

分析存在矛盾的4组干涉对，例如干涉对为151003—151015，通过目视判读形变条纹模糊，为不可检测情况，但带入模型为可检测形变梯度。其余3组分别为160307—160331，160331—160412和160412—160506，虽然相干性均较高，但目视判读形变条纹缺失，为不可检测情况。通过实地调查，形变区位于村落和农田之间，村落在1 a中均保持较高的相干性，使得3组干涉对具有较高的整体相干性，形变条纹在村落部分保持较清晰，而农田区域由于相干性较低，形变条纹缺失。表明对于地物类别复杂区域，尤其是具有建筑物等人工设施的区域，研究区局部相干性较高，整体平均相干性无法有效代表该区域的形变条纹清晰程度，因此，相干性的统计分布也是需要考虑的一方面。

4 结论

以暖温带半湿润季风气候条件下的河南省焦作市某矿区为研究区，本文采用Sentinel-1A数据研究形变监测中不同地表类型、季节和时间基线影响下相干性和可监测形变之间的关系，并采用水准数据验证形变梯度函数模型和目视判读分析结果，得出以下3方面结论：

1)裸地和村落类型全年保持较高相干性，夏季时间基线为24 d时，相干性均优于0.51，目视识别可以有效确定矿区的形变位置和范围，通过模型验证形变梯度也在可检测范围。

2)农田覆盖类型夏季时间基线为12 d时，相干性优于0.51，目视识别可以有效确定矿区的形变区域和范围，但受噪声影响条纹模糊，通过模型验证形变梯度位于可检测最小形变梯度上； 当时间基线为24 d时，无法有效确定形变区域和范围。

3)采用真实水准数据验证模型的适用性，16组验证数据中模型准确率为75%。由于在地物类别复杂区域，局部相干性较高，整体平均相干性无法有效代表该区域的形变条纹清晰程度。本文研究能初步帮助在InSAR应用中根据研究目标地物类型和监测季节选择合理的数据分布。