DInSAR数据的地面沉降研究
——西班牙穆尔西亚市土体参数识别和地面沉降预
测
[西班牙] R. Tomás, G. Herrera b, J. Delgado a等
田 芳 译；冯翠娥、段 琦 校译
地面沉降是世界上很多地区都遭受的一种灾害，每年都会引起巨大的经济损失。穆尔西亚市（Murcia）位于西班牙东南部，从上世纪90年代就受到地面沉降的影响。本文利用永久散射体干涉测量（PSI）技术对含水层过度开采引起的地面沉降进行了遥感监测。特别是，相关像元分析技术（Coherent Pixels Technique，简称CPT）已经应用到ERS和ENVISAT卫星获取的SAR图像上。利用1993～1995年的CPT位移时间序列对提出的一维沉降模型进行参数识别。因此，CPT时间序列已经成功地用于获取土体的物理参数。然后利用该模型预测了1993-2007年间的变形。通过比较模型预测值与1995～2007年间实际的沉降时间序列，发现平均绝对误差为3.2±2.5 mm。尽管使用的这个一维模型是简化的，但是研究结果显示了CPT获取的位移信息在识别与验证因地下水过度开采引发的地面沉降数值模型中的有效性，能够用来预测未来水位下降时的含水层响应。
1 简 介
含水层过度开采引发的地面沉降灾害正影响着我们的社会。很明显的现象就是，大面积地区的地面在多年内持续出现几毫米甚至几米的垂向位移。沉降问题对于城市地区而言是非常重要的一个问题，因为它能够破坏基础设施，引发严重的经济损失。据估计，全世界有150多个城市出现了因地下水过量开采而引发的严重的地面沉降问题（Hu等，2004）。一些熟知的例子包括意大利波河平原（Po Valley），墨西哥城，美国羚羊谷、圣塔克拉拉谷（Santa Clara Valley）和圣约魁谷（San Joaquin Valley），泰国曼谷以及中国上海。在西班牙东南部的大城市，因地下水过量开采而引发的地面沉降已经造成5000多万欧元的损失，在1992~1995年旱期之后产生了严重的社会影响（Martínez等，2004）。
在过去的十年间，SAR差分干涉测量（DInSAR）已经成为一种重要的地面位移监测遥感技术，拥有低成本和毫米级精度等优点。最简单的DInSAR技术是基于一对SAR图像产生的单个干涉图（Rosen等，2000），以此来分析地下水开采引起的沉降的文献可见Galloway等（1998），Hoffmann等（2003）和Hoffman（2003）。DInSAR结果的一次显著改进要归功于永久散射体干涉测量（PSI）这组算法，PSI基于一大组SAR图像获得的多个干涉图的同时处理（Ferretti等，2000；Berardino等，2002；Mora等，2003；Arnaud等，2003；Werner等，2003；Hooper等，2004）。应用PSI进行沉降监测的例子可参考Ferretti等（2004），Tomás等（2005），Casu等（2006），Zerbini等（2007），Bell等（2008）和Herrera等的著作（2009）。
面对社会日益受地面沉降影响的现实，文献中出现了各种地面沉降的预测方法。按照Xu等（2008）的标准，可将预测方法分为五类：（a）统计方法，如影响函数、灰色理论和回归分析；（b）一维数值模型；（c）准三维渗流模型；（d）三维渗流模型；（e）三维全耦合模型。
统计方法利用地面沉降与其他因素之间的简单关系（灰色理论模型），地面沉降与开采量或者水位的关系（回归分析方法），或者简单地建立起沉降在时间上的变化趋势（影响函数）。一维数值方法认为影响含水层和弱透水层固结的水位变化只发生在垂向上。这种方法只利用垂向的土体参数计算某一个点上的沉降。准三维渗流和三维渗流模型利用太沙基一维固结方程计算沉降（Terzaghi和Fr?hlich，1936）。这两个方法之间的差别体现在含水层和弱透水层中水的渗流方向上。准三维渗流模型认为含水层中的渗流是水平的，弱透水层中的渗流是垂向的，而三维渗流模型假设水流的方向是三维的。三维全耦合模型利用Biot的3D固结理论，通过模拟3D水流来计算沉降（Biot，1941）。
选择一个最方便的沉降预测模型依赖于几个复杂的因素，还取决于不同地区的地质条件（Hu等，2002）。
利用上述这些方法进行沉降预测的文献可参考Gambolati（1975），Helm（1975），Helm（1976），Gambolati等（1991），Monjoie等（1992），Mizumura等（1994），Shearer等（1998），Gambolati等（2001），Larson等（2001），Hu等（2002），Chen等（2003），Zhou等（2003），Don等（2005），Ferronato等（2007），Shi等（2007），Shi等（2008），Xue等（2008）和Wu等（2009）的著作。
本文将利用称为相关像元分析技术（Mora等，2003；Blanco等，2008）的PSI方法与ERS和ENVISAT装载的SAR传感器获得的图像，重点分析发生在
西班牙东南部穆尔西亚市的地面沉降。前人的工作已经用原位数据示范和证实了不同PSI技术（包括CPT）在该地区的适用性（Tomás等，2005；Herrera等，2008；Tomás，2009）。本文将致力于建立一个简单的沉降模型，并将其作为一种预测工具。利用不同时期的干涉测量时间序列数据对该模型进行识别和验证，以显示其在预测地面沉降未来发展方面的潜力。
本文的结构如下：第二节介绍研究区的地理和地质背景；第三节介绍并简单评述一下干涉测量获得的沉降数据。接下来的第四节主要是介绍地面沉降预测方法，包括公式、常数的识别，以及利用时间序列数据进行验证。第五节是主要的结论。
2 地理和地质背景
塞古拉河（Segura River）的Vega Media（简称VMSR）位于Betic山脉东部的Bajo Segura
盆地（Montenat，1977）。VMSR的南边界由盆地基底岩石（二叠系~三叠系）构成，北边界由沉积岩构成（泥灰岩、砂岩和砾岩）。山谷里是来自于塞古拉河和Guadalentín河的新近沉积物，其中地表是全新统，地下一定深度处是更新统到上新统。
从岩土工程角度来看，这些新近沉积物的压缩性很高，问题最大。Rodríguez Jurado等人（2000）和Mulas等人（2003）对VMSR所有沉积物进行了岩土工程描述。他们的模型显示，出露在山谷边界处的同一种沉积岩也在山谷中的某一深度处发现，其特点是压缩性很低，甚至可以忽略。而它们上方的表层新近沉积物的压缩性为中等~高。
VMSR也是所谓的“Guadalentín–塞拉古河第四系含水层系统No. 47”的一部分（IGME, 1986）。该含水层可以划分为两个单元（Cerón和Pulido，1996；Aragón等，2004）：浅部单元，为微承压含水层或者弱透水层，由粉质粘土、粘质粉土与砂夹层构成；第二单元，即“深部承压含水层”，为细砂和砾石组成的多层含水层。浅部单元的水位在地下几米。由于沉积物中细粒成分的含量极高，所以这部分含水层的渗透系数较低，几乎没有被开采。因此，水位季节波动较小，最大为几米。
第二单元的水平和垂向渗透系数的变化范围分别是10~100 m/d和1~50 m/d（Aragón等， 2004）。第二单元是包含几个砂砾石层的含水层组。开采程度最高的位于这个单元的顶部，埋深大约为20m，该层的测压水位在地表以下几米，在干旱期由于过度开采出现了时间上的显著变化。特别是在1992~1995年间，这种现象特别明显，测压水位的峰值下降了15m。因此，VMSR
发生了大面积的地面沉降，造成了建筑物和大型公共设施的破坏（Mulas等，2003；Martínez等，2004）。Tomás等人（2005，2006），Tomás（2009）和Herrera等人（2008，2009）利用SAR差分干涉测量法测量了穆尔西亚大城市地区在干旱期间的地面沉降，探测到的最大位移是12cm。从水文地质学角度来看，第二单元的更深处存在其它砂砾石层，但是对它们的了解很少。
Mulas等人（2003）提出了一个刻画和模拟1992~1995年旱期发生的沉降的模型。根据这个模型，当水从深部含水层最顶部的砂砾石层中抽出的时候，就会发生沉降。垂向梯度的产生使得地下水从上部含水层（浅部单元）朝着砂砾石层向下流动，引起水位下降。因此，浅部含水层的孔隙水压力降低，含水层发生固结。
超孔隙水压力的消散速率是地面沉降计算中的一个关键问题。对5个未扰动的淤泥和粘土样品进行了渗透性的实验室测试，结果表明渗透性非常低且具有变化性（8.20×10?11 ~6.24×10?8 m/s）。因此，可以推断固结过程是缓慢的。在这种认识下，为了测量弱透水层和含水层中的测压水位，穆尔西亚市安装了分层水压计，结果显示了这两个层的测压水位变化非常相似，极有可能是因为易于排水的砂层的存在。
3 沉降数据
本文的地面沉降测量是通过称为相关像元分析技术（CPT）的PSI技术获得的。可以在Mora等人（2003）和Blanco等人（2008）的文献中找到关于该技术的详细描述，本文在此只做简单描述。
3.1 方法概述
用两幅SAR图像得到差分干涉测量相位（Δψint）的公式为（Hanssen，2001）： Δψint = Δψflat + Δψtopo + Δψmov + Δψatmos + Δψnoise （ 1）式中，Δψflat是不考虑地形条件下与测距差有关的扁平球体部分，Δψtopo是地形相位，Δψmov是发生在两幅SAR图像获取期间的沿着卫星视线向（LOS）测量到的地面形变产生的相位，Δψatmos是由于大气影响产生的相位，Δψnoise是残余噪声源。式（1）中的前两项可以通过分析获得。特别是，Δψflat是很容易知道的，Δψtopo可以利用外部的DEM中计算得到。
CPT算法（Mora等，2003）假设，与变形相关的相位组成（Δψmov）可以分成两个新的相位项，一个是由于线性变形（Δψlinear）产生的，即整个时段内恒定速率的变形，另一个是由于非线性变形产生的（Δψnon-linear）。因此，CPT的应用分为两个主要的步骤，即提取线性和非线性变形项。线性项的提取包括估计平均变形速度和DEM误差。这个估计是通过调整一个模型函数来实现
的，该模型函数只针对那些显示出比较好的干涉时间相干性的像素而建。非线性项是应用空间-时间过滤器去除大气影响并识别非线性变形的低分辨率和高分辨率部分来提取的。与非线性变形部分比较的时候，大气隔离是可能发生的，在某种程度上，是由于大气影响具有时间和空间上的不同行为。
关于CPT处理方法的其他说明如下：首先，通过对同一传感器获取的图像计算干涉图的不同子集（没有交叉干涉图），并利用奇异值分解法将变形结果联系起来，可参考Blanco等著作（2005, 2008），这样就有可能综合使用来自不同传感器（ERS和ENVISAT）的图像。第二，缺少某段时间内的有效图像（如1994年和2001年）是可以解决的，一种方法是增加干涉图的最大容许时间基线（例如3年），目的是将数据空白之前和之后获取的图像联系起来，另一种方法是对数据空白扩大时间过滤器的时间窗。用这种方法的话，即使缺少某些时间间隔的测量数据，但整个时间序列还是完整的。这种方法在大部分为线性变形的假设下是有效的。
3.2 数据集和详细处理过程
本文分析的SAR数据总共由81幅图像构成，包括ERS-1（6幅）、ERS-2（56幅）和ENVISAT-ASAR（19幅），图像获取时间从1993年4月~2007年3月。为了进行处理，提取了许多大约10×10 km的单元格，对应于穆尔西亚的城市地区。根据SAR图像对形成的所有可能的干涉图中，只选出了185幅干涉图进行CPT处理，借助于空间基线、时间基线和多普勒质心差形成的3D空间的Delauney三角剖分。选择被限定在那些空间基线和时间基线分别小于250m和1000天，多普勒质心差在800Hz以下的干涉图。用于抵偿干涉相位的地形部分的外部DEM的精度为25 m×25 m，属于IGN（National Cartographical Service of Spain，西班牙国家地图服务）的地图数字数据库E20。利用多层相干标准已经选择出了具有可用的地面形变信息的像素，该标准就是针对大部分干涉图只选择那些相干性超过某一阈值的像素（本例中，至少40%的干涉图集分别是0.6，0.5和0.4）。这种多层方法的优势是使低质量像素也能参与到沉降计算中，因为它们很好地利用了那些以前处理过的高质量像素（Blanco等，2006，2008）。相干性计算意味着干涉图的一种空间平均，也称之为“多视（multi-looking）”，降低了SAR图像原始分辨率。本次工作利用“方位角中的15个像素×距离中的3个像素”的多视操作，因此最终的地面分辨率为60 m×60 m。
3.3 DInSAR处理结果和验证
利用CPT技术获得了穆尔西亚市1993~2007年间的总变形（LOS位移）图，该图的平均密度为每平方公里76个相干像素。从图中可以发现，该市南部和东部的地面沉降量较大，整个研究期内大于5cm，而北部地区的情况较为稳定。
基于DInSAR数据的地面沉降研究
