济宁兖州地区地面沉降InSAR监测试验研究

第一章 绪论

1. 任务来源

项目名称：山东省重点地区矿山开采遥感动态监测

课题名称：济宁兖州地区地面沉降InSAR监测试验研究

2. 工作区范围

工作区主要包括济宁市、兖州市、曲阜市和邹城市所辖的主要煤矿区，区内煤矿分布广泛，主要为兖矿集团下属各煤矿。如图1.1所示，工作区经纬度如下：

北纬：35°16′00″～35°45′00

东经：116°26′00″～117°15′00″

图1.1 工作区覆盖范围图

3. 目标任务

围绕济宁和兖州地区煤矿开采、地下水开采等多种因素作用下地面沉降（陷）的分布与发展状况，开展以InSAR技术为主要手段的地面沉降调查与监测工作。

主要工作任务为兖州矿区所属各煤矿开采沉陷区及非煤矿沉降区的沉降调查与监测试验，利用ENVISAT卫星数据应用InSAR技术开展主要塌陷区的分布范围、演化趋势调查，试验煤矿开采区地面沉陷的动态变化过程及监测。为研究煤矿开采诱发地面沉降（陷）的变化规律及特征，探索越界开采的可能关系提供基础资料。

4．完成情况

围绕工作区地面沉降和开采沉陷监测，针对缓慢变形和快速沉降两种不同的形变类型，分别应用D－InSAR技术和永久散射体InSAR技术开展了工作区4000平方公里范围2004年至2010年地表形变监测，获取了多期动态监测速率和累积沉降量结果，查明了工作区地面沉降（陷）的影响范围和分布区域。同时开展2012年工作区开采沉陷动态监测，获取了当前沉降的范围和特征。本项工作全面完成了既定目标任务，到达了预期目标。

5．取得的主要成果

1）完成济宁、兖州和邹城三个地区4000多平方公里煤矿开采区地面沉降历史状况调查，获取了2004年以来多期InSAR监测地面沉降速率，得到了不同时段的地面沉降速率和总体累积地面沉降量，综合处理工作区ENVISAT ASAR数据22景。

2）研究并应用了常规D－InSAR技术和永久散射体（PS）InSAR技术获取了短时间间隔下矿区快速沉降动态监测成果和长时间条件下全区域地面沉降速率。利用2004年12月-2006年2月和2008年12月-2010年10月2个时段的18景ASAR数据，综合应用多种InSAR处理技术提取了2个时段的地面沉降速率，查明了整个工作区地面沉降的分布现状，确定了区域性地面沉降分布特征和影响范围。

3）根据2004年12月-2010年10月的数据获取了近6年工作区累积地面沉降量，完成了沉降影响范围的统计分析；确定主要沉降中心位置和沉降量。

4）根据煤矿区开采沉陷的特点，获取了2012年4个时间段的矿区各开采点地沉降范围和强度动态变化过程。

第二章 工作区地面沉降（陷）基本状况

1. 工作区主要地面沉降（塌陷）状况

工作区以济宁、兖州、曲阜和邹城为主，地面沉降主要分为两类，分别为地下采水引起的缓慢地面沉降以及煤矿开采引起的地表沉陷。前者为缓慢微小变化，主要发生于济宁市与兖州市城区以及城区周边，而后者则多发生于煤矿区，与采矿条件、煤层厚度与赋存条件，开采方式等因素相关。

济宁市地面沉降主要发生在城区，发现于1988 年，1988 年7 月至1989 年12月，沉降中心位于牛市、后城、人民公园一带，15 个月内最大沉降量为65 mm， 沉降速率52 mm/a，发生地面沉降面积为61. 7 km²，其中大于20mm 的面积为22. 6 km²；1992 年8 月，地面沉降量大于20mm 的面积已达48. 0 km²，沉降中心最大沉降量达109mm。至1998年5月，城区已形成以南岱庄—师范学校为中心的沉降漏斗， 漏斗中心最大沉降量达202.3 mm，累计沉降量大于60mm的面积已达79km²。

自上世纪70年代大规模开采以来，济宁、兖州、邹城地区采煤塌陷问题越来越严重，塌陷区主要分布在邹城市、济宁市、微山县、任城区境内，鱼台县、金乡县和汶上县境内也有少量分布，土地塌陷总面积达到13518.2km², 其中邹城市境内采煤塌陷面积最大，占总面积的31.42%，任城区采煤塌陷面积占总面积的23.65%，微山县采煤塌陷面积占总面积的19.14%，其它县市塌陷面积占总面积的25.79%。

2. 济宁、兖州地区采煤塌陷状况

（一）煤炭资源分布及开采状况

      济宁市煤炭资源极为丰富，境内含煤面积3920平方公里，占全市国土面积的37%，煤炭资源储量270亿吨。境内主要煤田有兖州煤田、济宁煤田、滕南煤田、滕北煤田（济宁部分）等。煤炭资源主要分布在湖东冲积平原及黄泛冲积平原和泗河平原，除南四湖水域外多为高产稳产粮区，涉及兖州、曲阜、邹城、任城、中区、微山、鱼台、金乡、嘉祥、汶上、梁山11个县市区的近百个乡镇。境内煤层主要为厚煤层，大部分厚度达8-12米，较薄煤层也在2-3米之间，而且煤质优良，可开采利用价值较高。煤炭分布区域人均土地资源量少，耕地垦殖率高、村庄密集、人口密度大，工业产值及经济效益高。自上世纪70年代进入大规模开采以来，济宁市煤炭资源开采已进入全面开发和迅速发展期。境内现有兖矿集团、淄矿集团、枣矿集团、临矿集团、肥矿集团等10多家采煤企业。截至2009年底，济宁市现有生产矿井53对、在建矿井10对，设计生产能力8769万吨。2009年，济宁市原煤产量达到8115万吨，首次突破八千万吨，占到全省的56.4%，居全国第三位，煤炭资源的开发利用已成为济宁经济社会发展的支柱产业，在全市国民经济中有着举足轻重的地位。

（二）采煤塌陷地面积与分布

从20世90年代，随着济宁市煤炭开采规模的迅速扩大，土地塌陷逐渐成为严重的问题。目前，济宁市采煤塌陷地主要分布在邹城市、微山县、兖州市、曲阜市、任城区、鱼台县、汶上县和金乡县的28个乡镇， 308个村庄。

（三）采煤塌陷地类型

塌陷地的形成一般是在采空后3至5年逐步形成并稳沉的。济宁市地下煤层大部分是两层或两层以上，须在分层开采完后才能形成最后的稳沉。就采煤塌陷地而言，全市是采煤塌陷地形成的初期，当煤炭生产达到壮年期时，大面积塌陷会达到稳沉，原来的地质环境和地形地貌将改变。

塌陷地类型划分是认识塌陷地的重要手段，由于塌陷地类型地域差异明显，目前国际、国内还没有统一的划分标准和体系，依不同参照出现了多种划分方法。作者认为，从塌陷区治理的角度出发，应以塌陷地稳定程度、沉陷强度、塌陷地性质和形态(积水状况、土壤性质等)作为划分的主要指标，并要结合当地潜水位予以确立。据此，济宁市采煤塌陷地类型划分如表2。

表2　济宁市采煤塌陷地类型划分

截止2005年底，济宁市因采煤塌陷的22万亩土地中，浅层塌陷地占7. 9万亩，中层塌陷地占1. 6万亩，深层塌陷地占12. 5万亩。深层塌陷在济宁市表现较为突出，在全国具有代表性和典型性，特别是深层塌陷处于粮食高产稳产平原区，复垦难度大，对当地农业生产造成了严重破坏。

3. 济宁、兖州地区采煤塌陷区变化趋势

    2000-2005年采煤塌陷地呈增长趋势，增长速度各年份不同，大致成减缓趋势。出现这一现象的原因有以下几点:一是各年份采煤量不同；二是部分采矿企业开始回采第二层煤甚至是更多层煤，导致原来塌陷地区的多次沉降，下沉值增大，但面积扩张不明显；三是各县市积极开展了土地复垦活动，抑制了塌陷地增长速度。

在煤炭资源丰富的任城区、兖州市、邹城市金三角以及微山县，地面塌陷较突出，一些地方集中连片。

(1)邹城市境内采煤塌陷面积最大，占总塌陷面积的31.42%，主要位于邹城西部。兖州市采煤塌陷面积占15.22%，主要位于兖州南部。塌陷始于20世纪80年代，1995年进入高峰期，目前各矿均处于开采的后期，逐渐进入塌陷减缓期。

(2)任城区境内煤矿尽管开采历史不长，但塌陷面积却较多，占总塌陷面积的23. 65%。任城区境内均属新建矿井，塌陷在2000年开始发现，逐渐进入高发期，今后将是整个济宁境内的重点塌陷区。

(3)微山县是济宁市境内开采煤矿最多、开采历史最长的县，但其塌陷地面积仅占总塌陷面积的19.14%，原因是采煤主要位于微山湖下，对陆地破坏反而少。微山县西部及南部塌陷区最早发现于20世纪80年代初，80年代末至今为塌陷高峰期，随着境内资源的储量减少，逐渐进入塌陷减缓期。

这三个区域内下覆煤层厚达7-9m，甚至10m以上，地面塌陷深度最大可达8m，再加上该区地下水位较高，因此，易出现大面积的积水区，造成耕地绝产。

除此之外，鱼台县境内有鹿洼煤矿开采造成的塌陷区，金乡县金桥煤矿、汶上县唐阳煤矿开采对地面的影响也开始显现。但由于地下赋存煤层较薄，一般地面没有积水。

第三章 工矿区地面沉降（陷）InSAR监测技术

1. InSAR基本原理

InSAR（SAR Interferometry或Interferometric SAR），即合成孔径雷达干涉测量，是指利用SAR数据中的相位信息进行干涉测量处理，结合雷达参数和卫星位置信息反演地表三维及其变化信息的遥感手段。D-InSAR即差分合成孔径雷达干涉测量，是指对干涉相位进行差分处理提取变化信息的干涉测量手段，是InSAR技术的延伸。InSAR技术已被广泛用于大区域地形测图和地表信息提取。常见的InSAR方式有顺轨（Along track），交轨（Across track）和重复轨（Repeat pass）。由于地表形变监测是测量一定时间间隔内的地表形变量，因而重复轨是目前应用的主要方式，具体是指利用同一地区不同时刻获取的SAR数据进行差分干涉处理以提取地表形变信息。

2. D-InSAR基本数据处理方法

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

  (一) 干涉像对配准

在进行SAR干涉测量时，SAR图像必须精配准，从而保证干涉条纹具有良好的相干性。由于两幅图像的数据几乎来自空间的同一位置（星载SAR复数据，基线长度与天线到地表距离相比很小），两幅图像的相干像元主要的不同在于存在一定的偏移、小范围的拉伸及方位向轻微的旋转。图像配准一般经过粗配准，像元级配准、子像元级配准。影像配准的方法有基于卫星轨道参数的复图像对配准和利用x和y方向偏移量为变量，以相关系数作为匹配质量评价指标来进行粗匹配和精匹配。经过配准以后的影像对进行复共轭相乘就可以得到干涉纹图。

（二）地形相位消除

由干涉相位分析可以发现，要解算反映地形的地形相位，必须消除平地效应和地形的影响。平地效应消除的方法有：①基于轨道参数和成像区域中心点的大地经纬度计算平地效应；②根据图像能量计算平地效应；③通过测量距离向和方位向的占优势的干涉条纹频率来计算平地相位，然后消除平地效应。此后，根据已有DEM数据，结合卫星参数模拟地形相位，从干涉相位中减去该相位分量，由此得到差分相位图。   

（三）相干图生成

相干系数反映了相干性的高低，是指导相位解缠的一个重要指标。类似于幅度影像的相关系数经常用于影像匹配的测度，复影像之间用相干系数来衡量相似程度和干涉纹图的质量。相干性测度可以看作实数影像之间相关系数在复数域的推广。经过对应像元的相干性估算，得到了反映相干系数的相干图。

（四）噪声滤除

选择合理的滤波器来消除因大气波动所造成引起的相位变化。

（五）相位解缠

主要是将干涉相位进行解缠，得到能反映地形变化的地形相位。

（六）形变值计算和地理编码

    利用解缠相位和干涉测量的几何参数，进行差分相位到形变量的转换，将转换后形变量结果转换到特定坐标系下，生成形变图。

3. 地表形变D-InSAR监测技术

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

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

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

4. 永久散射体干涉测量（PSInSAR）技术及发展

永久散射体干涉测量（Permanent scatterer Interferometry，PSInSAR）技术的提出将InSAR技术向前推进了很大一步。PSI的核心思想是对永久散射体干涉相位进行时间序列分析，根据各相位分量的时空特征，估算大气波动，DEM误差以及噪声等，将其从差分干涉相位中逐个分离，最终获取每个PS的线性和非线形形变速率、大气延迟量（Atmosphere Phase Screen）以及DEM误差。经PS方法处理，获取的年度形变率的精度可以达到毫米级。与PSINSAR其具有相同处理思想的方法有CPT（Coherent Point Target）、IPTA（Interferometric Point Target Analysis）等。这些方法的共同点是基于大量的SAR数据（一般大于20甚至30景），从中筛选出具有稳定散射特性的相干点目标，构成离散点观测网络（较之常规的变形监测网密度更高），通过分析PS点目标相位变化获取地表形变状况。由于将永久散射体作为观测对象，降低了空间基线对相干性的影响，即使在临界基线的条件下，仍然可以通过分析PS差分干涉相位的变化反演形变信息。但该方法往往需要反映地表形变特征的先验模型，如线性形变速率模型。另外，为了提高散射体高程的估算精度，并进行大气校正，需要大量的SAR数据进行统计分析。

PS技术一般采用的线性形变模型提取点目标对应的形变量，如测量长时间下保持稳定移动速率的地表移动的现象。该方法的优点是能一次性的获取中尺度（约2000km²）范围内的地表形变信息。由于非线性形变可以用线性形变模型来模拟，因而一些非线性形变也可以通过线性形变测量得到。若观测对象表现出明显的非线性特征，并且形变量变化大，则在PS点目标覆盖的范围内出现了不连续的区域，产生不连续（空间和时间上的）原因是由于形变本身超出了所采用的模型的边界条件。这种情况下，若利用基于线性模型估算的形变速率来反演一定时间内的形变量，则必将与实际情况相差较大远。可以通过两种方法来弥补线性模型模拟非线性形变的不足，一种是采用非线性形变模型，另外一种是将长时间间隔分解为数个短时间段，利用函数模型模拟各个时间段内的形变量，进而求解非线性量。非线性模拟的处理过程相当复杂，而且非常耗时，限制了其用于大面积的形变测量，但随着处理技术的进步，处理时间将逐步缩小，处理的范围也可以进一步扩大。

PSInSAR对数据量要求较高，只有SAR影像个数达到一定的程度才能筛选出在整个时间跨度内具有稳定信号的PS点，对于相干目标较多的地区，如城区等，由于地物在长时间间隔下保持了很高的相干性，数据量较少的情况下也可进行处理。由于PSInSAR基于离散点相位解缠（时域或二维平面），当点的密度不能达到形变量对相位解缠的要求时，则很难得到准确的结果。一般情况下，至少必须满足5个/km² 才能完成相位解缠。对点密度要求的另一个原因是大气问题，鉴于大气的空间相关距离，因而要求PS间距尽可能小于1km，以减小大气的影响。

与序列干涉纹图方法比较而言，PSInSAR方法可以利用长基线距的SAR数据进行形变参数估计，在DEM误差、大气相位估计以及非线性形变等的解算上有明确的阐述。此外，通过对残余相位的处理，可以进一步提高估计精度和PS数量，在城市地面沉降监测，甚至重大基础设施的监测方面将发挥更大作用。尽管表现形式和处理方法不尽相同，但PS方法的一个共同点是基于大量的SAR数据（15-25景甚至更多），从中筛选出具有稳定散射特性的点目标（即PS，永久性散射体），构成离散点观测网络（较之常规的变形监测网密度更高），通过分析PS点目标不同时期的相位变化获取地表形变情况。PS点目标可以分为两种，一种是自然因素形成的，另一种是人为因素形成的，即使时间间隔长达数年，它们仍然具有稳定可靠的雷达后向散射特性。PS实地对应的地物类型为楼房角，房顶，路灯，井盖，裸露岩体以及其他在长时间条件下仍具有稳定可靠散射特性的散射体等。尽管在PS方法进行数据处理前无法得知某一地区存在的PS点的位置及密度等情况，但可以肯定的是城区PS点目标分布较多，一般在200-600个/km²。可见，PS方法更适合城市地表形变测量。

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

5. 主要雷达数据及参数 

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

表1所列为当前几种雷达传感器所提供的SAR数据。从雷达波段而言，目前主要的中等分辨率数据为C波段数据（ERS、ENVISAT和RADARASAT-1/2）。由于这些卫星具有较高的定位精度（可以达到5cm），空间基线解算结果可靠，数据相干性较好，因此，在多数地表形变监测中均可达到很好的监测结果。除此之外，其他波段的SAR传感器主要是ALOS、CosmoSkymed以及 TerraSAR-X。ALOS所搭载的SAR传感器为L波段，属于JERS的后续，所获数据的分辨率较前有所提高，重访周期为46天。

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

结合本次工作的要求，为了满足对工作区开采沉陷历史状况和现状的调查，满足一景图像对全工作区的覆盖，选择利用ENVISAT卫星接收的存档数据作为数据源（其它卫星无存档数据）。该数据轨道定位精度高、分辨率适中，其重放周期和覆盖范围均满足本次调查工作的目标要求。

第四章 主要工作内容与关键技术

1. 工作区地表形变InSAR监测内容

工作地区为典型的煤矿开采区，开采沉陷诱发的矿区环境破坏不仅出现在井田开采上覆地区，同时也影响到了城市和城镇的地面稳定性。比较而言，采区上覆地区的沉降强度和范围远大于沉陷盆地边缘地区，但由于上覆地区多为农田，引发的直接破坏程度较低；而矿业城市往往位于沉陷盆地边缘，受开采沉陷影响较小，由于城镇人口集中、建筑物密集，由此所造成的破坏程度往往表现明显。因而，矿山沉陷InSAR监测的首要任务是明确开采沉陷的变化特征，进而选择合理的数据处理方法和相关技术参数。

（1）开采沉陷InSAR动态监测

综合分析开采沉陷的基本特点和InSAR技术应用的主要条件，矿山开采沉陷InSAR监测需选择合理InSAR系统参数。包括雷达数据获取的周期、分辨率和波长大小。开采沉陷快速活动变化产生高的相位梯度，当相位梯度小于临界梯度时，即未发生非相干移动的前提下，对于同一种雷达数据，利用短基线（包括时间基线和空间基线）干涉处理方法可以实现开采沉陷动态监测。此外，相同的波长情况下，适宜选择高分辨率雷达数据，以此降低相位梯度，保证相位解缠的可靠性。

（2）矿业城市缓慢地面沉降监测

济宁、兖州等城市不仅受采空塌陷的影响，同时也受到地下水开采的影响。其中，济宁市区受以往老采空区缓慢沉降影响突出。这种缓慢沉降一般持续2-4年，甚至更长时间。连续缓慢沉降将直接影响城市建筑物的稳定性，特别老采空区上部开展开发和建设时。由于城市存在若干数量的相干目标，因而可利用永久散射体（或相干目标InSAR时间序列分析）技术实现区域性缓慢沉降的连续监测。

2. 工作任务

根据工作地区地面沉降发展变化基本状况，结合卫星数据获取情况，本项工作取得的主要成果分为2项：

第一项开展了济宁、兖州、曲阜和邹城地区地面沉降基本状况调查。本项工作的主要目标在于查明全区地面沉降分布范围和沉降强度，主要利用2004年12月至2006年2月和2008年12月至2010年10月间的ENVISAT雷达数据开展兖州、济宁等4个地区区域性地面沉降发育和发展特征调查，查明了2个时段的地面沉降速率和2004年至2010年6年的累积沉降量，编制相应的成果图件，查明了目前全区地面沉降的发育特征。

第二项主要针对矿山开采沉陷动态监测，利用相干性良好的季节获取的雷达数据，开展了多期矿山开采沉陷动态监测，由干涉图的变化状况反映了矿区开采沉陷的范围，编制相应的成果图件。

3. 工作方法

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

（一）工作区SAR数据获取

本项工作采用的主要雷达数据为欧空局自2004年12月至2010年10月间和2011年11月－2012年2月间接收的ASAR数据，其有效工作波段是C波段，数据获取周期为35天和30天（2010年10月轨道调整后周期变为30天），分辨率为20m。

（二）数据处理和分析

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

（三）地表形变信息提取

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

（四）InSAR监测结果后处理

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

4. 关键技术方法

针对兖州和济宁市区地面沉降等缓慢微小地表形变现象的监测，工作中采用永久散射体干涉测量技术（PSInSAR）。永久散射体干涉测量技术的核心是对于相干点目标干涉相位的时序分析，其出发点是为解决传统InSAR技术在低相干条件下地表形变信息的准确提取，降低大气波动等误差对监测结果的影响，通过利用累积SAR影像进行干涉相位时间序列分析。常规的时序分析方法利用相干系数作为相干目标识别的标准，由于受制于相干系数估算窗口的影响，噪声会引入到估计的相位中去，同时，滤波过程中也会降低干涉相位的真实性。因此，就产生了利用相干点目标InSAR技术，其基本思想是利用大量连续观测获取的SAR数据开展相干点目标的时序分析，提取地表形变速率和累积形变量（图4.2）。

对于干涉像对的选择，在数据量较少（小于10景）的情况下，选择基线较短（一般小于300m）的干涉图像，进行短基线干涉像对组合，增大单一差分干涉纹图可获取形变信息的时间跨度，尽可能多的获取观测值的个数，以此来求解形变速率。这就是短基线InSAR技术的核心思想。短基线InSAR技术中为恢复形变过程，将差分干涉相位视为观测值，而相干目标的形变序列为未知数，利用参数估计的方法求解每一相干目标的形变序列值。

PSInSAR技术的提出将InSAR技术向前推进了很大的一步，其核心思想是对永久散射体（Permanent Scatter）的干涉相位进行时间序列分析，根据各相位分量的时空特征，估算大气波动，DEM误差以及噪声等，将其从差分干涉相位中逐个分离，最终获取每个PS的线性和非线形形变速率、大气延迟量（Atmosphere Phase Screen）以及DEM误差。经PS方法处理，获取的年度形变率的精度可以达到毫米级。

PSInSAR技术一般采用的线性形变模型提取点目标对应的形变量，如测量长时间下保持稳定移动速率的地表移动的现象。该方法的优点是能一次性的获取中尺度（约2000km²）范围内的地表形变信息。由于非线性形变可以用线性形变模型来模拟，因而一些非线性形变也可以通过线性形变测量得到。若观测对象表现出明显的非线性特征，并且形变量变化大，则在PS点目标覆盖的范围内出现了不连续的区域，产生不连续（空间和时间上的）原因是由于形变本身超出了所采用的模型的边界条件。这种情况下，若利用基于线性模型估算的形变速率来反演一定时间内的形变量，则必将与实际情况相差较大远。可以通过两种方法来弥补线性模型模拟非线性形变的不足，一种是采用非线性形变模型，另外一种是将长时间间隔分解为数个短时间段，利用函数模型模拟各个时间段内的形变量，进而求解非线性量。非线性模拟的处理过程相当复杂，而且非常耗时，限制了其用于大面积的形变测量，但随着处理技术的进步，处理时间将逐步缩小，处理的范围也可以进一步扩大。

PSInSAR对数据量要求较高，只有SAR影像个数达到一定的程度才能筛选出在整个时间跨度内具有稳定信号的PS点，对于相干目标较多的地区，如城区等，由于地物在长时间间隔下保持了很高的相干性，数据量较少的情况下也可进行处理。由于PSI基于离散点相位解缠，当点的密度不能达到形变量对相位解缠的要求时，则很难得到准确的结果。一般情况下，至少必须满足5个/km² 才能完成相位解缠。对点密度要求的另一个原因是大气问题，鉴于大气的空间相关距离，因而要求PS间距尽可能小于1km，以减小大气的影响。

与序列干涉纹图方法比较而言，PSInSAR方法可以利用长基线距的SAR数据进行形变参数估计，在DEM误差、大气相位估计以及非线性形变等的解算上有明确的阐述。此外，通过对残余相位的处理，能进一步提高估计精度和PS数量，在城市地面沉降监测，甚至重大基础设施的监测方面将发挥更大作用。

5. 工作中所用的SAR数据

用于干涉测量的SAR主要为ENVISAT ASAR标准模式数据。根据本项工作的主要任务和工作区范围，该数据的空间分辨率为20m，能满足大范围地面沉降详细调查与矿区沉陷动态监测需要，是目前常用且应用效果显著的SAR数据。获取ASAR数据共计22景。

区域性地面沉降监测所使用的数据从2004年12月到2010年10月，计18景。

第五章 区域地面沉降与开采沉陷InSAR监测

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

工作中共接收ENVISAT ASAR数据22景，ENVISAT卫星变轨前18景，由于该卫星于2010年10月后轨道调整，在变轨后从2011年年底截止2012年3月接受数据4景。考虑到矿山塌陷的大变形过程对相干性的影响，以及季节性植被覆盖变化，因而选择冬季干涉像对，且间隔时间较短的数据进行D－InSAR数据处理，共获取变轨前8个能清晰反应矿山沉陷的干涉纹图，用于反映工作区矿区开采沉陷2010年之前历史状况的动态过程。

监测成果揭示了济宁、兖州、曲阜和邹城地区各个煤矿开采区在2004年11月到2005年3月、2005年12月到2006年2月、2009年12月到2010年2月以及2011年11月到2012年2月四时间段的地面塌陷基本情况，表明了采煤塌陷的快速变化过程。干涉条纹（每个干涉条纹色周表示3.1cm的沉陷变化）的密集程度反映了沉陷幅度的大小，按照干涉条纹图的空间分布特征可以确定沉降的幅度，并基本判别地下煤矿的开采方向。

全区总的开采沉陷区主要分布于兖州、曲阜和邹城沿市泗河沿岸以及济宁市区周边。这些的地区不同时段内差分干涉条纹密集，表明了各个煤矿在监测时段内开采沉陷的发生位置和范围。兖州市和曲阜市泗水河沿岸主要开采沉陷范围分布在张家村、小孔村、兴隆庄镇、中心店镇、平阳寺镇和北宿镇等地区，由雷达图像上发现这些地区沉降塌陷集水坑密布，因而由于沉降引起的干涉条纹密集。济宁市周边四个方向均有煤矿分布，开采沉陷的范围比较密集，但市区内部没有塌陷引起的干涉条纹变化。差分干涉条纹在济宁市区东北部分布于柳巷镇的南营，南部分布于前十里营、后里庙和南阳湖农场一带。市区西北部分布较为密集，分别为梁庄、天宝寺、二十里铺镇、南张镇和刘堤头。汶上县的石村和豆村也有开采沉陷出现。

在2004年11月到2005年2月3个月的时间内，共有4个差分干涉纹图（图5.1a、b、c、d）显示了该时间段的塌陷变化过程，各个采区的干涉条纹密度随着时间间隔的增加而整体上增大。由041121-050306差分干涉纹图（图5.1b）可知兖州市城区东部与曲阜交界（泗河北）的开采沉陷漏斗干涉条纹达5个色周， 沉陷量达到15cm，表明在105天内由于煤矿开采引起了最大15cm的沉降，沉降沉陷出盆地形状的分布特征，与传统开采沉陷的沉陷盆地分布相同。

位于济宁南张镇西北的沉陷漏斗在041121-050130（图5.1a）的干涉纹图超过2个色周，沉陷量超过6cm。050103和050306两景SAR数据相干性好，差分干涉纹图（图5.1d）各沉陷中心清晰，可以明显看到该沉陷漏斗中心向北偏移。

2005年12月到2006年2月有1个相干性好的干涉纹图（图5.1e），该干涉纹图表现了该时间段内各沉陷中心的沉陷状况，曲阜与兖州在泗河交界处的沉陷漏斗干涉条纹有5-6个色周，中心沉陷量达15-18cm；济宁南张镇-刘堤头间的沉陷漏斗条纹达4个色周，中心沉陷量超过12cm。

2009年12月到2010年2月相干性良好的3个干涉纹图（图5.1f、g、h）清晰地表现了该时间段内的开采沉陷状况。由091220-100124（图5.1f）和100124-100228（图5.1g）干涉纹图可知，位于曲阜与兖州在泗河交界处张家村煤矿一带的沉陷范围在35天的时间间隔内达到3个色周的沉降变化，表明在091220-100228共70天的时间内沉陷量达到近18cm，属于快速沉降阶段。

济宁主城区西北前王村南沉陷中心在091220-100124和100124-100228干涉纹图中分别达到了6个和3个以上色周，共计约10个条纹，表明最大沉陷量超过30cm。图5.1h为直接利用091220和100228两景数据干涉后的条纹图，反映了105天内整个沉降区的变化状况，是图5.1f和图5.1g反映的沉降量的总和。总体上，干涉图的空间位置和分布形态一致，但局部干涉图的条纹清晰程度降低。这种结果产生的原因是由于时间间隔进一步增大，相干性下降引起的。同时由于沉陷量增大，局部出现非相干的移动变化，这种快速变化限制了利用InSAR干涉图之间确定沉降量的精准程度。同时，位于楼庄和梁家营一带出现了局部回弹，干涉条纹与其它沉降引起的形态相反，变化量为1个色周，影响范围约2-3km。

2. 煤矿区开采沉陷2011年-2012年现状InSAR动态监测

2011年以后获取的数据为降轨后的ENVISAT ASAR数据，共4景，重访周期30天，在工作区范围内仍保持相干性，可用于矿区开采沉陷动态监测。与变轨之前的数据相比较，相干性较前较为明显，主要原因是空间基线长度略有增大。

利用2011年11月至2012年2月的获取的4景ENVISAT数据进行组合，根据时间间隔依次为111124－111224，111224－120123、120123－120222，进行差分干涉处理得到了3个时间间隔均为30天的干涉纹图（图5.2a、b、c）。

从111124-111224干涉纹图（图5.2a）可知：曲阜张家村和小孔村一带干涉条纹色周达1－2周，30天的时间内沉陷量为3cm以上，在之前（图5.1f）该地区沉陷漏斗中心位于西侧与兖州交界处。

从120123-120222干涉纹图可以清晰地看到济宁南张镇一带沉陷中心干涉条纹色周超过4个，60天的沉陷量超过12cm，而小孔村沉降中心色周为1个，沉陷量为3cm左右，可见该沉陷中心在2011年11月到12月沉陷较快。

通过对3个时间段的比较，可以看出各个矿区开采沉陷的时空变化特征，及各个时间段内随时间变化，开采沉陷的发展方向，一些老的开采沉陷区由于积水区域失去相干，但仍可监测其沉陷范围，如兴隆庄煤矿、邹城市辖西北与兖州交界处、微山湖周边等，依然保持着较快的地面沉陷速率。

   总体上，2012年初工作区内主要的开采沉陷区仍分布于济宁市区周边、兖州与曲阜交界沿泗河一线。由InSAR监测到的快速沉降变化为开采沉陷动态过程，间接反映了各个煤矿的采掘范围。

   对于InSAR形变监测而言，煤矿区开采沉陷属于非线性变化过程，一般需经历初始沉降阶段、快速沉降阶段和稳沉阶段。在相干性良好的情况下，目前可以利用缩短数据获取的时间间隔，即利用短周期的SAR数据，可很好的监测各个采区的沉降过程。但实际应用中，受制于采空区上覆植被变化的影响，相干性往往随着季节的变化而降低，影响了对每个煤矿开采区的动态连续监测。上述监测结果均为冬季获取，干涉条纹总体上非常清晰，反映的矿区开采沉陷状况真实可靠。由于矿区会出现很大的塌陷，随着塌陷的积累和水位的上升，出现集水坑，这影响到了直接测量累积塌陷量，不仅仅是InSAR手段。因而，对于开采沉陷对于城区和城市的影响，仍需要采用相干目标连续监测的方法，将能够监测到的信息提取出来，确定整个采区的影响程度。

3. 全区2004年-2010年区域性地面沉降InSAR连续监测

尽管受制于相干性的降低，难以应用常规D－InSAR技术进行采空区快速沉陷的连续监测，但对于采空区相对较为缓慢的形变过程，仍可以利用长时间序列进行PSInSAR处理，提取相干目标的沉降速率。这对于研究矿业城市地面形变具有重要意义。本次工作区内的济宁、兖州和邹城等城市，地面沉降主要分为缓慢沉降以及煤矿开采沉陷。前者为缓慢微小变化，主要发生于济宁市城区以及城区周边，而后者则多发生于煤矿区，与采矿条件、煤层厚度与赋存条件，开采方式等因素相关。因而，在煤矿区开采沉陷动态监测的基础上，实现全区地面沉降的连续监测。

4. 全区累积地面沉降量

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

通过对累积地面沉降量空间插值分析得到该区域累积沉降等值线图（图5.5和图5.6），济宁、兖州、曲阜和邹城四辖区约4060平方公里的面积，去除大面积水域共约3922平方公里，分析统计可得各沉降范围影响面积如表5.3所示，各占比例如图5.7所示。

对济宁、兖州、曲阜、邹城各沉降范围影响面积进行单独统计，各地区主要沉降主要分布如表5.4－5.7所示

5. 矿区范围与沉降范围比较

        为比较开采沉陷影响的实际范围，将济宁和兖州地区煤矿区矿权范围数据与实际开采沉陷范围进行叠加，得到图5.8。图中蓝色范围为主要煤矿区矿权划定范围，比较发现，两个主要沉降区分别对应着兖州和济宁地区两个主要矿区范围。

      从沉降发生范围上，济宁兖州两地区的沉降范围均对应于煤矿矿权范围，沉降漏斗的形态上与矿权范围均有非常好的一致性。其它无矿权范围地区均为较小沉降区。

对比沉降量的大小可知，兖州地区大于60mm沉降量的地区为煤矿开采区。从北向南依次对应着单家村煤矿、古城煤矿、星村煤矿、兴隆庄煤矿、杨村煤矿、鲍店煤矿和横河煤矿。沉降最大处位于兴隆庄煤矿。整个沉降量为60mm等值线所包围的区域均为煤矿矿权覆盖范围。

     济宁地区开采沉陷范围主要集中在葛亭煤矿、运河煤矿、岱庄煤矿、唐口煤矿和许厂煤矿。沉降量为60mm以上的等值线均为上述矿区。济宁南部的东山王楼煤矿沉降突出，中心区大于80mm，影响范围较广。

第六章 总结与进一步的工作

1．工作总结

本次工作针对兖州和邹城地区区域性地面沉降和煤矿开采区快速沉陷两种类型地表形变，利用2004年－2012年多期EVSAT ASAR雷达数据，应用PSInSAR处理技术与常规D－InSAR差分干涉处理方法，获取了全区2个不同时间段的区域性地面沉降速率和6年时间内全区累积沉降量结果，并以冬季相干性良好的差分干涉纹图实现了开采沉陷动态过程监测，得到了多个时段矿区动态沉陷监测成果。监测结果查明了济宁、兖州、曲阜、邹城四市地面沉降分布状况，统计得到了全区累积沉降量和影响地区范围，全面掌握了四个地区目前地面沉降的总体分布态势和局部沉降特征。

区域性地面沉降结果表明，济宁、兖州、曲阜、邹城的地面沉降主要与矿山开采相关，城市活动等其他因素带来的地面沉降甚微，工作区内地面沉降显著区有两个，一个为兖州市东南、曲阜市西部以及邹城市西北部沉降区；另一个明显的沉降区域为济宁市任城区北部。沉降区内2004年12月至2006年2月最大年沉降速率达到50mm/a，2008年9月到2010年10月最大年沉降速率超过40mm/a。比较两期监测成果，虽然最大沉降速率有所减小，但沉降严重区域（年沉降速率>20mm/a）有扩大的趋势。

利用常规D-InSAR技术实现了矿区开采沉陷的动态监测。根据ENVISAT卫星从2004年至2012年间获取的多景雷达数据，获取了不同时期工作区内各个煤矿地面开采引起的地面变形状况，差分干涉图中干涉条纹密度的增加表明了开采沉陷的动态积累过程。考虑到矿山塌陷大变形以及季节性植被覆盖变化对相干性的影响，工作中主要以冬季获取的数据进行D－InSAR数据处理，共获取2010年之前8期，2011-2012年3期监测成果，反映了工作区内各矿区开采沉陷2010的动态过程。区内沉煤矿如古城煤矿、兴隆庄煤矿、杨村煤矿、东滩煤矿、鲍店煤矿、许厂煤矿、葛亭煤矿、运河煤矿等监测效果显著。2011-2012年的结果显示沿泗河分布的兴隆庄煤矿、南屯煤矿等，地面沉陷依然较快。

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

2．进一步的工作

由于ENVISAT 卫星轨道调整，2011年11月至2012年2月获取的ASAR数据相干性较之前差。为了保证本项工作的连续性，需选择其他雷达数据作为替代数据源。建议下一步工作中采用RADARSAT－2卫星获取的精细模式SAR数据开展相应的监测工作（图6.1所示）。该数据分辨率为8m，优于ENVISAT数据，覆盖范围为50*50km²，获取周期为24天，应用PSInSAR时序分析方法，依连续获取数据既可满足矿区动态监测，也可满足项目全区缓慢沉降监测的需要。