基于DInSAR数据的地面沉降研究.docx

1、基于DInSAR数据的地面沉降研究基于DInSAR数据的地面沉降研究 西班牙穆尔西亚市土体参数识别和地面沉降预测西班牙 R. Toms, G. Herrera b, J. Delgado a等田 芳 译；冯翠娥、段 琦 校译地面沉降是世界上很多地区都遭受的一种灾害，每年都会引起巨大的经济损失。穆尔西亚市（Murcia）位于西班牙东南部，从上世纪90年代就受到地面沉降的影响。本文利用永久散射体干涉测量（PSI）技术对含水层过度开采引起的地面沉降进行了遥感监测。特别是，相关像元分析技术（Coherent Pixels Technique，简称CPT）已经应用到ERS和ENVISAT卫星获取的SAR

2、图像上。利用19931995年的CPT位移时间序列对提出的一维沉降模型进行参数识别。因此，CPT时间序列已经成功地用于获取土体的物理参数。然后利用该模型预测了1993-2007年间的变形。通过比较模型预测值与19952007年间实际的沉降时间序列，发现平均绝对误差为3.22.5 mm。尽管使用的这个一维模型是简化的，但是研究结果显示了CPT获取的位移信息在识别与验证因地下水过度开采引发的地面沉降数值模型中的有效性，能够用来预测未来水位下降时的含水层响应。1 简 介含水层过度开采引发的地面沉降灾害正影响着我们的社会。很明显的现象就是，大面积地区的地面在多年内持续出现几毫米甚至几米的垂向位移。沉降

3、问题对于城市地区而言是非常重要的一个问题，因为它能够破坏基础设施，引发严重的经济损失。据估计，全世界有150多个城市出现了因地下水过量开采而引发的严重的地面沉降问题（Hu等，2004）。一些熟知的例子包括意大利波河平原（Po Valley），墨西哥城，美国羚羊谷、圣塔克拉拉谷（Santa Clara Valley）和圣约魁谷（San Joaquin Valley），泰国曼谷以及中国上海。在西班牙东南部的大城市，因地下水过量开采而引发的地面沉降已经造成5000多万欧元的损失，在19921995年旱期之后产生了严重的社会影响（Martnez等，2004）。在过去的十年间，SAR差分干涉测量（DIn

4、SAR）已经成为一种重要的地面位移监测遥感技术，拥有低成本和毫米级精度等优点。最简单的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）。应用PS

5、I进行沉降监测的例子可参考Ferretti等（2004），Toms等（2005），Casu等（2006），Zerbini等（2007），Bell等（2008）和Herrera等的著作（2009）。面对社会日益受地面沉降影响的现实，文献中出现了各种地面沉降的预测方法。按照Xu等（2008）的标准，可将预测方法分为五类：（a）统计方法，如影响函数、灰色理论和回归分析；（b）一维数值模型；（c）准三维渗流模型；（d）三维渗流模型；（e）三维全耦合模型。统计方法利用地面沉降与其他因素之间的简单关系（灰色理论模型），地面沉降与开采量或者水位的关系（回归分析方法），或者简单地建立起沉降在时间上的变化趋势（

6、影响函数）。一维数值方法认为影响含水层和弱透水层固结的水位变化只发生在垂向上。这种方法只利用垂向的土体参数计算某一个点上的沉降。准三维渗流和三维渗流模型利用太沙基一维固结方程计算沉降（Terzaghi和Frhlich，1936）。这两个方法之间的差别体现在含水层和弱透水层中水的渗流方向上。准三维渗流模型认为含水层中的渗流是水平的，弱透水层中的渗流是垂向的，而三维渗流模型假设水流的方向是三维的。三维全耦合模型利用Biot的3D固结理论，通过模拟3D水流来计算沉降（Biot，1941）。选择一个最方便的沉降预测模型依赖于几个复杂的因素，还取决于不同地区的地质条件（Hu等，2002）。利用上述这些方

7、法进行沉降预测的文献可参考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方法与

8、ERS和ENVISAT装载的SAR传感器获得的图像，重点分析发生在西班牙东南部穆尔西亚市的地面沉降。前人的工作已经用原位数据示范和证实了不同PSI技术（包括CPT）在该地区的适用性（Toms等，2005；Herrera等，2008；Toms，2009）。本文将致力于建立一个简单的沉降模型，并将其作为一种预测工具。利用不同时期的干涉测量时间序列数据对该模型进行识别和验证，以显示其在预测地面沉降未来发展方面的潜力。本文的结构如下：第二节介绍研究区的地理和地质背景；第三节介绍并简单评述一下干涉测量获得的沉降数据。接下来的第四节主要是介绍地面沉降预测方法，包括公式、常数的识别，以及利用时间序列数据进行

9、验证。第五节是主要的结论。2 地理和地质背景塞古拉河（Segura River）的Vega Media（简称VMSR）位于Betic山脉东部的Bajo Segura盆地（Montenat，1977）。VMSR的南边界由盆地基底岩石（二叠系三叠系）构成，北边界由沉积岩构成（泥灰岩、砂岩和砾岩）。山谷里是来自于塞古拉河和Guadalentn河的新近沉积物，其中地表是全新统，地下一定深度处是更新统到上新统。从岩土工程角度来看，这些新近沉积物的压缩性很高，问题最大。Rodrguez Jurado等人（2000）和Mulas等人（2003）对VMSR所有沉积物进行了岩土工程描述。他们的模型显示，出露在山

10、谷边界处的同一种沉积岩也在山谷中的某一深度处发现，其特点是压缩性很低，甚至可以忽略。而它们上方的表层新近沉积物的压缩性为中等高。VMSR也是所谓的“Guadalentn塞拉古河第四系含水层系统No. 47”的一部分（IGME, 1986）。该含水层可以划分为两个单元（Cern和Pulido，1996；Aragn等，2004）：浅部单元，为微承压含水层或者弱透水层，由粉质粘土、粘质粉土与砂夹层构成；第二单元，即“深部承压含水层”，为细砂和砾石组成的多层含水层。浅部单元的水位在地下几米。由于沉积物中细粒成分的含量极高，所以这部分含水层的渗透系数较低，几乎没有被开采。因此，水位季节波动较小，最大为几

11、米。第二单元的水平和垂向渗透系数的变化范围分别是10100 m/d和150 m/d（Aragn等， 2004）。第二单元是包含几个砂砾石层的含水层组。开采程度最高的位于这个单元的顶部，埋深大约为20m，该层的测压水位在地表以下几米，在干旱期由于过度开采出现了时间上的显著变化。特别是在19921995年间，这种现象特别明显，测压水位的峰值下降了15m。因此，VMSR发生了大面积的地面沉降，造成了建筑物和大型公共设施的破坏（Mulas等，2003；Martnez等，2004）。Toms等人（2005，2006），Toms（2009）和Herrera等人（2008，2009）利用SAR差分干涉测量法

12、测量了穆尔西亚大城市地区在干旱期间的地面沉降，探测到的最大位移是12cm。从水文地质学角度来看，第二单元的更深处存在其它砂砾石层，但是对它们的了解很少。Mulas等人（2003）提出了一个刻画和模拟19921995年旱期发生的沉降的模型。根据这个模型，当水从深部含水层最顶部的砂砾石层中抽出的时候，就会发生沉降。垂向梯度的产生使得地下水从上部含水层（浅部单元）朝着砂砾石层向下流动，引起水位下降。因此，浅部含水层的孔隙水压力降低，含水层发生固结。超孔隙水压力的消散速率是地面沉降计算中的一个关键问题。对5个未扰动的淤泥和粘土样品进行了渗透性的实验室测试，结果表明渗透性非常低且具有变化性（8.2010

13、11 6.24108 m/s）。因此，可以推断固结过程是缓慢的。在这种认识下，为了测量弱透水层和含水层中的测压水位，穆尔西亚市安装了分层水压计，结果显示了这两个层的测压水位变化非常相似，极有可能是因为易于排水的砂层的存在。3 沉降数据本文的地面沉降测量是通过称为相关像元分析技术（CPT）的PSI技术获得的。可以在Mora等人（2003）和Blanco等人（2008）的文献中找到关于该技术的详细描述，本文在此只做简单描述。3.1 方法概述 用两幅SAR图像得到差分干涉测量相位（int）的公式为（Hanssen，2001）：int = flat + topo + mov + atmos + noi

14、se （1）式中，flat是不考虑地形条件下与测距差有关的扁平球体部分，topo是地形相位，mov是发生在两幅SAR图像获取期间的沿着卫星视线向（LOS）测量到的地面形变产生的相位，atmos是由于大气影响产生的相位，noise是残余噪声源。式（1）中的前两项可以通过分析获得。特别是，flat是很容易知道的，topo可以利用外部的DEM中计算得到。CPT算法（Mora等，2003）假设，与变形相关的相位组成（mov）可以分成两个新的相位项，一个是由于线性变形（linear）产生的，即整个时段内恒定速率的变形，另一个是由于非线性变形产生的（non-linear）。因此，CPT的应用分为两个主要的

15、步骤，即提取线性和非线性变形项。线性项的提取包括估计平均变形速度和DEM误差。这个估计是通过调整一个模型函数来实现的，该模型函数只针对那些显示出比较好的干涉时间相干性的像素而建。非线性项是应用空间-时间过滤器去除大气影响并识别非线性变形的低分辨率和高分辨率部分来提取的。与非线性变形部分比较的时候，大气隔离是可能发生的，在某种程度上，是由于大气影响具有时间和空间上的不同行为。关于CPT处理方法的其他说明如下：首先，通过对同一传感器获取的图像计算干涉图的不同子集（没有交叉干涉图），并利用奇异值分解法将变形结果联系起来，可参考Blanco等著作（2005, 2008），这样就有可能综合使用来自不同传

16、感器（ERS和ENVISAT）的图像。第二，缺少某段时间内的有效图像（如1994年和2001年）是可以解决的，一种方法是增加干涉图的最大容许时间基线（例如3年），目的是将数据空白之前和之后获取的图像联系起来，另一种方法是对数据空白扩大时间过滤器的时间窗。用这种方法的话，即使缺少某些时间间隔的测量数据，但整个时间序列还是完整的。这种方法在大部分为线性变形的假设下是有效的。3.2 数据集和详细处理过程本文分析的SAR数据总共由81幅图像构成，包括ERS-1（6幅）、ERS-2（56幅）和ENVISAT-ASAR（19幅），图像获取时间从1993年4月2007年3月。为了进行处理，提取了许多大约10

17、10 km的单元格，对应于穆尔西亚的城市地区。根据SAR图像对形成的所有可能的干涉图中，只选出了185幅干涉图进行CPT处理，借助于空间基线、时间基线和多普勒质心差形成的3D空间的Delauney三角剖分。选择被限定在那些空间基线和时间基线分别小于250m和1000天，多普勒质心差在800Hz以下的干涉图。用于抵偿干涉相位的地形部分的外部DEM的精度为25 m25 m，属于IGN（National Cartographical Service of Spain，西班牙国家地图服务）的地图数字数据库E20。利用多层相干标准已经选择出了具有可用的地面形变信息的像素，该标准就是针对大部分干涉图只选择

18、那些相干性超过某一阈值的像素（本例中，至少40%的干涉图集分别是0.6，0.5和0.4）。这种多层方法的优势是使低质量像素也能参与到沉降计算中，因为它们很好地利用了那些以前处理过的高质量像素（Blanco等，2006，2008）。相干性计算意味着干涉图的一种空间平均，也称之为“多视（multi-looking）”，降低了SAR图像原始分辨率。本次工作利用“方位角中的15个像素距离中的3个像素”的多视操作，因此最终的地面分辨率为60 m60 m。3.3 DInSAR处理结果和验证 利用CPT技术获得了穆尔西亚市19932007年间的总变形（LOS位移）图，该图的平均密度为每平方公里76个相干像素

19、。从图中可以发现，该市南部和东部的地面沉降量较大，整个研究期内大于5cm，而北部地区的情况较为稳定。在进一步应用上述沉降估算之前，已经将其同原位变形仪的测量数据进行了比较。从2001年2月至2007年3月，有15个监测深度大约为15m的变形仪提供了沉降测量数据。通过对SAR数据分析获得的形变信息与沿着LOS投影的变形仪网络测量数据的比较来评价CPT的准确度。为了使这种比较更容易，对LOS投影的变形仪时间序列在与SAR图像相同的时间间隔内进行插值，即20002007年。从与变形仪的比较中计算出两个质量参数：（1）CPT和LOS投影的变形仪变形时间序列之间的绝对差的平均值和标准差（4.54.1 m

20、m），（2）CPT和LOS投影的变形仪变形时间序列之间的差的平均值和标准差（2.64.7 mm）。这些值与其他一些研究人员在相似的试验中获得的数据是一致的（参见Strozzi等（2001），Hanssen（2003），Colesanti等（2003），Casu等（2006），Crosetto等（2008）， 及Herrera等（2009）。4 沉降计算本文利用一个一维模型来模拟测压水位变化产生的沉降。这个模型假设变形直接由来自于测压水位变化产生的垂向有效应力的变化引起，弱透水层孔隙水压力平衡与砂砾石含水层水位变化是瞬时完成的。值得注意的是，这种假设意味着当测压水位稳定的时候，固结尚未发生。正如

21、第2节论述的那样，这种假设在本研究中可以接受。穆尔西亚市面积有几平方公里，研究深度大概为15m。利用现有数据无法进行全三维的水流与沉降分析。这也是为什么已经选择了几个一维模型代替含水层弱透水层系统的全局模型进行分析的主要原因。接下来将介绍：（a）利用19931995年期间的DInSAR数据计算的弱透水层变形参数的模型识别过程；（b）模型对几个井的沉降历史预测情况；（c）沉降预测结果的验证，比较沉降模拟值与19952007年的DInSAR数据。4.1 模型识别地面沉降评价在岩土工程中属于常规工作。本文考虑三种参数来评价沉降：（i）潜在可压缩土层的厚度，（ii）应力状态的变化，及（iii）将前两个

22、参数关联起来的系数。我们认为，应力状态的变化是由水位变化引起的。因此，这个关系可以表达为： = h Ssk D = h Sk （2）式中，是沉降量（m），D是可压缩土层的厚度（m），h是水位下降量（m），Ssk是弱透水层的储水率（m1），Sk是弱透水层的储水系数（无量纲）。Sk或Ssk代表着弱透水层的可塑性，它们随着应力变化而变化。根据测压水位（H）以及测压水位是高于还是低于有记录的最大的水位降（Hp），而采用不同的值。因此，Hp相当于众所周知的岩土工程中所说的前期固结压力，代表着土体承受的最大有效应力。这个值将弹性的可恢复的变形与非弹性的、不可恢复的变形区分开： （3）值得注意的是，弹性储水

23、率的定义表明水是不可压缩的。对于松散的冲积含水层系统而言，这个观点是成立的，因为这样的系统其储水系数一般都高于因孔隙压力引起压缩形成的储水系数（Galloway和Hoffmann，2007）。利用24个现有观测井的水位数据计算了弱透水层的非弹性储水系数，这些井所在位置的DInSAR变形量也是已知的（图1和图2）。利用这些数据可以绘制出表示测压水位变化和弱透水层变形关系的应力应变曲线图（图1和图2）。利用Riley（1969）提出的图示法来计算非弹性储水系数（图1B和表1），以此校正模型。这个方法已经被使用多次（Hoffman，2003；Hoffmann 等，2003；Bubey，2003；To

24、ms等，2006；Galloway和 Hoffmann，2007；Zhang 等，2007a，b），用来确定应力应变曲线（弹性或者非弹性）上的一段的斜率： （4）式中，D是由于测压水位下降h之后引起的沉降。19931995年间的水位是VMSR历史上发生过的首次下降幅度最大的。因此，认为这个时期计算出的储水系数是非弹性的。这一假设被证明是恰当的，因为测压水位下降会显著减少孔隙水压力，从而使弱透水层经受历史上的最大有效应力。利用穆尔西亚市已有观测井资料计算出的非弹性储水系数（Skv）为4.900104 3.058103，平均值为1.5131030.680103。将弹性储水系数（Ske）作为非弹性储

25、水系数（Skv）的一个百分比来求取。为此，用回弹指数（Cs）与压缩指数（Cc）的比值Cs/Cc等同于Ske/Skv。按照UNE 103-405-94 （UNE，1994）标准进行了139次固结试验从而获得了这个比值，估计Cs/Cc为15%。利用穆尔西亚市观测井资料计算出的弹性储水系数的范围为7.3501054.587104，平均值为2.2691041.020104。图1 A：沿LOS投影的沉降量和测压水位，用于模型识别和预测B：应力应变曲线和非弹性储水系数的计算示意图图2 测压水位变化与沿LOS投影的DInSAR沉降量；几个井的应力应变曲线表1 非弹性储水系数（Skv）和弹性储水系数（Ske）

26、，前期固结压力水头（Hp），模型1995-2007年的误差，CP：相干像素4.2 模型运行每个井的Skv和Ske计算出来之后，将数据外推到其他测压水位时间序列上，以检验模型的模拟效果。因此，这种外推相当于对未来情境的一种预测，但又具有关键的数据优势能够去验证预测结果。鉴于对弱透水层压缩变形的水文地质概化，因此选择用一维方法进行研究。首先，确定前期固结水头（Hp），以此来区分测压水位变化引起的弹性变形和非弹性变形。将季节变化中的最低水位作为前期固结水头。本模型假设每个井中整个变形土体具有相同的Hp值。这个假设与Toms等人（2007）计算出的前期固结压力值是一致的，证实靠近地表处，潜水位上面的土

27、体是高度超固结的，潜水位之下为正常固结。将每个井已经计算出来的参数带入到式（2）中来估算沉降量。计算得到了19932007年的沉降演化（t）情况，它是每次水位变化（hi）引起的变形（i）的总和。要注意的是，式（2）假设弱透水层的水头与下部含水层的水头是瞬时平衡的。这种假设对于薄的弱透水层而言是一种合理的概化（Leighton和Phillinps，2003；Galloway和Hoffmann，2007）。图3显示了利用提出的公式获得的沉降预测结果。图3 模型预测结果虚线：测压水位；实线，DInSAR变形量；点线，模拟变形量。4.3 模型验证本文提出的模型已经被用于模拟19932007年之间发生的

28、地面沉降，这段时间的SAR图像是可用的。图3显示了几个井所在位置的模拟结果。这些结果已经被用来测试模型的自相容性。将DInSAR和模拟时间序列之间的平均误差和偏差作为两个序列之间拟合程度的质量指示标志。绝对误差相当于相同时间序列数据之间的平均绝对误差。表1显示了模拟值和实际DInSAR数据之间的比较结果。以本次研究中的24个井为例进行分析，19952007年的模型平均误差为0.83.6mm，极值为3.2和5.6mm（表1）。将这个误差作为所分析时间序列上每个日期的沉降预测值和DInSAR沉降测量值之间误差的平均值来计算。所有井的平均绝对误差是由模拟变形量和DInSAR测量变形量之间的绝对误差计

29、算出来的，为3.42.5 mm，极值（最大值和最小值）分别为5.7和1.4mm。这些结果与Herrera等人（2009）获得的结果相似，他们的数值模拟结果和DInSAR位移时间序列之间的绝对误差的平均值为5.54.7 mm。5 结 论本文综合使用了SAR差分干涉测量时间序列沉降值与19931997年的水位，确定了穆尔西亚市含水层系统的土体变形参数。获得的参数是非弹性储水系数（Skv）和弹性储水系数（Ske），已被用来进行简单的沉降预测。所提出的沉降预测模型分别将土层厚度和水位变化作为沉降的条件因素和触发因素，模拟了19932007年的沉降变化。通过与DInSAR19952007年的位移测量值进行比较，对沿LOS投影的结果进行了验证。本研究涉及了穆尔西亚市的24个井，它们的模拟值和DInSAR变形时间序列之间的绝对误差的平均值是3.22.5 mm。这些结果证明了PSI DInSA技术可以提供用于沉降计算的土体变形参数，从而成为研究含水层过度开采引起测压水位变化最终引发地面沉降的非常有用的工具。 译自Engineering Geology（2010）