杨泽发 马泽林 乔思宇

(中南大学地球科学与信息物理学院,湖南 长沙 410012)

地下采煤容易导致地表产生移动变形,从而破坏建(构)筑物等基础设施。我国村庄下采煤较为常见,因此,科学开展开采损害鉴定是解决矿方和建筑物所有权人赔偿纠纷的重要途径。在开采损害鉴定中,客观真实地划定沉陷边界至关重要[1-2]。一般而言,待鉴定建(构)筑物坐落于沉陷边界以外,则可基本判定其损害不是由地下开采导致;反之,则需要根据开采沉陷量级及建(构)筑物结构和空间变化特征进一步研判其损坏情况,以便据此划分损坏等级。传统开采沉陷边界划定方法主要分为两大类[3]:① 根据地面实测沉降划定;② 在缺乏地面实测沉降的矿区,结合地下工作面开采进度与概率积分法[4]等理论模型预测地表沉降,并据此划定开采沉陷边界。然而,这两类传统方法均需依赖矿区地表实测资料(包含实测沉降或地下开采参数),对于缺乏实测资料的矿区则难以适用。

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)是一种微波遥感技术,具有大范围、高精度、高空间分辨率的地表形变监测能力。目前,该技术已成为矿区形变监测的重要手段[5-7]。倘若能基于InSAR形变观测值划定沉陷边界,则至少可带来两大优势:① 相较于传统矿山“线”状沉降观测值,InSAR空间分辨率更高,因而能更直观地划定沉陷边界;② InSAR技术可利用存档数据“回溯”矿区历史变形(最早可追溯到1991年),可为缺乏实测资料的矿区沉陷边界划定提供客观、精确的形变资料。

然而,受SAR斜视成像的限制,InSAR仅能监测地表真实变形矢量在卫星视线方向(Line-of-sight,LOS)的一维投影[8],而非真实沉降。此外,受观测噪声的影响,不同地区的InSAR观测误差存在差异。因此,如何基于InSAR观测值可靠地划定沉陷边界是InSAR支持下开采损害鉴定的关键。事实上,国内外学者在基于InSAR观测值的沉降盆地边界圈定方面已开展了部分探索。这些工作大都基于忽略东西和南北方向水平移动的InSAR沉降估计值,并以下沉1 cm为阈值划定沉陷边界[9-12]。这些方法用于大致圈定沉陷边界尚可,但若用于辅助开采损害鉴定则存在以下局限:① 忽略矿区东西和南北方向水平移动可能导致较大误差,降低InSAR沉降估计值的精度;② 在现有规程中,1 cm阈值作为沉降边界是以水准观测垂直沉降为前提,并根据不同等级水准监测精度综合确定。然而,与水准等具有严密检核条件的观测技术不同,InSAR形变监测精度受众多因素影响(例如时空失相干、大气延迟、噪声、处理参数等),导致其在不同地区精度可能存在显著差异。若直接将1 cm InSAR沉降等值线作为开采沉陷边界则不能充分考虑不同区域、不同InSAR数据监测精度差异,故而难以保证沉陷边界的可靠性。

为了克服上述局限,本研究提出了基于InSAR技术的开采沉陷影响边界划定方法。首先,仅忽略InSAR视线向形变不敏感的南北向移动分量,并通过融合升降轨InSAR观测值获取矿区历史沉降,提高InSAR沉降监测精度;然后,以非形变区监测值为样本自适应评估InSAR沉降估计值的整体精度,并以置信区间为准则划定矿区沉陷边界。该方法有助于提升InSAR矿区沉陷边界划定精度,推动InSAR技术更科学地指导开采损害鉴定。

1 研究方法

本研究提出的基于InSAR技术的开采沉陷影响边界划定方法主要包括3个步骤:首先,利用时序In-SAR技术获取矿区地表视线向一维形变;然后,联合升降轨InSAR视线向形变观测值估计地表开采沉降值;最后,通过非形变区InSAR沉降估计值评估其观测精度,并据此自适应划定沉陷边界等值线。

1.1 SBAS-InSAR技术基本原理

时序InSAR算法较多[13],本研究以矿区较为常用的小基线集InSAR(Small Baseline Subset InSAR,SBAS-InSAR)[14]为例,分析其基本原理。首先假定N+1景同一轨道SAR影像覆盖待监测区域,其观测时间按照先后顺序可表示为T0,T1,…,TN;然后根据观测时间和空间基线居中的原则从SAR影像中选取一景影像作为超级主影像,并将剩余影像与其配准;最后设置合理的时空基线阈值,根据配准的SAR影像组成小基线集InSAR干涉对(其数量为M)。

在对构建的M对InSAR干涉对都进行差分干涉处理后,即可获得M个LOS向观测值。根据小基线集时空基线网络,构建地表LOS向形变速率与观测值的观测方程组:

式中,矩阵A为由不同InSAR干涉对时间差组建的设计矩阵,由根据构建的短基线InSAR对和SAR成像时间组成;V为估计的相邻SAR成像期间的形变速率V=[V1,V2,…,VN]T;ΔdLOS为M个InSAR干涉对获取的视线向观测值ΔdLOS= [ΔdLOS1,ΔdLOS2,…,ΔdLOSM];ε为观测误差。根据最小二乘或奇异值分解方法求解形变速率V,在此基础上,根据形变速率与时间的关系便可估计SAR成像时刻的时序形变:

1.2 联合升降轨InSAR的矿区开采沉陷监测

InSAR一维视线向形变是地表真实三维形变分量的投影[15],即:

式中,θ为SAR传感器入射角;α为飞行方位角;dv、de和dn分别为垂直、东西和南北向形变。

在现有的InSAR沉陷边界圈定方法中,直接忽略东西向和南北向水平移动估计沉降[7,16],简称NOV(Neglecting Horizontal Movements)方法,即:

考虑到煤矿区水平移动量级较大[17],NOV方法估计的沉降值精度不高。在目前SAR卫星近南北向飞行轨道下,南北向水平移动分量对InSAR一维视线向形变贡献较小,因此,本研究仅忽略南北向水平移动分量[18-19],而非传统NOV方法中同时忽略南北和东西方向,则垂直沉降可联合升降轨InSAR视线向形变观测值进行估计,即:

式中,θasc和θdes分别为升降轨SAR入射角;αasc和αdes分别为升降轨SAR飞行方位角;dasc和ddes分别为升降轨InSAR视线向形变观测值。

相较于NOV方法,本研究采用的垂直沉降估计方法仅忽略了InSAR视线向形变观测值不敏感的南北向水平移动,因此,理论上精度优于NOV方法。

1.3 基于非形变区统计的沉陷边界划定

在获得地表沉降盆地后,阈值的选取在很大程度上决定了沉陷边界划定的可靠性。传统方法大都基于水准观测值,并采用地表下沉1 cm为阈值。然而,InSAR形变观测值受大气延迟、解缠噪声等多误差源影响,不同地区InSAR观测值的精度存在差异。若使用统一的阈值划定沉陷边界,则无法保证结果精度,但现有研究极少讨论InSAR沉陷边界划定阈值。

本研究假定不受地下开采影响的区域不存在移动变形,因此,在这些区域中,升降轨融合解算的沉降观测值理论上是由于InSAR观测误差导致的。基于该假设,沉降观测值的整体标准差δ(dv)可近似表示为[20]

式中,dv(j)、为非形变区(可通过目视大致圈定,其像素个数为N)第j个升降轨InSAR沉降解算值和垂直向形变均值。

根据随机误差理论,一倍标准差(Standard Deviation,STD)表示的置信水平仅为67%左右,若直接使用一倍STD值划定沉陷边界可能出现边界过大的现象。为此,本研究将InSAR沉陷边界划定阈值dv_boundary设定为

式中,β为置信区间对应的STD倍数,例如在90%、95%和99%的置信区间下,β可分别取1.65、1.96和2.58。

2 试验验证

2.1 模拟试验

2.1.1 升降轨InSAR形变观测值模拟

本研究首先通过模拟试验验证方法的可行性和可靠性。具体地,首先利用概率积分法(我国最常用的开采沉陷预计模型之一)模拟地下工作面开采导致的地表三维形变;然后,将模拟的三维形变分别投影至升轨(θasc=39°,αasc=350°)和降轨(θdes=34°,αdes=190°)视线方向。为了使模拟试验更符合实际情况,本研究在投影的升降轨视线向形变中加入了大气紊流延迟和噪声导致的观测误差,以此模拟含误差的升降轨InSAR视线向形变观测值,结果如图1(a)和图1(b)所示。

2.1.2 沉陷边界InSAR划定

为了能够划定沉陷边界,本研究首先根据模拟的升降轨InSAR视线向形变观测值,并采用式(5)估计了地下开采导致的垂直沉降值;之后,以非形变区In-SAR沉降估计值为样本,采用式(6)计算了沉降估计值的整体STD值为0.45 cm;最后,分别令β为1、1.65、1.96和2.58,以此获得置信区间分别为67%、90%、95%和99%的沉陷边界,其结果如图1(c)所示。

图1 采煤沉陷区地表形变等值线Fig.1 Surface deformation contour in coal mining subsidence area

由图1(c)可知:受观测误差影响,采用1倍STD值划定的沉陷边界呈现出较大波动,且比参考边界范围大很多。采用1.65倍和1.96倍STD值作为沉陷边界阈值时虽然比1倍STD值划定的边界更接近于参考边界,但部分区域(例如沉降盆地的西南处)与参考边界的差异仍然较大。相较于其他3个阈值,根据2.58倍STD值划定的沉陷边界与参考边界最为接近。

2.2 真实试验

2.2.1 研究区概况及InSAR数据处理

本研究选取河北省唐山市某煤矿(图2中五角星符号标志处)作为真实案例再一次验证提升方法的可行性和可靠性。该矿区位于唐山市郊区,地表主要为村庄和农田。为了尽可能地降低农田区植被季节性变化对监测结果的影响,收集了研究区2017年11月—2018年5月(主要跨越冬季和春季)16景升轨和16景降轨Sentinel-1 SAR影像数据集,升降轨Sentinel-1影像的覆盖范围如图2所示,两者获取时间仅相差1 d。

图2 研究区位置及SAR数据覆盖情况Fig.2 Location of study area and SAR data coverage

采用SBAS-InSAR技术分别处理收集的升轨和降轨Sentinel-1影像集,从而获得研究区升降轨In-SAR视线向形变观测值。具体地,首先设定时空基线阈值为30 d和200 m,并根据升降轨SAR影像时空基线分布选取小基线InSAR干涉对。之后,采用距离向和方位向5∶1多视处理削弱噪声,并采用差分InSAR干涉流程处理获得小基线集解缠相位图。最后,采用时空域滤波分别削弱大气延迟影响,并采用奇异值分解方法获得研究区升降轨InSAR视线向形变观测值。

2017年11月—2018年5月研究区升降轨视线向累计形变分别如图3(a)和图3(b)所示。分析可知:2017年11月—2018年5月研究区地表出现了两个明显的沉降盆地,其中,区域1的沉降盆地能够被较为清晰地监测到,区域2的沉降盆地中心部分受地表水体影响(导致InSAR完全失相干)则未能被监测出。结合工作面采掘进度情况可知:研究区共有两个工作面开采,区域1的工作面于2017年6月开始采动,直至2019年1月结束,而区域2的工作面在2017年6月前已开采结束。据此可知,区域1的沉降盆地由地下工作面开采导致,区域2则可能由于残余沉降或邻近开采引起的上覆岩层“活化”所致。

2.2.2 沉陷边界InSAR划定

为了便于划定沉陷边界,首先融合升降轨Sentinel-1 InSAR监测的视线向形变观测值,在忽略南北向水平移动分量的基础上解算开采沉陷盆地(图3(c));然后以非形变区沉降值为样本,统计了非沉降区估计值的均值和标准差(STD)分别为0.1 cm 和0.9 cm;之后,以1、1.65、1.96和2.58倍的STD值为阈值划定了置信区间分别为67%、90%、95%和99%的沉陷边界(图3(c))。由于缺乏该时间段内监测区域的实测沉降值,因此难以科学地评估不同阈值划定的沉降边界的可靠性。鉴于此,本研究以唐山市其他地质采矿条件类似矿区的实测边界角(即56°)为依据,结合研究区域工作面参数(平均倾角为10°,平均采深为785 m)划定了参考沉陷边界。需要指出的是,受煤层倾角影响,边界角在工作面走向、倾向上山和倾向下山方向可能存在差异,但由于研究区域的工作面为近水平煤层,其在走向、倾向上山和倾向下山方向的边界角差异通常较小。因此,本研究忽略了该差异,直接将3个方向的边界角全部设定为56°。

图3 2017年11月—2018年5月研究区形变特征Fig.3 Deformation characteristics of the study area from November 2017 to May 2018

由图3(c)进一步分析可知:1倍STD值作为等值线阈值划定的沉陷边界远大于参考边界。通过不断增大等值线阈值,划定的沉陷边界逐渐接近于参考边界。当等值线阈值为2.58倍的STD值(即99%置信区间)时,划定的沉陷边界在移动盆地西部与参考边界较为接近。然而,需要说明的是,该图中东南处划定的沉陷边界与参考边界存在较大差异的原因在于该处为老采空区(图3(a)区域2),受区域1处邻近工作面开采的影响导致区域2处地表再次下沉,造成该处地表沉降较大,故而使得该处划定的沉陷边界大于参考边界。

3 讨 论

3.1 等值线阈值对于InSAR沉陷边界划定的影响

现有研究大都采用下沉1 cm InSAR沉降等值线作为划定沉陷边界的阈值,并未考虑不同地区InSAR沉降观测值误差的差异。本研究中,模拟试验部分1.96倍STD值(图1(c))和真实数据部分的1倍STD值(图3(c))所对应的阈值约1 cm。从图中可以看出,采用下沉1 cm等值线阈值划定的沉陷边界均与参考边界存在显著差异,且该差异的空间分布与InSAR解算的沉降误差有关。即采用恒定的沉降等值线阈值划定沉陷边界的可靠性难以保证。

理论上,受地表覆盖变化差异、InSAR数据监测能力差异以及数据处理方法可靠性差异等众多因素影响,不同地区InSAR监测的沉降值精度也不同。因此,难以采用统一的InSAR沉降等值线阈值划定沉陷边界。通过模拟和真实案例分析可知,根据In-SAR观测数据自身包含的近似精度信息(本研究选用的是统计标准差),通过设定合适的置信区间能够自适应地判断InSAR沉降等值线阈值,并据此划定沉陷边界。本研究模拟试验和真实案例试验中,选用99%的置信区间(即2.58倍STD)划定的沉陷边界与参考边界较为接近。

3.2 忽略InSAR水平移动贡献对沉陷边界划定的影响

如前所述,忽略InSAR视线向形变观测值中东西向和南北向二维水平移动贡献,从而利用InSAR视线向形变观测值估计开采沉陷值是目前较为常用的方法(即NOV方法,式(4))。因此,现有研究在划定沉陷边界时大都基于该方法估计开采沉降值。事实上,该方法估计的沉降值可能含有较大误差,影响了据此划定的沉陷边界精度,从而“误判”或“错判”了沉陷边界。

为了分析该问题,本研究首先以2.1.1节模拟的升降轨InSAR视线向形变为基础,使用NOV方法分别估计了升轨和降轨InSAR垂直沉降。然后,选取非形变区沉降值为样本,分别计算NOV方法获取的升降轨InSAR垂直沉降的STD值。最后,采用99%的置信区间划定了沉陷边界,如图4(a)所示。从图中可以看出,在升轨观测条件下,基于NOV沉降值划定的沉陷边界在盆地东南处出现了明显的向内收缩现象。在降轨观测条件下,以NOV沉降值为基础划定的沉陷边界则在盆地西北处出现了明显的向内收缩现象。这就意味着,若采用NOV方法划定的沉陷边界辅助开采损害鉴定,则可能“漏判”沉降盆地东南(升轨)或西北(降轨)区域的建(构)筑物损害程度。相较于NOV方法,基于升降轨InSAR观测值融合的垂直沉降估计方法划定的沉陷边界基本与参考边界吻合,未出现显著的系统性偏移。类似的现象在真实数据试验中也存在,如图4(b)所示。由该图可知:以降轨NOV沉降估计值为基础采用99%置信区间为指导划定的沉陷边界在西北处存在显著的收缩现象;然而,受老采空区上覆岩层“活化”的影响,导致以升轨NOV沉降估计值为基础采用99%置信区间为指导划定的沉陷边界在东南处未出现显著收缩现象。采用本研究方法以99%置信区间为指导划定的沉陷边界(融合升降轨InSAR观测值估计的沉降值为基础,按照2.58倍STD值划定的沉陷边界,即图4(b)中升降轨边界 )与参考边界在西部区域吻合较好,而同样受老采空区上覆岩层“活化”的影响导致其在东南处划定的沉陷边界较大于参考边界。

图4 以99%置信区间划定的沉陷边界Fig.4 Subsidence boundary defined with 99% confidence interval

4 结 论

本研究提出了一种基于InSAR技术的开采沉陷影响边界划定方法,并通过模拟数据和真实案例试验验证了该方法的可行性和可靠性。所得结论如下:

(1)以升降轨InSAR观测值融合解算的矿区沉降为基础,通过合理设置等值线阈值能够实现不依赖地面实测资料的沉陷边界划定。

(2)现有研究中采用下沉1 cm的恒定等值线阈值未考虑不同地区、不同InSAR数据监测的垂直沉降精度差异,因而难以保障沉陷边界的可靠性。

(3)本研究以非形变区沉降值为统计样本评估InSAR沉降估计值的整体STD值,并选用2.58倍STD值(即99%置信区间)为阈值划定的沉陷边界与参考边界吻合度较高。

(4)以NOV方法估计的沉降值为基础划定的沉陷边界会出现向内收缩现象,可能会“误导”开采损害鉴定。