基于D-InSAR 技术的煤矿区地面沉降监测研究*

宋继德，邸志众

(1．陕西省地质矿产勘查开发局 第二综合物探大队，陕西 西安 710016;2．北京中勘迈普科技有限公司，北京 100086)

0 引言

InSAR(SAR Interferometry 或Interferometric SAR)，即合成孔径雷达干涉测量，是指利用SAR 数据中的相位信息进行干涉测量处理，结合雷达参数和卫星位置信息反演地表三维及其变化信息的遥感手段［1，2］。D-InSAR 即差分合成孔径雷达干涉测量，是指对干涉相位进行差分处理提取变化信息的干涉测量手段，是InSAR 技术的延伸。InSAR 技术已被广泛用于大区域地形测图和地表信息提取。常见的InSAR 方式有顺轨(Along track)，交轨(Across track)和重复轨(Repeat pass)。由于地表形变监测是测量一定时间间隔内的地表形变量，因而重复轨是目前应用的主要方式，具体是指利用同一地区不同时刻获取的SAR 数据进行差分干涉处理以提取地表形变信息。由于InSAR具有全天候、全天时、覆盖面广和高精度获取地表信息的能力，因而在地球科学的诸多领域如地震、火山、矿区塌陷、地面沉降、滑坡、冰川等灾害中逐渐得到了广泛的应用［1，5］。

本文以山东省重点地区矿山开采的遥感动态监测为研究内容，展开济宁、兖州地区地面沉降InSAR 监测试验研究，围绕济宁和兖州地区煤矿开采、地下水开采等多种因素作用下地面沉降(陷)的分布与发展状况，介绍以D-InSAR 技术为主要手段的地面沉降调查与监测工作的研究情况与成果。

1 差分合成孔径雷达干涉测量原理与数据处理方法

1．1 基本原理

重复轨干涉测量的几何原理［2］，如图1 所示。图中O1、O2分别为两次获取SAR 信号的天线位置，B 为两天线间空间距离，称为空间基线。

图1 重复轨D-InSAR 几何示意图Fig．1 Geometric schematic diagram of repeat pass D-InSAR

设基线与水平方向的夹角为α，H 为平台高度，地面点P到两天线间的距离分别为R1和R2，θ 为第一副天线的侧视角，P 点高程为h。将基线沿雷达视线向分解，垂直于视线向分量为B⊥，称为垂直基线。雷达对地物观测时，接收到地物的后向散射回波信号可以表示为:

两次观测所获取的SAR 影像分别称为主影像(Master)和辅影像(Slave)，接收到P 点的雷达信号分别为:

雷达后向散射回波信号中的相位信息包含两部分信息，一是与往返路径相关的相位，另一部分是目标后向散射相位，可以表示为:

两次获取的相位分别为:

式中:ΔR 为路程差，与雷达与目标的距离相关。若两次雷达获取的相位中的散射相位保持稳定，即:

则经过共轭相乘处理后得到的干涉相位为:

干涉相位φint的贡献值(Contributions)由以下5 个因素构成:平地相位δφflat;地形相位δφtopo;两次观测期间观测目标沿雷达视线向(LOS)移动引起的相位变化δφmov;电磁波传输过程中大气波动引起的相位延迟δφatm;其他噪声相位δφnoise。则干涉相位可以表示为:

1．2 数据处理方法

根据式(8)，要得到反映地形信息的干涉相位δφtopo，应当首先从干涉相位和中消除平地相位的影响，其次从中减去其他各项得到地形相位，再进行相位解缠和相位到高程的转换，将生成的高程信息转换到特定的坐标系统下，生成DEM 数据［1，6］。由于地形形变相位、热噪声相位、大气波动相位以及轨道偏差引起的相位波动可以通过合理的处理方法来消除或者降低。对于ERS-1/2，其特殊的“串行模式”(Tandem Mode)SAR 数据是生成DEM 的最佳选择［6］。因而，利用InSAR 获取DEM 数据的处理流程，见图2。

图2 InSAR 数据处理流程Fig．2 Data processing flow of InSAR

1．3 D-InSAR 地表形变监测技术及应用

随着D-InSAR 技术研究的不断深入，利用其进行地表形变信息获取的数据处理方法也在不断进步。针对不同形变体的变形特征，产生了不同的监测方法。按照监测目标空间分布的不同和数据处理方式的差别，可将InSAR 信息提取方式分为两类:基于点目标的D-InSAR 方法和基于连续目标的D-InSAR 监测方法。

基于单个点目标的方法以角反射器方法为代表，主要针对的是形变区域无相干目标或存在极少数相干目标的情况，监测区域较小，适合于局部变形监测，如滑坡体、危岩活动监测等［1，7，8］。

连续目标方法指的是在空间大面积范围内，长时间条件下依然存在较多相干目标的情况［1］。常规差分干涉测量方法难以解决多时相动态监测问题。缓慢连续变形体的移动监测依赖于长时间连续监测，对于重复轨监测，长时间周期会引起失相干，这种情况下，对于大区域连续形变体的移动监测则比较困难。由于监测对象多是与人们生产和生活紧密相关的城市或城镇，这些区域在长时间间隔下依然保持了很好的相干性，因而，可以通过多时相分析，获取这些区域的地表形变信息。由此，则产生了累积纹图法(Stacking Interferograms)方法和永久性散射体(Permanent Scatterers)技术。这两种方法通过长时间序列分析，弥补了常规重复轨干涉测量的不足，提高了时间采样率，能更准确的获取形变随时间变化而产生的形变。可以将形变特征分为线性和非线性分别进行处理［3，4，8］。但由于永久性散射体技术对SAR 数据数目要求较高，一般情况下难以完成大面积形变。

2 研究区域地表形变InSAR 监测的技术方法

山东省兖州地区为典型的煤矿开采区，开采沉陷诱发的矿区环境破坏不仅出现在井田开采上覆地区，同时也影响到了城市和城镇的地面稳定性，因而矿山沉陷InSAR 监测的首要任务是明确2004 ～2010 年开采沉陷的变化特征，进而选择合理的数据处理方法和相关技术参数。

针对缓慢变形和快速沉降两种不同的形变类型，采用常规D-InSAR 技术与永久散射体干涉测量相结合的工作方法、利用永久性散射体干涉测量(PSI)技术相较其他技术，在剔除地形误差和大气波动对相位信号的影响方面的优势，且目前可获得的SAR 数据在精度和分辨率方面也足以满足PSI 对数据的要求。

2．1 永久散射体干涉测量(PSInSAR)技术

永久散射体干涉测量(Permanent Scatterers Interferometry SAR，PSInSAR)的核心思想是对永久散射体干涉相位进行时间序列分析，根据各相位分量的时空特征，估算大气波动，DEM 误差以及噪声等，将其从差分干涉相位中逐个分离，最终获取每个PS 的线性和非线形形变速率、大气延迟量(Atmosphere Phase Screen)以及DEM 误差。经PS 方法处理，获取的年度形变率的精度可以达到毫米级。其基本思想是基于大量的SAR 数据(20-30 景甚至更多)，从中筛选出具有稳定散射特性的相干点目标，构成离散点观测网络(较之常规的变形监测网密度更高)，通过分析PS 点目标相位变化获取地表形变状况。由于将永久散射体作为观测对象，降低了空间基线对相干性的影响，即使在临界基线的条件下，仍然可以通过分析PS 差分干涉相位的变化反演形变信息、提取地表形变速率和累积形变量［8］(见图3)。

PS 技术一般采用线性形变模型提取点目标对应的形变量，非线性模拟的处理过程复杂，而且非常耗时，限制了其用于大面积的形变测量，但随着处理技术的进步，处理时间将逐步缩小，处理的范围也可再进一步扩大。

PSInSAR 对数据量要求较高，只有SAR 影像个数达到一定的程度才能筛选出在整个时间跨度内具有稳定信号的PS 点，对于相干目标较多的地区，如城区等，由于地物在长时间间隔下保持了较高的相干性，数据量较少的情况下也可进行处理。一般情况下，至少必须满足5 个/km2才能完成相位解缠。

图3 PSInSAR 数据处理的基本方法Fig．3 The basic method of data processing of PSInSAR

PS 方法在处理过程中通过分析每个点目标对应的相位，将其中包含的因地表形变，高程以及大气波动引起的相位分量加以区分，最终得到各个点对应的形变相位。通过分析PS 数据集，估算每个干涉纹图中由大气波动引起的相位分量，此外可以确定出每一点精确的高程信息。完成上述过程后，在干涉处理时将二者对相位的贡献量加以去除。经PS 方法处理，获取的年度形变率的精度可以达到±0．1 mm/a。

2．2 可获取的主要雷达数据及参数分析

当前用于InSAR 技术研究和应用较多的是ERS -1/2、ENVISAT、RADARSAT -1/2 以及ALOS PALSAR 获取的SAR 数据，其空间分辨多为20 ～30 m，可归为中等分辨率雷达数据。与之相比较而言，近几年发射的雷达卫星以高分辨率为主，包括Cosmo- Skymed 及TerraSAR - X，其最高分辨率可以达到1 ～3 m。表1 所列为当前几种雷达传感器所提供的SAR 数据。

表1 SAR 数据及相关参数Tab．1 SAR data and relevant parameters

从雷达波段而言，目前主要的中等分辨率数据为C 波段数据(ENVISAT)。由于这些卫星具有较高的定位精度(可以达到5 cm)，空间基线解算结果可靠，数据相干性较好，因此，在多数地表形变监测中均可达到较好的监测结果。高分辨率雷达卫星TerraSAR-X 重复周期为11 d，采用X 波段，对地表形变等变化信息更为敏感，可以应用于形变速率极快的监测，矿山开采沉陷等非常适用。Cosmo-Skymed 与TerraSAR -X 相似，也属于高分辨率雷达数据，最高分辨率达1 m，扫描模式获取数据为3 m地面分辨率。在空间基线不大的条件下，L 波段的SAR 数据能够在很长的时间内保持较高的相干性，能满足变化幅度较大的形变监测。特别对于快速的地表移动监测，应用L 波段数据也较为合适。同时，由于这些数据分辨率的提高，空间采样单元的个数将会增大，能更为准确的描述地表形状。因此，在很大程度上满足了从短时间，变化速率高以及长时间，变化幅度大的形变监测。

2．3 工作方法与技术

针对兖州地区地表形变调查和监测的工作重点，主要目标为查明工作区地面沉降(陷)的历史状况和现状，调查兖州、济宁和邹城下属的多个矿区开采沉陷的基本特征和分布范围，开展兖州市地面沉降InSAR 监测。因而，采用常规D-InSAR 技术与永久散射体干涉测量相结合的工作方法进行。

2．3．1 工作区SAR 数据获取

采用的雷达数据为欧空局自2004 年11 月至2010 年10 月间接收的ENVISAT SAR 数据，其有效工作波段是C 波段，数据获取周期为35 d，空间分辨率为20 m，是目前常用且应用效果显著的SAR 数据。获取ASAR 数据共计24 景，如表2 所示。覆盖范围，见图4。

表2 ENVISAT 卫星数据列表Tab．2 List of ENVISAT satellite data

图4 ENVISAT 卫星数据覆盖范围Fig．4 Coverage of ENVISAT satellite data

区域性地面沉降监测实际使用的数据从2004 年12 月到2010 年10 月计20 景，根据区域地面沉降长时序监测需求形成干涉像对。

2．3．2 数据处理和分析

数据处理和分析系统以专业InSAR 数据处理软件，如GAMMA、Doris 和SARSCAPE 为主，选择利用ArcGIS 平台为数据后处理系统，针对工作区环境特点及矿区开采沉陷关键In-SAR 处理分析技术，分别利用常规D-InSAR 技术和永久散射体InSAR 技术分别针对快速沉陷和缓慢沉降开展监测。InSAR 数据处理过程主要包括最优干涉数据选择、干涉像对的基线估算、平地效应去除、相干性分析、大气效应分析、数据配准、离散点相位解缠，时间序列分析方法等。

2．3．3 地表形变信息提取

采用永久散射体InSAR 时间序列分析技术提取地表形变信息，并结合若干实地观测结果获取不同时间间隔的地面沉降速率图。针对矿山开采沉陷的发生特征，利用常规D-InSAR 技术开展数据处理，城市地表缓慢变化和活动断裂移动监测利用PSInSAR 技术开展数据处理与分析。

2．3．4 监测结果后处理

在多期地面沉降信息分析与综合的基础上编制工作区地面沉降(陷)速率图和累积沉降量图。

3 区域地面沉降与开采沉陷InSAR 监测与分析

3．1 煤矿区开采沉陷历史状况InSAR 动态监测

考虑到矿山塌陷的大变形过程对相干性的影响，以及季节性植被覆盖变化，选择冬季干涉像对，且间隔时间较短的数据进行D-InSAR 数据处理，共获取8 个能清晰反应矿山沉陷的干涉纹图，用于反映工作区矿区开采沉陷2011 年之前历史状况的动态过程。

全区总的开采沉陷区主要分布于兖州、曲阜和邹城沿市泗河沿岸以及济宁市区周边。这些地区不同时段内差分干涉条纹密集，表明了各个煤矿在监测时段内开采沉陷的发生位置和范围。在2004 年11 月到2005 年2 月3 个月的时间内，共有4 个差分干涉纹图(图5a、b、c、d)显示了该时间段的塌陷变化过程，各个采区的干涉条纹密度随着时间间隔的增加而整体上增大。由041121 -050306 差分干涉纹图(图5b)可知兖州市城区东部与曲阜交界(序号2)的开采沉陷漏斗干涉条纹达5 个色周，沉陷量达到15 cm，表明在105 d内由于煤矿开采引起了最大15 cm的沉降，沉降呈现出盆地形状的分布特征，与传统开采沉陷的沉陷盆地分布相同。位于济宁南张镇西北(序号1)的沉陷漏斗在041121 -050130(图5a)的干涉纹图超过2 个色周，沉陷量超过6 cm。050103 和050306 两景SAR 数据相干性好，差分干涉纹图(图5d)各沉陷中心清晰，可以明显看到该沉陷漏斗中心向北偏移。

2005 年12 月到2006 年2 月有1 个相干性好的干涉纹图(图5e)，该干涉纹图表现了该时间段内各沉陷中心的沉陷状况，曲阜与兖州在泗河交界处的沉陷漏斗干涉条纹有5 -6 个色周，中心沉陷量达15 ～18 cm;济宁南张镇-刘堤头间(序号3)的沉陷漏斗条纹达4 个色周，中心沉陷量超过12 cm。2009 年12 月到2010 年10 月相干性良好的3 个干涉纹图(图5f、g、h)清晰地表现了该时间段内的开采沉陷状况。由091220 -100124(图5f)和100124 -100228(图5g)干涉纹图可知，位于曲阜与兖州在泗河交界处张家村煤矿一带(序号4)的沉陷范围在35 d的时间间隔内达到3 个色周的沉降变化，表明在091220 -100228共70 d的时间内沉陷量达到近18 cm，属于快速沉降阶段。济宁主城区西北前王村南沉陷中心(序号5)在091220 -100124 和100124 -100228 干涉纹图中分别达到了6 个和3 个以上色周，共计约10 个条纹，表明最大沉陷量超过30 cm。图5h 为直接利用091220 和100228 两景数据干涉后的条纹图，反映了105 d内整个沉降区的变化状况，是图5f 和图5g 反映的沉降量的总和。

总体上，干涉图的空间位置和分布形态一致，但局部干涉图的条纹清晰程度降低。这种结果产生的原因是由于时间间隔进一步增大，相干性下降引起的。同时由于沉陷量增大，局部出现非相干的移动变化，这种快速变化限制了利用InSAR 干涉图之间确定沉降量的精准程度。同时，位于楼庄和梁家营一带出现了局部回弹，干涉条纹与其它沉降引起的形态相反，变化量为1个色周，影响范围约2 ～3 km。

3．2 全区2004 ～2010 年区域性地面沉降InSAR 连续监测

利用长时间序列进行PSInSAR 处理，提取相干目标的沉降速率。研究区内的济宁、兖州等城市，地面沉降主要分为缓慢沉降以及煤矿开采沉陷。前者为缓慢微小变化，主要发生于济宁市城区以及城区周边，而后者则多发生于煤矿区，与采矿条件、煤层厚度与赋存条件，开采方式等因素相关。因而，在煤矿区开采沉陷动态监测的基础上，实现全区地面沉降的连续监测。

图5 工作区煤矿开采沉陷InSAR 动态监测结果图Fig．5 Dynamic monitoring results for applying InSAR in subsidence area caused by coal mining

利用2004 年12 月～2006 年02 月和2008 年9 月～2010 年10 月两个时间段的ENVISAT ASAR 数据进行PSInSAR 时序分析，得到两个时间段的平均沉降速率，如图6、图7 所示。监测结果表明:工作区总的沉降范围与矿区分布一致，分别位于济宁市任城区北部、兖州市东南部、曲阜市西部以及邹城市西北部;工作区内有两大明显的沉降区域，兖州市东南、曲阜市西部以及邹城市西北部沉降连接成片，沉降速率大于20 mm/a的区域超过180 km2，另一个明显的沉降区域为济宁市任城区北部。

2008 年09 月～2010 年10 月InSAR 监测地面沉降结果(见图7)表明两大沉降区沉降范围有所扩大，兖州—曲阜—邹城沉降区内沉降中心沉降速率有所减缓，但西部整体沉降加剧，东部有所减缓;济宁任城区沉降区与汶上县、嘉祥县交界处沉降明显减缓，主城区西北郊南张镇一带有所加剧，北部岱庄—栗乡一带有所减缓，但任城区东部和南部出现明显沉降中心。

3．3 全区累积地面沉降量分析

通过2004 年12 月至2010 年10 月共24 景ASAR 数据的PSInSAR 时序分析处理，得到6 年间工作区累积地面沉降量，如图8 所示。由InSAR 监测累积地面沉降量所示可知，在近6 年的时间内，全区最大累积沉降量超过250 mm，两个沉降区分别为济宁市北部煤矿区和曲阜、兖州交界沿泗河一线。

图6 2004．12 ～2006．02 InSAR 监测地面沉降速率图Fig．6 Land subsidence rate by applying InSAR from December 2004 to February 2006

济宁、兖州辖区约4 060 km2的面积，去除大面积水域共约3 922 km2，通过对累积地面沉降量空间插值分析统计可得各沉降范围影响面积，如表3 所示。

图7 2008．09 ～2010．10 InSAR 监测地面沉降速率图Fig．7 Land subsidence rate by applying InSAR from September 2008 to October 2010

图8 2004．12 ～2010．10 InSAR 监测累积地面沉降量Fig．8 The cumulative land subsidence amount by applying In-SAR from December 2004 to October 2010

表3 累积地面沉降影响范围统计表Tab．3 The statistical table of influence range about cumulative land subsidence

对济宁、兖州各沉降范围影响面积进行单独统计，各地区主要沉降主要分布，如表4、表5 所示。

表4 济宁市累积地面沉降影响范围统计Tab．4 The statistical table of influence range about cumulative land subsidence in Jining city

表5 兖州市累积地面沉降影响范围统计Tab．5 The statistical table of influence range about cumulative land subsidence in Yanzhou city

3．4 矿区范围与沉降范围比较

通过对累积地面沉降量空间插值分析得到该区域累积沉降等值线图，将济宁和兖州地区煤矿区矿权范围数据与实际开采沉陷范围进行叠加，得到如图9 所示的累积沉降范围与开采范围。图中左下角深灰色范围为主要煤矿区矿权划定范围，比较发现，两个主要沉降区分别对应着兖州和济宁地区两个主要矿区范围，其它无矿权范围地区均为较小沉降区。对比沉降量的大小可知，大于60 mm沉降量的地区为煤矿开采区、整个沉降量为60 mm等值线所包围的区域均为煤矿矿权覆盖范围;东山王楼煤矿沉降突出，中心区大于80 mm，影响范围较广。

图9 工作区累积地面沉降范围与矿区开采范围Fig．9 The mining range and the cumulative land subsidence range in study area

4 结论

本次研究工作针对兖州和邹城地区区域性地面沉降和煤矿开采区快速沉陷两种类型地表形变，利用2004 ～2010 年多期ENVISAT ASAR 雷达数据，应用PSInSAR 处理技术与常规D-In-SAR 差分干涉处理方法，获取了全区2 个不同时间段的区域性地面沉降速率和6 年时间内全区累积沉降量结果，并得到了多个时段矿区动态沉陷监测成果。监测结果查明了济宁、兖州地区地面沉降分布状况，统计得到了全区累积沉降量和影响地区范围，全面掌握了该区域目前地面沉降的总体分布态势和局部沉降特征。区域性地面沉降结果表明，济宁、兖州的地面沉降主要与矿山开采相关，城市活动等其他因素带来的地面沉降甚微;工作区内地面沉降显著区有两个，一个为兖州市东南、曲阜市西部以及邹城市西北部沉降区，另一个明显的沉降区域为济宁市任城区北部。比较两期监测成果，虽然最大沉降速率有所减小，但沉降严重区域(年沉降速率＞20 mm/a)有扩大的趋势。

煤矿开采沉陷是上述区域沉降的主要因素，由于煤矿开采和开采过程的输排水等，改变了工作区局部地质环境，出现地面破坏。这种破坏与大规模开采地下水引发的区域性大范围沉降漏斗不同，从空间分布上沉降中心比较集中，沉降速率比较大，局部甚至出现地裂缝和塌陷。因而，从监测手段上InSAR 技术实现面上控制和城市全覆盖，对于单个点上的连续快速沉降，还需进一步利用其它地面观测的方式开展工作。

为了保证本研究工作的连续性和发展性，下一步工作中拟采用RADARSAT-2 卫星获取的精细模式SAR 数据开展相应的监测工作。该数据分辨率为8 m，优于ENVISAT 数据，覆盖范围为50 ×50 km2，获取周期为24 d，应用PSInSAR时序分析方法，连续获取数据既可满足矿区动态监测，也可满足项目全区缓慢沉降监测的需要。

［1］ 刘国祥．利用雷达干涉技术监测区域地表形变［M］．北京:测绘出版社，2006:148．

［2］ 舒宁．雷达影像干涉测量原理［M］．武汉:武汉大学出版社，2003:99．

［3］ 程滔，单新建． CR、PS 干涉测量联合解算算法研究［J］． 地震，2007，27(2):64 -71．

［4］ 范景辉，李梅，郭小方，等．基于PSInSAR 方法和ASAR 数据监测天津地面沉降的试验研究［J］．国土资源遥感，2007(4):23 -27．

［5］ 张洁．干涉雷达数据处理技术与应用［D］． 武汉:中国地质大学(武汉)，2005．

［6］ 韩春林．SAR 成像技术和SAR 图像处理算法研究［D］．成都:电子科技大学，2004．

［7］ 杨成生． 基于D-InSAR 技术的煤矿沉陷监测［D］． 西安:长安大学，2008．

［8］ 董玉森．干涉雷达中的永久散射体和交叉干涉技术的研究与应用［D］．武汉:中国地质大学(武汉)，2007．