地表沉降是一种缓慢变化的地质灾害，可形成一系列次生灾害链，与人类社会生产活动、经济建设、城市发展等重要因素密切相关.

我国的华北平原、长三角地区和汾渭盆地在过去几十年中是遭受严重地表沉降灾害的三个主要地区(He，2005).为跟踪地表沉降发展趋势，制定相应地表沉降防控措施，逐步开展了地下水动态监测、常规精密水准测量、基岩标和分层标监测.这些监测手段精度高，但空间分布十分有限，不能获得与大范围地表沉降空间分布一致的监测信息.GPS虽然能够获取大空间尺度地表运动信息(江在森等，2003；张跃刚等，2006；杨国华等，2007)，但由于多种因素的限制，离散的GPS监测点不可能密集布设.因此，对于大范围地表沉降监测，以上手段均难以获得空间连续分布的地表形变场信息.

近年来，研究人员基于干涉点目标分析思想陆续提出了多种新的DInSAR 地表垂直形变场探测方法，其中，永久散射体差分干涉技术(permanent scatterers，PSInSAR)应用较为广泛(Funning et al.，2007；屈春燕等，2011；罗三明等，2012b，2014).该技术在相干散射体上跟踪干涉相位随时间的演变过程，进而获取地表长期缓慢的变化信息.意大利学者Ferretti等通过实验证明，PS目标在经过8年的时间间隔后仍能保持稳定的散射特性(Ferretti et al.，2001，2002).

PSInSAR技术已在火山研究(Hooper et al.，2004，2007)、形变探测(Decriem et al.，2010；Chang et al.，2010；罗三明等，2012b)、地表沉降(郭炳跃等，2011；罗三明等，2014)等领域取得了大量研究成果.对已有研究结果的不完全分析显示，该方法一般基于一个片区(Frame)的地表垂直形变场空间分布 信息(如ASAR卫星的标准范围为100 km×100 km)，往往不能满足大空间尺度形变场研究.如华北平原、 长三角和汾渭盆地地表沉降面积分别达到48550 km²(何庆成等，2006)，10000 km²(郭坤一等，2003)和7000 km²(地球杂志，2012).

如何对PSInSAR技术进行进一步延伸与改进，以适应更大空间尺度地表形变场研究需要，中外学者已作了探索性工作.比如，Ge等利用PSInSAR方法分析了渤海湾大范围地表形变场(Ge et al.，2010)；Casu等反演了美国内华达州600 km×100 km 地表垂直形变场(Casu et al.，2008)；在2011年的Fringe研讨会上，Adam展示了利用PSInSAR技术获得的大时空尺度地表形变信息的研究成果(Adam，2011)；荷兰和意大利学者利用PSInSAR技术分别分析了荷兰全境和整个意大利半岛地表形变场变化规律(Cuenca et al.，2011).但这些研究成果目前仅限于会议报告，主要技术流程还不可获取.

本研究采用PSInSAR技术解析了华北平原3个轨道12个独立片区共272景ASAR重轨数据，通过对不同片区的PSInSAR技术结果进行基准转换与数据拼接，获取了该区域2006—2010年约12万km²地表垂直连续形变速度场(注：华北地区的地表形变主要在垂直方向，而ASAR的成像几何也主要反映的是垂向分量，因此本文的LOS形变主要反映的是垂向形变)，并利用同期地表水准测量结果 对本文提出的数据处理策略获得的结果进行了比较. 2 方法

采用PSInSAR方法进行多轨干涉数据的处 理，目前主要还是按照不同片区分别处理.PSInSAR 方法需要对同一片区重轨数据确定一期主影像，把其余影像作为辅影像与之配准.由于不同轨道(Track)主影像选取的参考基准不一致，不同片区PS目标反演的地表形变场因时间参考点的不同而在基准上必然存在差异；在相位解缠过程中，每个独立片区具有自己的参考PS目标，使得各片区反演的地表变化信息在空间上存在一定误差；不同片区的干涉数据因获取时间不同、地理位置不同，其大气效应影响程度也就不同；不同轨道的垂直基线不同而产生的轨道误差也不相同等等.这些主要因素使得不同片区地表形变场呈分块状态，相邻片区在空间上存在“错位”(图 1b)，各片区不能构成一个连续场.基准转换与数据拼接就是要消除以上多项误差对不同片区的基准产生的“错位”现象，将不同轨道表示的形变场转换到同一个参考基准上，形成一个连续的曲面(图 1c). 2.1 同名点的提取

所谓同名点是指同一地物在重叠区内不同影像对应的PS点.理论上，相邻片区重叠区内的PS点均应是同名点，位于相邻片区同名点间的距离应为0.但由于影像上单个象元反射类型及反射角度的不同，此时的PS点在不同影像里的位置不一定完全相同，甚至是不同的地物对应的PS点，它们之间的距离一般不为0.因此，本文采取迭代最近点法则提取同名点，并剔除离群值.

设相邻两个轨道的PS目标构成的点集分别为A和B，其重叠部分的点集分别为U和V(图 1a)，点r ₁(x₁，y₁，z₁)和点 r ₂(x₂，y₂，z₂)为位于点集U和V中的同名点，两点之间的距离为

图 1

Fig. 1

图 1 相邻两轨道PSInSAR结果在重叠区域产生的差异(a)相邻两个轨道PS目标的俯瞰图;(b)两轨道重叠部分侧视图;(c)配准后两轨道重叠部分侧视图，红色点集与黑色点集分 别为相邻两片区PS点，绿色点集为(b)中黑色PS点经过配准 后生成的新的PS点集. Fig. 1 Differences between adjacent two tracks of PSInSAR in the overlapped region(a)Aerial view of two adjacent tracks of PS target;(b)Lateral view to the overlapped region;(c)Lateral view of the overlapped region after registration.

根据公式(1)，通过迭代计算可以生成两个新的点集P={p_i}和X={x_i}，即点集A和B的公共部分.设点集P和X中PS的个数分别为N_p和N_x，这里有N_p=N_x=N，根据迭代最近点原则，点集P和X须满足条件

利用公式(1)生成的点集P和X中的点在满足公式(2)的情况下，但有可能大于单个象元的分辨率，这样的点不视为同名点，应予以剔除.因此，在识别同名点的过程中，在考虑点间距离最小的同时，还需要考虑SAR影像的分辨率. 2.2 基准转换参数的计算

(1)点距d_i权的确定

不同轨道影像的成像时间不同，地表同一地物的反射类型和反射角度可能发生变化，导致该地物在不同轨道影像上对应的PS点的位置可能存在差异.因此，根据2.1节介绍的方法获得的点距集d_i的大小可能不同，在计算转换参数时它们的贡献存在差异.为了提高计算转换参数的精度，需要引进d_i的权w_i，并定义w_i为

式中，d(p_i，x_i)由公式(2)求得，d_max为点距集d_i中的最大值，即d_max=max(d_i).

(2)计算点集P与X的协方差阵

设点集P和X的质心分别为μ_p和μ_x，则：

式中，.

根据μ_p和μ_x构建点集P与X的协方差阵为

式中，I ₃为3×3阶单位矩阵，Δ =[C₂₃，C₃₁，C₁₂]^T，C_ij=(Σ _px－ Σ _px^T)_ij.

(3)转换参数的计算

首先计算 Q(Σ _px)的特征向量和特征值，其中最大特征值对应的特征向量即为旋转参数 q _R，它的表达式为

而平移参数 q _T为

根据旋转参数 q _R和平移参数 q _T，构建点集P与X的转换参数 q ：

为了提高点集A与B基准拼接的精度，需要获得点集P与X的最佳转换参数.为此，根据下式计算最小均方差：

设点集A为参考点集，B^*为点集B经过配准后生成的新的点集，则有：

位于华北平原的研究区大范围处于沉降状态，北京、天津、廊坊、沧州等老沉降中心地带沉降速率处于较高发展态势，以城市为中心已扩展至整个华北平原.在降雨并不充沛的华北平原，大范围抽汲地下水是地表沉降的主要因素.本研究选取华北平原(115.32°E—118.79°E，36.81°N—40.58°N)作为多轨PSInSAR技术进行地表垂直形变场监测实验区(图 2)，采用StaMPS(Stanford method for persistent scatterers)方法(Hooper et al.，2007)解析PS目标.

图 2

Fig. 2

图 2 实验区PS片区空间分布图(T为轨道号，F为片区号)Fig. 2 Distribution of PS frames in the experimental area(T is track number. F is frame number)

实验数据是从欧空局申请的Envisat ASAR干涉数据，时间跨度为2006—2010年，涉及175、218和447三个轨道12个片区共272景数据(图 2).所有数据均为降轨数据，VV极化.3个轨道数据覆盖华北平原约12万km²，干涉参数如表 1所示.

表 1(Table 1)

表 1 ASAR雷达数据:垂直基线、时间长度、多普勒质心频率 Table 1 Interferometric parameters(perpendicular baseline，date length，Doppler centroid)of ASAR data 序号 No 轨迹号 Track 轨道号 Orbit 帧号 Frame 成像日期 During 主影像 Master 垂直基线 B⊥(m) 多普勒质心 频率(Hz) 1 2 3 4 175 20770~44818 2799 2006—2010 20090419 -30~511 259~281 20770~44818 2817 2006—2010 20090208 64~826 247~296 20770~44818 2835 2006—2010 20090524 -308~251 237~283 20770~44818 2853 2006—2010 20090524 -305~500 238~270 5 6 7 8 218 23318~44861 2799 2006—2010 20081029 -452~793 255~287 23318~44861 2817 2006—2010 20080716 -553~695 247~281 21815~44861 2835 2006—2010 20080611 -405~382 250~273 21815~44861 2853 2006—2010 20081029 -323~457 234~267 9 10 11 12 447 26052~45090 2799 2006—2010 20090227 -306~457 265~289 26052~45090 2817 2006—2010 20080801 -285~553 258~288 26052~45090 2835 2006—2010 20080801 -286~544 245~273 26052~45090 2853 2006—2010 20080801 -326~534 245~267

表 1 ASAR雷达数据:垂直基线、时间长度、多普勒质心频率 Table 1 Interferometric parameters(perpendicular baseline，date length，Doppler centroid)of ASAR data

表 1显示，12个独立片区重轨数据垂直基线大部分在500 m左右，多普勒质心频率在250 Hz上下.根据PSInSAR方法分析结果(图 3)，从北京—廊坊—天津到任丘、沧州经泊头、衡水、德州到南宫，12个片区重轨数据在5年时间内保持了良好的干涉性，较为完整地揭示了研究区地表垂直形变场空间分布特征.

图 3

Fig. 3

图 3 华北平原PS目标速率空间分布图Fig. 3 Spatial distribution of PS target rates in the North China Plain

华北地区一等精密水准测量的周期为3~5年 甚至更长，因而我们无法通过同期水准观测成果对 整个研究区的PSInSAR方法获得的地表形变场结果进行验证.从上世纪80年代开始，天津地区于每年10—12月进行一次地表沉降一等水准测量，具有良好的周期性和连续性，与本文采用的SAR数据具有较好的同步性.因此，我们对本文提供的数据处理策略反演的地表形变场的验证采用的是天津地区2006—2010年的同期一等水准测量结果. 3.2 处理结果

利用StaMPS方法(Hooper et al.，2007；罗三明等，2012b)，分别对12个独立片区完成PS目标的识别，根据本文提供的数据处理策略，对12个片区的PS目标进行基准连接，构成由12个片区组成的连续地表垂直形变速度场(图 3).在研究区约12万km²范围内，共识别出1202958个PS目标，每km²包含约10个PS，该密度足以全面揭示整个研究区地表垂直形变场空间分布特征及其与该区内断裂的关系(图 4).

图 4

Fig. 4

图 4 华北平原PSInSAR方法反演的地表垂直形变场速率与断裂空间分布(红色实心圆为5级以上地震)Fig. 4 Rates of vertical ground deformation derived from PSInSAR inversion，active faults and epicenters of earthquakes M5 or greater(red solid circles)in the North China Plain

图 5为天津地区的PS目标与水准测量点在空间分布上的套合图.天津地区一等水准观测里程每年约2000 km，观测数据的平差计算以位于蓟县的基岩点为参考基准，采用分期平差(罗三明等，2012a)数学模型对2006—2010年每期观测值进行观测误差的消除或最大限度的弱化处理，共获得290个具有同期观测成果的水准点.

图 5

Fig. 5

图 5 天津地区PSInSAR与同期水准测量结果套合图(蓝色实心圆为水准点空间分布，AB与CD为两个剖面的位置)Fig. 5 Nested map of PSInSAR and coeval leveling measurements(Blue solid circles denote benchmarks. AB and CD are locations of two profiles)

由于PS目标的平均密度远大于水准点的密度，为了对二者的结果进行比较，我们采取了如下处理方法：首先将水准测量平差结果投影到雷达视线向(Envisat ASAR的视线向入射角为23°)，再假设每个水准点周围100 m范围内地表形变梯度是一致的，该范围内的PS目标与地表水准测量结果具有相同的形变梯度.根据该原则，在290个具有同期观测成果的水准点中筛选出了56个满足该条件的水准点，二者的比较见表 2.

表 2(Table 2)

表 2 水准观测结果(LOS)与PSInSAR方法反演结果对比(2006—2010)(单位：mm) Table 2 Comparison of PSInSAR and leveling measurements in 2006 to 2010(unit:mm) 序号 水准结果 PS结果 误差 序号 水准结果 PS结果 误差 序号 水准结果 PS结果 误差 1 -3.3 -2.3 -1.0 20 -49.5 -43.2 -6.3 39 -30.3 -39.2 8.9 2 -10.1 -10.9 0.8 21 -55.1 -46.2 -8.9 40 -29.1 -37.3 8.2 3 -9.4 -11.9 2.5 22 -50.6 -43.1 -7.5 41 -80.1 -76.0 -4.1 4 -10.8 -14.2 3.4 23 -33.6 -33.2 -0.4 42 -20.7 -16.7 -4.0 5 -6.2 -1.7 -4.5 24 -40.5 -44.8 4.4 43 -16.0 -18.8 2.8 6 -8.0 -10.0 2.0 25 -15.2 -16.5 1.2 44 -60.2 -65.2 5.0 7 -10.1 -10.7 0.6 26 -16.2 -17.4 1.2 45 -66.0 -71.9 5.9 8 -33.8 -34.9 1.1 27 -20.1 -22.1 2.1 46 -39.3 -36.5 -2.8 9 -24.9 -21.7 -3.3 28 -20.5 -21.4 0.8 47 -44.1 -46.3 2.2 10 -53.4 -55.3 1.9 29 -17.4 -17.6 0.2 48 -30.8 -36.5 5.7 11 -56.7 -58.0 1.3 30 -12.5 -16.6 4.2 49 -41.9 -33.7 -8.2 12 -45.7 -47.6 1.9 31 -12.7 -18.5 5.8 50 -42.4 -33.7 -8.7 13 -12.1 -13.6 1.5 32 -20.6 -24.6 3.9 51 -18.9 -25.0 6.1 14 -16.5 -14.7 -1.7 33 -30.7 -29.7 -1.0 52 -19.2 -24.3 5.1 15 -26.8 -26.4 -0.4 34 -33.4 -37.7 4.3 53 -35.3 -41.2 5.9 16 -30.0 -32.8 2.8 35 -23.5 -28.0 4.5 54 -11.9 -20.3 8.3 17 -62.8 -59.4 -3.4 36 -22.4 -26.7 4.3 55 -30.5 -38.0 7.5 18 -64.3 -60.5 -3.8 37 -10.9 -14.3 3.4 56 -26.8 -23.9 -2.9 19 -51.3 -43.2 -8.1 38 -38.6 -47.2 8.6

表 2 水准观测结果(LOS)与PSInSAR方法反演结果对比(2006—2010)(单位：mm) Table 2 Comparison of PSInSAR and leveling measurements in 2006 to 2010(unit:mm)

为进一步分析PSInSAR方法反演结果的可靠 性，我们在图 5中的AB和CD位置绘制了两个剖 面，图 6为PSInSAR和水准二者反演的地表垂直形变速度场在两个位置上的剖面图.图 6显示，水准测量结果和PSInSAR方法反演结果在趋势上具有良好的一致性.

图 6

Fig. 6

图 6 PSInSAR结果与水准结果剖面AB和CD的比较Fig. 6 Comparison of PSInSAR and leveling measurements on profiles AB and CD

已有研究成果显示，在植被稀疏、岩石裸露的区域(如我国西部)，PSInSAR结果的精度可以达到2 mm；在城市地区，虽然因人工建筑物较多而可以识别出较高密度的PS目标，但相应地干扰源也较多，因而降低了PS目标反演的地表形变信息的精度，一般情况下可达到5 mm；在植被较厚、人工目标较少的地区，因可识别的、能满足良好干涉条件的相干体既少而且离散，PSInSAR的精度一般在1 cm左右.

为了对12个片区拼接后的可靠性进行评估，我们计算了各重叠区同名点的标准差(表 3).由表 3可以看出，对各片区进行拼接后，17个重叠区同名点的标准差在0.83~3.34 mm之间，其中超过3 mm的有3个，说明各片区识别出的PS目标具有较高的质量，由此获得的大空间尺度地表形变场信息是可靠的.其中，精度最高的为研究区东部的175轨道，4个片区的拼接精度均在1 mm以内.根据表 3还可以看出，无论是同轨相邻片区(Frame)间的拼接(南北向拼接，表中前3列)，还是不同轨道相邻片区间的拼接(东西向拼接，表中后2列)，研究区东部的拼接精度都要高于西部，这可能与研究区东、西部地表植被覆盖程度存在差异有关.

表 3(Table 3)

表 3 重叠区同名点标准差(单位：mm) Table 3 St and ard deviations of corresponding PS in overlapped regions(unit: mm)

表 3 重叠区同名点标准差(单位：mm) Table 3 St and ard deviations of corresponding PS in overlapped regions(unit: mm)

根据表 2统计，在56个水准点中有37个与PSInSAR方法反演结果的误差小于5 mm，占总点数的66％，有2个水准点结果与PS方法结果相差8.9 mm.进一步计算可得中误差为4.72 mm，满足我国对城市地面沉降观测5 mm/km的偶然中误差和10 mm的全中误差的精度指标要求(北京市测绘设计研究院，1999).

表 2中，少数水准点与PSInSAR方法反演结果的差值较大，这可能与二者在空间上的分布位置和密度有关.一般情况下，在SAR图像上无法解析出水准测量点，即PS目标和水准测量采集的数据不在同一点上，尤其在地表垂直形变梯度较大的地区，位置偏差直接导致形变信息的不同，二者的差值就可能产生较大差异.由剖面图 6可以看出，由于水准点的密度要远小于PS目标的密度，速度场剖面在经过沉降中心地带时二者的差值明显偏大.这也说明：对于不同的对地观测手段获得的观测结果，一般 宜进行趋势间的相互比较，如果将其中一种结果作 为标准去衡量其它观测手段获得的结果可能有失严谨. 3.5 结果分析

华北平原分布广泛的第四纪松散沉积物构成了该区地下水含水系统，经历了第四纪以来的新构造运动，形成一个水文地质结构复杂而又相互联系的地下水系(何庆成等，2006).本研究结果显示(图 4)，沉降区主要沿两个方向分布，一个沿廊坊—天津市区—滨海新区呈北西向分布，另一个沿天津市区—静海—沧州—泊头—德州呈北北东向分布，而天津位于两个沉降带的交汇部位.天津地区地表主沉降区集中在工业用水相对集中的津南区—滨海区一带，两地中心沉降区平均沉降速率分别达到-61.6 mm/a和-66.7 mm/a，这可能与近10多年来天津地区工业重心向津南和滨海转移导致地下水需求扩大有关.而位于沧州—泊头—德州一带，年平均沉降速率分别达到-34.7 mm/a、-37.7 mm/a、-29.3 mm/a，地表沉降速率仍处于较高发展态势. 北京东部沉降区2006—2010年沉降量为173.5 mm，年平均沉降速率达到-34.7 mm/a，表明北京东部地区地表沉降趋于加速发展趋势.

从图 4可已看出，两个地表沉降带恰好与北向的廊坊—武清断裂，以及北北东向的沧东断裂位置吻合.特别是在天津城区、沧州、泊头—德州等地由多个沉降漏斗区构成了一条北北东向的主沉降带.华北地区自有记载以来构造活动一直很强烈，是中国大陆东部主要地震活动区，从历史地震来看，该区发生过多次5级以上地震，除唐山老震区外，其余地震大多发生在断裂附近，说明研究区地震活动与断裂的分布确实存在一定联系.地下断层及周边发育着众多裂隙，而地下水沿裂隙存储、分布和迁移，地下水位的变化可能通过改变断层变形行为方式而影响着区域构造变形图像.地下水长时间抽取，造成局部地下水位下降，地表下沉，而周边地下水会沿着裂隙快速迁移补充，这可能在一定程度上反映了地表沉降带与断裂带吻合的现象，同时也反映出由地下水开采导致的地表沉降受断层分布的控制. 4 结论

本研究选取的ASAR干涉数据，其重轨周期最短为35天，采用PSInSAR和迭代最近点算法处理，可以获得大面积地表垂直形变速度场，从信息的时间频度和空间分布上大大弥补了常规对地观测时空信息不足的缺点.通过在华北平原应用，并得到以下认识和结论：

(1)提出了一种基于基准连接的多轨PSInSAR获取大空间尺度地表形变场缓慢变化过程的有效方法，扩展了PSInSAR技术在空间尺度应用方面的局限性，使大面积获取PS点成为可能.

(2)对17个重叠区同名点误差统计显示，标准差在0.83~3.34 mm之间，说明各片区识别出的PS目标具有较高的质量，获得的大空间尺度地表形变场具有较高的可靠性.进一步分析表明，研究区东部的拼接精度都要高于西部，说明东部的PS目标的质量要高于西部，这可能与东西部地表植被存在差异有关.

(3)采用ASAR干涉数据，通过本研究提出的处理策略，获取了华北平原约12万平方公里2006—2010年间的地表垂直形变速度场及其空间分布特征.结果表明，以北京、廊坊、天津、沧州、泊 头—德州等城市为中心，地表沉降已扩展至整个华 北平原，沉降中心区5年内的地表沉降量分别达到174、132、321、173 mm和188 mm，表明实验区地表沉降处于较高发展态势.

(4)所获得地表垂直形变速度场与同期一等水准测量结果进行了验证，精度达到4.72 mm，为大时空尺度地表垂直形变场和跨地区、跨区域大型交通枢纽工程变形等监测提供了一种行之有效的观测新途径.

(5)两个地表沉降带分布与北向的廊坊—武清断裂和北北东向的沧东断裂位置吻合.这可能与地下断层及周边发育着众多裂隙有关.地下水沿裂隙存储和迁移可能通过改变断层变形行为方式而影响着区域构造变形图像，反映出由地下水开采导致的地表沉降受断层分布的控制.