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误差。
这个估计是通过调整一个模型函数来实现的，该模型函数只针对那些显示出比较好的干涉时间相干性的像素而建。