全分辨率与多视法DInSAR矿区地面沉降监测研究

李宾， 陶秋香， 牛冲， 刘晓帅， 张亚凤

1. 山东科技大学 测绘与空间信息学院,山东 青岛 266590； 2. 山东省地质测绘院,济南 250014

地面沉降主要是一种在自然和人为因素的作用下,使地下松散土层及岩层压缩而导致地面高程逐渐降低的工程地质现象.随着社会经济飞速发展,人们在过度开采煤炭、石油、天然气、地下水等资源以满足自身需要的同时[1-2],地面沉降的程度和范围正在进一步加深和加大[3-4],这已成为经济社会可持续发展的重要制约因素.合成孔径雷达差分干涉测量(Differential Interferometric Synthetic Aperture Radar,DInSAR ) 技术具有全天候、全天时、高分辨率和监测范围广等优势,能够应用于矿区地面沉降监测[5].1989年Grabriel等[6]利用Seasat卫星L波段的SAR影像首次论证了 DInSAR可用于厘米级的地表形变监测.而后,国内外研究人员开展了DInSAR监测矿区地面沉降的相关理论与应用研究,取得了一系列研究成果[7-12].现有研究大多在差分干涉测量之前进行了多视处理,以抑制SAR影像的斑点噪声,进而获取矿区地面沉降信息,即多视法DInSAR.然而,多视处理会降低影像的分辨率,不可避免地会降低InSAR监测到的最大沉降梯度和沉降量[13-14].全分辨率DInSAR则对SAR影像不进行多视处理,在保持原有SAR影像分辨率的基础上进行差分干涉处理,从而提高DInSAR技术监测最大沉降梯度和沉降量的能力．

本研究以山东省某矿区为研究区,选取2018年10月18日-2019年1月22日时间段内6景哨兵数据为数据源,分别用全分辨率和多视法(距离向与方位向之比为5∶1)进行差分干涉处理,并结合矿区同时间段的工作面开采信息和水准数据监测,系统对比分析和研究了全分辨率和多视法DInSAR在矿区地面沉降监测中的应用效果.研究结果表明,由于受斑点噪声的影响较为严重,全分辨率DInSAR在沉降梯度和沉降量都较大的沉降区域,其差分干涉条纹的连续性变差,相位解缠和形变解算的精度降低,对比多视法DInSAR在矿区沉降盆地中心的地面沉降监测结果,并未有明显改善．

1 研究区概况

研究区位于华北平原的南缘,黄河北岸,地理坐标范围为116°59′00″-117°06′00″E,36°47′00″-36°50′00″N.图1给出了矿区的具体位置及范围.矿区东西向长约10 km,南北向宽约5 km,面积约50 km2.矿区南距济南市区8 km,北距济阳县城20 km,有国道、高速穿过井田,交通运输便利.矿区地面标高+22～+27 m,地势平坦,局部地势低洼,雨季易积水.黄河从井田南部1 km处流过,井田中东部有邢家渡引黄总干渠,中西部有大寺河通过,为本井田主要地表水系．

审号图：GS(2019)1822号．图1 研究区地理范围

2 数据源

2.1 SAR影像数据

在矿区形变监测实际应用中,夏季植被覆盖茂密,失相干严重,对于矿区形变区域的定位和形变监测反演造成很大阻碍[15-16].研究中采用冬季的SAR影像数据进行DInSAR处理,获取矿区地面沉降信息.本研究选取了覆盖矿区的时间跨度为2018年10月18日-2019年1月22日的6景Sentinel-1A升轨SAR影像,组成5组干涉对进行差分处理.SAR影像参数如表1所示．

表1 DInSAR差分干涉对

2.2 DEM数据

航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission,SRTM)由美国太空总署(National Aeronautics and Space Administration,NASA)和国防部国家测绘局(National Imagery and Mapping Agency,NIMA)以及德国与意大利航天机构共同合作完成联合测量[17-18].2000年2月11日,美国发射“奋进”号航天飞机,其载SRTM系统,进行了11 d的数据采集工作,获取了覆盖地球80%以上的陆地数字高程模型(Digital elevation model,DEM).SRTM数据是迄今为止现势性最高、分辨率最高、精度最高的全球数字地形数据.为了相位校正和地理编码,本研究采用了地面分辨率约90 m的SRTM3 DEM,其绝对高程精度是±16 m ．

2.3 矿区部分工作面及其水准数据

2018年10月至2019年1月,研究区开采主要集中在第三、四采区,根据6景Sentinel-1A SAR影像的时间跨度(2018年10月18日-2019年1月22日),本研究收集、整理和分析了三、四采区5个工作面的开采情况和232个水准点监测数据,用于比较、分析和验证全分辨率和多视法DInSAR矿区地面沉降监测的准确性和可靠性.图2显示了5个工作面和232个水准点的分布和范围,表2提供了5个工作面的基本开采信息.其中,各水准点的编号、倾向线按自南到北的方向顺序进行编号,走向观测线按工作面推进的方向顺序进行编号．

图2 三、四采区部分工作面及其水准点布设

表2 三、四采区5个工作面的基本信息

3 矿区沉降监测原理与数据处理流程

首先选取覆盖完整研究区形变前后的两幅SAR影像,分别作为主、辅影像,根据二者之间的成像几何关系进行复图像配准处理,进而计算同名点偏移量,使主、辅影像的像元对应同一地面单元[19-20].对配准后的影像数据作复共轭相乘,获取干涉条纹图.此时的干涉图中既包含目标地物的形变信息,又包括一系列干扰相位信息.为抑制各种噪声对干涉图质量的影响,对其进行滤波,继而引入外部DEM,模拟地形相位,与SAR干涉结果进行差分,消除地形相位,得到由地表形变引起的相位.但由于三角函数的周期性,此时的差分干涉相位缠绕在[-π,π]之间,因此需要进行相位解缠,得到代表真实地表形变的相位Ødefo,进而计算雷达沿视线方向的地表形变Δr：

(1)

式中：λ为波长；Δr为雷达视线向上的地表形变位移．

雷达视线向上的地表形变位移Δr：

(2)

雷达在垂直方向上的地表形变位移Δrvertical：

(3)

式中：θ为雷达视线的入射角．

技术具体流程如图3所示．

图3 数据处理流程图

4 结果及分析

将表1中5组干涉对的主、辅影像配准,利用外部DEM去除地形相位,再分别利用全分辨率和5∶1多视法进行干涉处理,得到差分干涉图.然后对差分干涉图做Goldstein滤波处理,得到滤波后如图4所示的差分干涉图．

图4 研究区全分辨率、多视处理滤波后的差分干涉图

分析图4可以看出：

全分辨率与多视法DInSAR所得到的差分干涉条纹的位置和分布非常一致.但从视觉效果上来说,全分辨率的差分干涉条纹在距离向(干涉图列方向上)被拉伸,可以更好地展示出地表形变的细节信息,从而也造成该方法得到的沉降中心的差分干涉条纹不如多视处理后的差分干涉条纹完整、连续．

图4d和4e的左侧均存在有规律的、类似系统误差的干涉条纹,这是由SAR影像自身造成的,因为2018年12月17日获取的SAR影像上相应位置就存在这种相同的条纹.全分辨率与多视法DInSAR均无法将其从差分干涉图上去除．

用最小费用流法对图4所示的差分干涉图进行相位解缠,由差分干涉相位值恢复为真实相位,然后根据公式(2)和(3),计算雷达视线向和垂直向上的地表形变量,最后对其地理编码,得到如图5所示的矿区垂直向地面沉降图.需要说明的是,为了便于对比和分析,图5中的沉降图都进行了距离向和方位向(列方向和行方向)上分辨率为15m×15m的重采样处理．

图5 全分辨率与多视法DInSAR监测的矿区地面沉降对比

结合表2,分析图5可以看出：

2018年10月18日至2018年11月11日(图5a),工作面4701的开采工作刚刚结束,全分辨率DInSAR在该工作面上监测到零星的-10～0 mm的沉降区域,多视法DInSAR 则监测到了较大的、较明显的-10～0 mm的沉降区域； 工作面3107的开采工作已经结束2个月,而4702的开采工作进行了3个月,全分辨率与多视法DInSAR均监测到这两个工作面上存在较大的、较明显的-10～0 mm的沉降区域； 工作面4707的开采工作进行了一个月,全分辨率DInSAR在该工作面上监测到零星的-10～0 mm的沉降区域,多视法DInSAR 则监测到了较大的、较明显的-10～0 mm的沉降区域； 工作面4705的开采工作进行了7个月,全分辨率与多视法DInSAR均监测到该工作面是主要的沉降区域,已形成沉降盆地,最大沉降的位置分别在a1和b1处,最大沉降值分别为-47 mm和-45 mm.但由于两幅SAR影像的空间基线(106 m)和时间基线(24 d)都较大,在该沉降盆地,全分辨率与多视法DInSAR均未得到连续的干涉条纹,监测的沉降盆地也比较凌乱、不连续.总之,全分辨率与多视法DInSAR监测的矿区地面沉降与矿区工作面的开采进度相符,但二者监测的沉降值和沉降范围并不相同．

2018年11月11日至2018年11月23日(图5b),全分辨率与多视法DInSAR均监测到矿区主要的沉降区在5个工作面上方,矿区的地面沉降仍在持续加重,二者监测到的主要沉降区的位置、沉降量级也比较一致,例如,二者监测的最大沉降位置a1和b1非常一致,而且都监测到该处存在比较明显的沉降盆地,但沉降量的大小略有不同； 全分辨率与多视法DInSAR都监测到了整个三、四采区都存在-20～0 mm的地面沉降,这有可能是由于DInSAR的两幅SAR影像数据质量的影响,也有可能是数据处理过程中有些误差难以去除造成的．

2018年11月23日至2019年1月22日(图5c-5e),整体而言,全分辨率与多视法DInSAR监测的矿区地面沉降的分布、形态和变化趋势等与矿区工作面的开采情况相符.但从细节上来看,二者监测到的最大沉降值、最大沉降的位置、范围的大小等并不相同.分析图5e发现,由于两幅SAR影像的时间基线(36 d)较大,两幅SAR影像的相干性较低,而且矿区沉降较为严重,全分辨率DInSAR所得到的差分干涉条纹的连续性较差,在工作面4702上方、4705右上方的沉降较为严重的沉降盆地,相位解缠和形变反演结果的连续性较差,沉降监测结果的精度变低.相比较而言,多视法DInSAR的形变监测结果较好,在4705右上方的沉降盆地得到了较为连续、精度较高的地面沉降监测结果.因此,全分辨率与多视法DInSAR监测到的最大沉降值出现的位置相差非常大．

为定量对比分析5个监测时间段内全分辨率与多视法DInSAR监测的矿区地面沉降情况,表3中对全分辨率和多视法DInSAR监测的各时段的最大沉降量、最大沉降位置进行了统计．

表3 全分辨率与多视法DInSAR监测的最大沉降值及其位置统计表

结合图5,分析表3可以看出：前4个时间段内全分辨率与多视法DInSAR监测的最大地面沉降位置基本一致,全分辨率DInSAR监测的地面沉降最大值均大于多视法DInSAR监测的地面沉降最大值,最大绝对差值为12 mm.在第5个时间段内,全分辨率DInSAR在b5处的沉降盆地没有得到连续的差分干涉相位和正确的相位解缠结果,从而导致了监测能力和精度降低,其监测到的最大沉降位置a5与b5并不在一个沉降盆地,二者相距甚远,DInSAR监测的a5处的地面沉降最大值小于多视法DInSAR监测的b5处的地面沉降最大值.总体上看全分辨率方法监测到的沉降最大值更有效．

5 精度验证与对比分析

为了验证全分辨率与多视法DInSAR监测的矿区地面沉降结果的精度与可靠性,将5个时间段(2018年10月18日至2019年1月22日)全分辨率与多视法DInSAR监测的矿区地面沉降累计结果与济阳矿区插值后同时间段经过精密水准测量所获取的地面沉降结果进行了对比分析.叠加水准点的沉降图如图6所示,沿观测线水准点的沉降量对比如图7所示,其均方根误差如表4所示．

图6 水准点与2018年10月18日至2019年1月22日第三、四采区矿区地面沉降

图7 沿观测线沉降结果对比图

表4 各观测线上全分辨率与多视法DInSAR监测的地面沉降的均方根误差

结合图6,分析图7和表4可以看出,全分辨率与多视法DInSAR监测的矿区三、四采区地面沉降与水准监测的地面沉降趋势在空间上基本一致,其精度也非常接近.图7a,7d中Q和L观测线上全分辨率DInSAR监测的沉降曲线波动要比多视法DInSAR监测的沉降曲线的波动明显,但均方根误差较小,图7b,7c中G和H观测线上两种DInSAR监测的沉降曲线比较平滑,全分辨率DInSAR比多视法DInSAR监测结果的均方根误差大．

6 讨论

全分辨率与多视法DInSAR监测的矿区地面沉降累计结果与济阳矿区插值后同时间段的精密水准测量所获取的地面沉降结果在沉降盆地中心位置有明显偏移,最大沉降量差异非常显著,这主要由InSAR相位干涉测量和侧视成像的固有局限性引起的．

InSAR相位干涉测量的固有局限性主要是指InSAR干涉测量监测形变的能力和精度很大程度上取决于波长、入射角、地面分辨率、重访周期等SAR成像固有参数.当同一像素单元内发生的形变量超出干涉波长的一半时,干涉测量结果无法反映地表真实形变[21].据此,Massonnet等[22]提出形变梯度的概念,并给出了差分InSAR监测到的最大形变梯度的理论公式,即雷达波长的一半与像素大小的比值.在实际应用中,由于地表形变信息极易被噪声掩盖,InSAR监测到的最大形变梯度真实值远小于理论值[23]．

InSAR侧视成像的固有局限性主要是指单轨道(升/降轨)DInSAR技术仅能监测得到监测区域一维地表形变信息,它仅仅代表了传感器成像方向上的真实空间位移,即雷达视线(Line of Sight,LOS)方向上的投影[24].由于现有在轨SAR 卫星的飞行轨道一般为近极地轨道,且多为右侧视成像,这种观测方式导致了DInSAR技术在垂直方向和水平方向(东西向和南北向)的形变监测敏感度存在差异,一维视线LOS向形变场并不能完全反映监测对象的真实形变状态,即存在视线向模糊问题[24].因此,直接由一维LOS向上的形变量转换成常用垂直向的形变量时,就会存在由视线向模糊引起的转换误差的问题,当研究区水平形变较大时,由单轨道(升/降轨)DInSAR技术监测的垂直向形变量值与实际地表垂直向的形变量值的差异也就增大,而且会出现形变中心偏移的问题[25]．

7 结论

本研究利用全分辨率和多视法DInSAR分别对山东济南某矿区5个Sentinel-1A SAR影像干涉对进行处理,结合矿区5个工作面开采信息和4条水准观测路线的水准监测数据进行精度验证和对比分析,得到如下主要结论：

2018年10月18日至2019年1月22日,全分辨率与多视法DInSAR均监测到矿区5个工作面上方的地面沉降逐渐加重,两种方法监测得到的矿区地面沉降的分布、形态和变化趋势比较一致,与矿区工作面的开采情况相符.但二者监测到的最大沉降值、最大沉降的位置、范围的大小等并不相同．

理论上来说,全分辨率DInSAR由于没有进行距离向上的多视滤波处理,其监测大梯度地面沉降的能力应该更强,但实际上,由于大梯度地面沉降区域的相干性较低,全分辨率DInSAR所得到的差分干涉条纹的连续性变差,从而使相位解缠和形变解算的精度受到较大的影响,监测得到的地面沉降最大值与多视法DInSAR监测的地面沉降最大值并没有较大差异．

在水准观测线上,全分辨率与多视法DInSAR监测得到的地面沉降分布曲线变化趋势非常一致,沉降量级相差较小,与水准监测的地面沉降分布曲线基本一致,二者监测的沉降盆地中心位置和最大沉降量也比较一致,但与水准监测的沉降盆地中心位置有明显偏移,最大沉降量也有明显差异.全分辨率DInSAR并未明显提高多视法DInSAR监测大梯度地面沉降的能力．