2012, Vol. 55
2. 台湾大学地质科学系, 台湾 台北 106;
3. 香港理工大学土地测量与地理资讯学系, 香港 九龙
2. Department of Geoscience, National Taiwan University, Taiwan Taipei 106, China;
3. Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Hong Kong, China
监测许多地球物理现象(如地震运动、火山喷发、地面沉降等)所引发的地表形变三维位移，对于准确认识地球物理现象的发生机制、运动过程和发展趋势，具有重要的地学数据支撑作用[1-4].目前，对于大范围区域地表形变的监测，常采用全球定位系统(Global Positioning System, GPS)或卫星雷达干涉测量(Interferometric Synthetic ApertureRadar, InSAR)方法.
GPS测量的优势在于对水平方向的地表位移具有较高的测量精度和较高的敏感度[5-7]，采用较长的观测时段，GPS 能探测出地表mm 级的水平位移，在垂直方向的形变测量精度稍低，但对于大尺度的垂直形变仍然具有较高的测量敏感度.其不足之处是基于地面稀疏点的离散测量，受野外地形条件、设备安装等因素限制，GPS 测站分布通常较为稀疏，尤其对于大范围的地表形变监测来说，单一采用GPS方法测量地表形变，往往难以获得空间连续的形变分布场.
InSAR 观测的优势在于其基于影像面的大范围测量，具有空间分辨率高、覆盖范围广、空间无接触遥感等技术特点，尤其对于大尺度的地表形变，如火山运动、地震位移等，InSAR 能够获得空间连续的地表形变场.对于随时间缓慢累积的小尺度地表形变，如城市地面沉降、震前缓慢位移等，基于永久散射体(Persistent Scatterers, PS)的InSAR 能够弥补干涉相位受时空失相关影响的不足，从稳定的PS单元上获得高密度的地表形变场[8-14].
将高精度的GPS 水平位移测量与高密度的PS-InSAR 垂直形变观测技术相结合，有望极大提升大地测量技术监测地表形变的可靠性与完整性，两者联合形成技术优势互补[7]，为地球物理研究和形变灾害监测提供完整的三维地表位移场，这对于探查区域地表形变的空间分布规律与发展趋势，并进一步认识相关地球物理现象的形成机理与诱发因素，具有重要的数据支撑作用.
由于InSAR 是沿雷达波视线方向(LOS)的一维位移测量，而对许多地球物理现象(如地震形变、火山喷发、城市地面沉降等)的解译或反演，需要获得地表在水平方向和垂直方向的三维位移，因此需要将雷达在LOS方向的形变转换为在垂直方向的位移.如果地表仅存在垂直方向位移而没有水平位移时，依据雷达波侧视角易于将LOS形变转化为垂直方向形变，但对于地表发生三维位移形变来说，则必须借助其它观测数据如GPS 观测值，将InSAR在LOS方向的形变转换为垂直方向的位移.
为解决GPS 与PS-InSAR 技术联合中面临的形变投影转换问题，并联合GPS与PS-InSAR 技术监测地表三维位移，本文将地面GPS 观测数据与PS-InSAR 相结合，通过建立雷达LOS 形变量与GPS观测量的几何投影转换模型，构建GPS台站与雷达PS目标(类似于天然GPS 测点)监测网，采用区域网参数平差算法，提取研究地区高精度高密度的地表形变场.以地质构造极为活跃的台湾岛及其西南屏东和高雄地区为例，采用该地区1995-1999年间获得的ERS-1/2卫星雷达影像和地面GPS观测数据，通过形变投影转换和联合建网数据处理，提取该地区由于板块构造运动和地下水抽取而导致的显著三维地表形变场.
2 GPS与InSAR的形变投影转换模型采用GPS接收站测量得到的地表形变量通常是沿东(E)、北(N)、垂直(H)三个方向的位移，而InSAR 测量所获得的位移是沿雷达视线方向(LOS)的一维形变[5， 15]，这两种方法存在形变测量方向的差异.对于许多地球物理现象(如地震形变、火山位移等)的解译和判读，需要获得地表在水平方向(包括东向、北向)和垂直方向的位移量，因此，对于InSAR 在雷达视线方向的地表形变，需借助相关数学模型把LOS方向形变量投影转换为易于理解的水平或垂直方向的形变量[15].
设地表的实际形变矢量为V，在沿东、北、垂直三个方向的形变分量为(ΔE，ΔN，θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH)，如图 1 所示.采用GPS观测得到其在东、北方向的水平形变分量为(ΔEGPS，ΔNGPS)，令VGPS=(ΔEGPS，ΔNGPS).假设SAR 卫星运行轨迹方向(SAR 影像方位向)与北方向的夹角为Φ(见图 1)，借助Φ 角将GPS水平位移量(ΔEGPS，ΔNGPS)转换为其在SAR 影像的距离向(Range)和方位向(Azimuth)的形变投影分量(ΔR，ΔA)，令VSAR=(ΔR，ΔA)，两者之间的转换关系为:
|
(1) |
|
图 1 雷达LOS与GPS形变投影几何转换 Fig. 1 Geometric transformation between InSAR LOS and GPS displacement |
如图 1所示，方位向相对于距离向和LOS方向均为垂直关系[15]，因此方位向的形变量ΔA在LOS方向形变投影为零，而距离向(Range)的形变量ΔR与LOS方向的夹角为θ，因此ΔR在LOS方向的形变投影分量为(ΔEGPScosΦ +ΔNGPSsinΦ)·sinθ，垂直方向的形变量θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH在LOS 方向的形变分量为θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH·cosθ，两者之和即为地表三维形变量投影到雷达视线方向(LOS)的形变Δr，因此有:
|
(2) |
上式中，Δr是InSAR 测量得到的LOS方向的形变量，θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH是地表在垂直方向的位移量，θ 为SAR 卫星的成像侧视角.
把(2)式改写为如下形式:
|
(3) |
上述的第(3)式建立了雷达LOS形变量与地表垂直方向形变量的转换关系.从该式可以看出，如果已测得地表的水平位移量(ΔEGPS，ΔNGPS)以及基于卫星轨道参数计算出的夹角Φ、θ，使用(3)就可以把InSAR 在雷达视线方向(LOS)形变量Δr转换为地表在垂直方向上的形变量θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH.
为定量评估水平位移观测误差可能对LOS 形变转换为垂直形变的影响，对(3)式中的ΔEGPS和ΔNGPS分别进行偏微分，获得如下的误差传播定量关系:δΔH=
|
(4) |
式中，mΔH为垂直位移中误差，mΔEGPS和mΔNGPS分别为东西向和南北向的位移中误差.据(4)式表明，水平位移误差对垂直形变转换的影响系数是角度(卫星轨迹方位角Φ 和SAR 侧视角θ)的三角函数，在θ小于45°时(如ERS卫星影像中心约为23°)，系数值小于1，水平位移误差对LOS的垂直形变转换误差是缩小作用(而非放大效应).若水平位移值在东西向和南北向分别有2 mm 的误差(总的水平位移误差为2.8 mm)，引起InSAR LOS 形变转换为垂直位移的误差量为0.8 mm.需要提及的是，随着卫星侧视角θ 的增大(如PALSAR 侧视角为34°)，水平位移误差对LOS位移向垂直形变转换的误差也会随之增大，超过45°则可能呈现误差放大效应.
3 卫星轨迹方位角和SAR 侧视角的计算模型如前述(3)式，要实现InSAR-LOS形变量向地表垂直方向形变量的变换，需要计算卫星轨迹的方位角Φ(即SAR 影像的方位向与北方向的夹角)、以及影像像素所对应的SAR 侧视角θ.
根据图 1所示几何关系，采用SAR 成像期间的卫星轨道状态矢量(即位置矢量和速度矢量)计算卫星轨迹方位角.对于每一景SAR 影像，均附带有获取该幅SAR 影像的卫星轨道数据，如ERS-1/2 卫星附带有5 个历元时刻的轨道点三维坐标(X，Y，Z)，也可以从国际上有关的轨道数据分发中心(如荷兰DEOS、德国D-PAF)获取SAR 影像成像期间的精密轨道矢量[16].
以ERS-1/2卫星轨道为例，假设Si、Sj为一幅SAR 影像成像期间卫星轨道的2个历元时刻，将空间轨道点的三维空间直角坐标(X，Y，Z)转换为大地经纬度坐标系(B，L，H)，其中B为大地维度，L为大地经度，H为大地高，再采用高斯-克吕格投影方法，将大地经纬度坐标(B，L，H)数学投影为高斯平面坐标(x，y)，采用如下公式计算卫星轨迹的方位角Φ:
|
(5) |
通过上述的两个投影变换过程，得到卫星轨迹在平面上的投影，如果投影点均位于一条直线上，则整幅SAR 影像可以采用同一个Φ 值进行计算.为考察卫星轨迹方位角在每幅影像成像期间内的变化情况，选取后续实验所用到的其中5 幅ERS-1/2SAR 影像的轨道资料进行测试计算，基于不同时间跨度(1995-1999年)的5幅SAR 影像及其DEOS精密轨道数据[16]，在其成像时间内分别取5对历元轨道位置的矢量组合，计算卫星轨迹与北方向之间的夹角Φ，计算结果如表 1所示.
|
表 1 ERS-1/2卫星轨迹与北方向方位角计算表 Table 1 The azimuth angle between trace of ERS-1/2 satellites and the north |
从表 1可以看出，在不同时间跨度的5幅影像成像期间内，卫星轨迹线与北方向的夹角较为稳定，均值都为北偏东12.5°左右.顾及方程式(3)中的(ΔEcosΦ+ΔNsinΦ)·tanθ，计算值均为角度Φ的正弦或余弦关系.因此，通过对理论公式的分析和实际轨道数据的测试计算，表明对整幅SAR 影像可以采用同一个Φ 值(如ERS-1/2 卫星取Φ=N12.5°E)，以此简化GPS 与InSAR 形变投影转换的计算量，且不损失形变计算精度.
除了计算卫星轨迹方位角之外，方程式(3)中还需要计算SAR 影像像素对应的侧视角θ.如图 1所示，GPS 测站的大地坐标和SAR 影像经地理编码后像素点大地坐标为已知，卫星的轨道位置可基于5个时刻(或多个时刻)的坐标采用如下多项式进行轨道拟合.
|
(6) |
上式中，(Xt，Yt，Zt)表示卫星的位置矢量，t为飞行时刻，Ai，Bi，Ci为拟合系数.根据上述的轨道矢量拟合方程，可计算出SAR 成像期间任意时刻的轨道位置.对于SAR 影像中任意一行像元的方位向时刻ti，利用SAR 影像起始成像时间t0、雷达脉冲重复频率(PRF)，计算其成像时刻:
|
(7) |
将SAR 影像像元的方位向成像时刻ti带入轨道方程(6)，则可以计算像元所对应的卫星轨道坐标，然后基于图 1所示的θ 所在三角形几何关系，即可计算θ 值.
4 PS-InSAR与GPS联合建网数据处理与基于单个干涉像对的D-InSAR 处理方法不同，PS-InSAR 是使用覆盖同一地区的多幅SAR 影像开展时序相位差分干涉，选取相位质量稳定的像素子集(PS)作为主要处理对象[8， 9， 17]，通过主从影像之间相位求差，即干涉技术，并顾及大气状态在近距离范围内的高度相关特性，通过邻域PS 相位再次差分，逐步分离和提取地表形变、地形误差和与大气贡献等信号分量.
利用前述的InSAR 与GPS的形变投影转换模型，可以实现PS-InSAR 在LOS方向的形变与垂直方向形变的投影转换.雷达PS 与GPS 联合建网用于提取地表形变的主要计算过程如图 2 所示，主要包括:InSAR 时序差分干涉、空间基线解算与网络约束平差，具体如下.
|
图 2 PS与GPS差分网络干涉数据处理流程 Fig. 2 The data processing flowchart of PS-InSAR and GPS networking technique |
假设有覆盖同一地区的N+1 幅时序SAR 影像，以其中一幅作为公共主影像[18]，其余N幅影像均配准并重采样到与主影像一致的像素空间，总共形成N个干涉对.将所有干涉对分别进行一次相位差分处理，形成N幅初始干涉图.干涉图中每一像素的初始干涉相位是多个贡献分量的总和，包括参考面相位、地形相位、形变相位和随机噪声等.
为从干涉相位中提取地表形变信号，需要分离和提取各信号分量.借助卫星精密轨道状态矢量并依据空间解析几何关系，可去除参考面相位[19]；借助精密轨道状态矢量和已有的数字高程模型(DEM)，可部分去除地形相位(由于DEM 可能存在高程误差)，得到二次差分干涉图[19]，余下的相位主要包括4个贡献分量，即
|
(8) |
上式中，第一项为残留高程相位，第二项为地表形变引起的相位.其中，λ表示雷达波长，Rm为雷达至地面目标的斜距，θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δH为DEM 高程修正值，Δr为地表沿LOS方向的形变量，ΦAtmk 为大气相位，nk为噪声相位和可能的非线性形变.
将前述的GPS水平位移与雷达LOS形变投影转换模型替换(8)式中的Δr，并令Tk为干涉对的时间基线，v为地面线性沉降/抬升速度，则(8)式可改写为:
|
(9) |
上述去除DEM 相位分量的差分干涉相位中仍包含有大气相位延迟的负面影响，顾及同一影像内大气状态的自相关特性[8， 9， 17]，从空间尺度上对邻域PS干涉相位再次差分，则可在一定程度上减弱大气相位延迟的影响.
对满足近距离条件的邻近PS目标建立相位差分方程，建立邻域PS点对的空间连接基线，其差分参数为邻近PS 的沉降速度差和高程修正差.采用复相关系数最大化方法求解空间基线差分参数[8， 9， 14]，满足相关性最好的一组解即为模型的匹配解.
在求解出邻域PS 的差分参数后，每条基线连接了基线两端PS 目标的形变速度差Δv与高程修正差ΔδH，这类似于大地测量中的水准观测网，大量的多余观测值(即PS连接差分参数)存在观测误差，出现闭合差ΣΔv≠0的情况.为探寻PS网络基线参数的全局合理性，以GPS 观测值作为约束条件，采用区域网参数平差法消除基线观测的几何不符值，并计算各PS的线性沉降速度量.
设vi、vj分别为PSi与PSj的线性沉降速度值，Δv为解算出的PS线性沉降速度差，以此作为平差处理的观测值.对PS 沉降速度网采用参数平差法处理，列出其误差方程:
|
(10) |
式中，ωv为模型的形变速度差与基线差分参数的残差，${{{\hat{\upsilon }}}_{j}}$、${{{\hat{\upsilon }}}_{i}}$为PS的形变速度估算值.
对于PS网中的每一条基线均建立上述误差方程，由所有基线的误差方程组成方程组，对误差方程组进行法化求解[14， 17]，即可解算获得每个PS 点在所研究时间跨度内的线性沉降速度v，由此获得地表随时间在垂直方向上的沉降/抬升形变时间序列以及在空间尺度上的形变分布.
5 实验数据与结果分析位于太平洋西北的台湾岛是世界上地质构造活动最为活跃的地带之一[2， 5， 15]，台湾地处菲律宾海板块(PhilippineSea Plate, PSP)与欧亚大陆板块(EurasianPlate, EUP)的碰撞活动带[15]，系由欧亚大陆板块、菲律宾海板块挤压而隆起的岛屿.以相对稳定的欧亚大陆板块作为参考，PSP每年以82mm/a的速度向西北方向挤压靠近(见图 3)，使得台湾岛上极易发生剧烈的地壳变形[5， 15].台湾西南沿海的屏东、高雄地区除受地壳构造运动影响外，其地下水的大量抽取也导致地表发生显著的地面沉降，根据屏东县沿海地区的航空影像调查显示，超过10347公顷的地区受地表沉降严重影响[5]，并导致重大的经济损失.
|
图 3 台湾地质构造框架图灰色方框为研究区域，浅灰色箭头表示PSP 板块向EUP 板块的挤压速度和方向，深灰色箭头表示构造逃脱，带三角形的线条为主要断层走向. Fig. 3 Taiwan geodynamic framework and structural mapThe gray oblique rectangle indicates the :^tudy area.Large light gray arrow indicates the direction and rate of the plate convergence of PSP and EUP.Small gray arrows indicate the tectonic escape.Major thrust faults with triangles are on the upthrust side. |
为定量监测台湾西南沿海屏东、高雄地区的地壳形变与陆地沉降，相关部门已于1995年建立了覆盖屏东高雄地区的“屏东GPS 网络"，包括由48 个监测站和1 个稳定基准站(S01R，位于稳定大陆框架的澎湖岛上，图 3中的三角形)组成的监测网[5]，监测站之间的距离为3~21km, 其分布如图 3 所示.从1995年8月至1999年8月，每年对这些GPS站进行精密测量，采用瑞士Bernese软件和IGS 精密星历处理GPS观测数据，以位于澎湖列岛上的连续运行跟踪站(白沙站，点名为S01R)作为稳定基准站，解算其余观测点相对于该点的形变，得到了台湾西南沿海地区的GPS地表形变水平位移量，如图 4所示，图中的矢量和误差圆标识了GPS观测位移的大小、方向及其标准差.
|
图 4 GPS测站的分布及其水平形变位移矢量 平面矢量代表水平位移，矩形方框为SAR 影像覆盖范围. Fig. 4 Distribution of GPS stations with horizontal movement velocity Horizontal displacement indicated by the vectors.The ground coverage of SAR acquisitions used in this study indicated by the oblique rectangle. |
图 4中的GPS 位移矢量展现了该地区水平位移形变场分布，该地区中央和西部沿海发生了显著的向西位移，越往西部测站位移速度有减少的趋势，年均水平位移量为30~50mm/a, 然而到了高雄沿海地区测站位移方向发生了明显的逆时针西偏南的位移，位移速度有明显增加的趋势.根据已有相关研究资料表明[5， 20-21]，台湾岛受菲律宾海板块和欧亚大陆板块的长期碰撞挤压，陆地表面发生剧烈的变形缩短，尤其是西南部的屏东高雄等地区存在构造逃脱现象，导致该地区地表发生显著的西向水平位移.
对于该地区的垂直形变场(如地表沉降或抬升)，由于受GPS测站数量与空间分布的限制，仅从现有的GPS观测数据难以获得该地区详细的垂直形变场空间分布及其发展趋势，需要联合长时间序列的卫星SAR 影像PS-InSAR 技术予以补充.
PS-InSAR 数据处理采用了ERS-1/2 卫星从1995-1999年间覆盖该地区的19幅单视复数影像(SLC)为数据源，选取其中的1998 年12 月3 日(Frame3145)影像作为主影像，其余18幅影像均被配准并重采样到主影像格网空间，18个干涉对的时间基线与空间基线分布如表 2 所示，最长的时间间隔约4年，最大的轨道空间距离超过700m.
|
|
表 2 研究地区的SAR影像干涉对 Table 2 Interferometric pairs used in the study |
对上述所有干涉对进行相位干涉预处理，采用Doris软件计算获得初始干涉相位图[19]，然后借助DEOS精密轨道状态矢量去除平地相位分量[16， 19]，采用JPL/NASA 的SRTM 数字高程模型(DEM)并依据干涉几何关系移除地形相位分量，得到各干涉对的差分干涉相位图.顾及InSAR 超短空间基线对形变提取的高敏感度因素，为考察短基线对地表垂直方向沉降测量的敏感能力，从19幅影像的组合中选取了垂直基线较小的干涉对进行差分干涉处理，其中干涉对31/01/1996-16/05/1996 形成的空间基线最小(B⊥ =19.2m, T=106days)，对应的差分干涉相位图如图 5所示.
|
图 5 差分干涉相位图(a)和解缠后的LOS方向形变图(b) Fig. 5 Wrapped differential interferogram (a) and unwrapped deformation map (b) |
图 5所采用的干涉对空间基线较短，垂直基线分量仅为19.2 m, 模糊度高约为450 m, 即使差分干涉所用的SRTMDEM 有误差，20m的DEM 误差引起的相位数值仅为0.28rad.因此，对于图 5a中的差分干涉相位图与图 5b中解缠得到的LOS方向形变图，其数值主要表现为地表形变引起的相位分量.图 5中的形变图较为清晰地表明，该地区屏东平原南部呈现明显的地表沉降，越往平原北部沉降量越小，直至出现微量的地表抬升，而在平原西侧的高雄地区则呈现明显的地面抬升.
上述干涉图是在有效空间基线极短的条件下得到的，对于其余更多的干涉像对，由于空间基线较长或时间失相关严重，差分干涉图中均未能呈现出较为清晰的干涉条纹，从这些干涉图很难获得地表在其余时期的形变资料，需要进一步借助基于稳定散射体(PS)的干涉测量方法探查地表的垂直形变.
采用相关系数阈值方法从时序相关系数图中选取均值γ≥0.3 的像素单元作为PS 目标，在SAR影像的公共覆盖范围内，共选取了3433 个PS 单元，相当于影像区域内GPS 数量密度的150 倍.该地区的GPS台站多为钢筋混凝土结构，在SAR 影像获取期间内具有较为稳定的雷达波散射特性，因此位于SAR 影像覆盖范围内的26个GPS台站(见图 4)多数被识别为PS单元，有5个台站由于受周围地物附加散射相位干扰较为严重，未能被识别为有效的PS单元.
将识别出的所有地面PS目标(包括21个GPS台站)进行邻域相位差分干涉，构建该地区的形变监测网.图 6 为采用Delaunay 方法构建的三角形网络，为增强网络的图形结构强度，在此基础上，采用距离限制原则增加了邻近PS 之间的连接基线，形成距离限制条件下的PS自由连接网，这使得PS连接关系的数量极大增加，这有助于提高PS 监测网的形变参数估算精度.
|
图 6 PS与GPS联合构建地表形变监测网 Fig. 6 The generated deformation monitoring network from PS and GPS stations |
图 6显示了该地区的天然PS目标(绿色圆点)与21 个GPS 测站(红色三角形)构成的形变监测网，该网络中的每一条连接基线表示了邻近PS 干涉相位的差分关系.基于前述的时序复相关目标函数解算基线差分参数，在搜索高程修正差的解空间中，顾及SRTM 高程误差约为10~20m, 因此引入高程修正差参数解空间的搜索范围为(-40 m, 40m)，高程修正差的增量(即解的步距)为0.5 m, 以在解空间中能够搜索到可靠的高程修正差数值，获得邻近PS 之间的差分参数，即形变速度差与高程修正差.基于稀疏测站的GPS水平位移采用东西向和南北向位移分量分别内插的方法，获得研究区域密集的水平二维形变场，然后使用前述的雷达LOS形变投影转换模型将LOS方向的形变速度差转换为垂直方向的形变速度差.
在求解出网络中每条基线的形变速度差数值后，由于存在大量多余观测值，使得形变速度差的闭合差不为零，为消除闭合差问题并求解出每个目标上的形变速度量，采用参数平差方法进行全局网络的联合解算[14， 17].此外，为消除干涉相位中残留的空间相关相位分量，并验证PS 目标形变速度的解算精度，将GPS在垂直方向的形变观测值作为全局网络平差的控制约束与精度检核.形变测量网络经平差处理后，得到各PS单元的形变速度量，如图 7所示.
|
图 7 研究地区PS的垂直方向形变速度 Fig. 7 The vertical deformation velocity at PS in the study area |
为验证PS形变速度的求解精度，以GPS 在垂直方向的形变观测量作为参考数据，计算PS 网络解算的形变速度值与GPS观测值之间的差异.选取其中均匀分布的7个GPS点作为检核点，其余点作为网络的约束控制点，得到各检查点的形变速度差异值，如表 3所示.
|
|
表 3 PS-InSAR网络与GPS检查点的垂直形变差异 Table 3 Comparison between GPS- and PS-derived vertical deformation velocitiesGPS PS-InSAR |
从表 3 中的形变数据可以看出，PS-InSAR 与GPS联网方法得到的垂直沉降量与GPS 测量得到的沉降量具有较高的一致性，两者的形变速度差在±2 mm 以内，采用如下均方根误差(Root Mean SquareError, RMSE)指标，评价两者的差异程度:
|
(11) |
式中，mv_def 为垂直形变测量的均方根误差，θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δHGPS与PS-InSAR网络为采用PS-InSAR 与GPS 联合建网得到的垂直形变值，θ 为雷达波侧视角，B⊥k 为空间垂直基线分量，δHGPS 为GPS网络测量得到的垂直形变值，n为检查点的个数.采用上述7个检查点数据进行计算，得到垂直形变的均方根误差mv_def =±1.5mm, 以GPS观测值作为参考数据，则PS-InSAR监测地表垂直形变的精度约为±2mm.
为进一步考察该地区地表在垂直方向上的形变空间分布，基于前述PS目标得到的形变速度值，采用内插方法获得该地区密集的垂直形变场，如图 8所示.
|
图 8 台湾西南屏东高雄地区的垂直形变场 Fig. 8 The coloured contour of vertical deformation field over SW Taiwan area |
从图 8 显示的垂直形变场空间分布来看，在1995-1999年期间，从屏东平原南部到高雄地区分别呈现出较为显著的地表沉降与地面抬升，形变量级从-15mm/a到+10mm/a.屏东平原南部(紫色区域)出现较为严重的地面沉降，主要包括林边、东港、新园等乡镇地区，沉降量从南部的-15 mm/a到北部的-5 mm/a, 越到南部沿海的林边地区，沉降量越大，而往北部地区，沉降量逐渐减小，到了屏东平原北部，其沉降量并不十分显著.此外，在西南部的高雄地区，反而出现了显著的地面抬升，抬升量从+5mm/a到+10 mm/a, 越到沿海地区，形变量越大，而到北部地区，形变量相对较小.据台湾工业技术研究院在该地区同期的水准测量成果表明[15， 22]，屏东平原南部呈现出显著的地面沉降，年均沉降速度为8~20 mm/a, 平原北部和高雄地区年均抬升速度约为14mm/a, 本文解算成果与水准观测形变量基本一致.
根据已有相关研究资料表明[5， 15， 20-21]，该地区南部的林边、东港出现显著的地面沉降，主要是由于该地区地下水的大量抽取而形成.近年来该地区水产业大规模发展，工业用水逐渐增多，对地下承压水的开采量呈急剧增加的趋势，抽取地下水作为地表水的补充资源，地下水的过量抽取导致其地下水位下降，原含水层中由于水被抽出而变成松散沉积物，在重力作用下，沉积物被逐渐压实，从而导致该地区发生显著的地面沉降.结合前述的GPS水平位移测量与PS-InSAR 垂直形变场可以发现，该地区呈现显著的地表三维形变，一方面来自于地下水抽取导致的地面沉降，另一方面由于台湾岛受两大板块(EUP与PSP)的相互挤压构造运动而导致的显著西向水平位移.
屏东平原北部和高雄地区呈现显著的地表抬升，可能的主要原因在于平原北部和高雄靠近台南台地，该地区属于地学推测的地表形变最前缘[2， 15， 22]，台南台地分别跨越西倾的后甲里断层(HouchialiFault)和东倾的台南断层(TainanFault)，这些断层形成了台南台地向上拱(PopUp)的地质构造，上升区域与台南台地的地形相符，台南台地这种上拱的特殊地质构造特征引发该地区显著的地表抬升.已有InSAR 观测数据表明[15， 22]，台南台地中央平均每年以大约12.5mm/a的速度沿雷达视线方向抬升，已施测的水准测量数据表明[15， 22]，该地区地表向上抬升速度约为14mm/a, 本文解算成果与已有地质构造分析和水准测量数据吻合性较好.
6 结论本文顾及GPS 具有高精度水平位移测量、PS-InSAR 具有高密度PS 分布的互补性优势，使用雷达LOS方向形变与GPS观测量之间的形变投影几何转换关系，联合GPS 测站与PS 目标共同构建区域形变监测网，基于参数平差方法求解区域地表三维形变场.
以地质构造活动极其活跃的台湾岛及其西南屏东高雄地区为例，基于26 个GPS 测站的水平位移测量与PS-InSAR 垂直位移监测，探测出该地区完整的三维地表形变场:由于板块构造挤压运动，造成该地区中央和西部沿海发生显著的向西位移，年均水平位移量为30~50mm/a, 高雄沿海地区发生明显的逆时针西偏南的水平位移；由于沿海地区地下水的大量开采，导致屏东平原南部呈现显著的地面沉降，从南部的-15mm/a到北部的-5mm/a, 而在西南部的高雄地区，出现了从+5mm/a到+10mm/a的地面抬升，到沿海地区抬升量有逐渐增大的趋势.
联合高精度的GPS 水平位移测量与高密度的PS-InSAR 垂直位移监测，共同构建区域地表形变监测网，一方面解决了仅仅采用稀疏GPS台站监测地表形变存在测量密度过小的问题，另一方面通过地面GPS精密观测值约束和控制雷达PS 监测网，极大地改善了PS-InSAR 技术监测地表形变的精度与可靠性.两种方法相互融合，形成技术优势互补，可为地球物理研究和形变灾害监测提供一种行之有效的三维地表形变空间观测新途径.
| [1] | Pritchard M E, Simons M A. Satellite geodetic survey of large-scale deformation of volcanic centres in the central Andes. Nature , 2002, 418(6894): 167-171. DOI:10.1038/nature00872 |
| [2] | Hu J C, Yu S B, Angelier J, et al. Active deformation of Taiwan from GPS measurements and numerical simulations. Journal of Geophysical Research , 2001, 106(B2): 2265-2280. DOI:10.1029/2000JB900196 |
| [3] | Ding X L, Liu G X, Li Z W, et al. Ground subsidence monitoring in Hong Kong with satellite SAR interferometry. Photogrammetry Engineering & Remote Sensing , 2004, 70(10): 1151-1156. |
| [4] | 单新建, 张国宏. 孕震区震前D-InSAR干涉形变场动态演化图像分析. 地震地质 , 2006, 28(3): 441–446. Shan X J, Zhang G H. An analysis of dynamic evolution of preseismic interferometric deformation fields in seismic Area. Seismology and Geology (in Chinese) , 2006, 28(3): 441-446. |
| [5] | Hu J C, Chu H T, Hou C S, et al. The contribution to tectonic subsidence by groundwater abstraction in the Pingtung area, southwestern Taiwan as determined by GPS measurements. Quaternary International , 2006, 147(1): 62-69. DOI:10.1016/j.quaint.2005.09.007 |
| [6] | 张勤, 赵超英, 丁晓利, 等. 利用GPS与InSAR研究西安现今地面沉降与地裂缝时空演化特征. 地球物理学报 , 2009, 52(5): 1214–1222. Zhang Q, Zhao C Y, Ding X L, et al. Research on recent characteristics of spatio-temporal evolution and mechanism of Xi'an land subsidence and ground fissure by using GPS and InSAR Techniques. Chinese Journal of Geophysics (in Chinese) , 2009, 52(5): 1214-1222. |
| [7] | Ge L L. Integration of GPS and radar interferometry. GPS Solutions , 2003, 7(1): 52-54. DOI:10.1007/s10291-003-0048-4 |
| [8] | Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2001, 39(1): 8-19. DOI:10.1109/36.898661 |
| [9] | Ferretti A, Prati C, Rocca F. Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2002, 38(5): 2202-2212. |
| [10] | Mora O, Mallorqui J J, Broqustas A. Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images. IEEE Transactions on Geoscience and Remote Sensing , 2003, 41(10): 2243-2253. DOI:10.1109/TGRS.2003.814657 |
| [11] | 王艳, 廖明生, 李德仁, 等. 利用长时间序列相干目标获取地面沉降场. 地球物理学报 , 2007, 50(2): 598–604. Wang Y, Liao M S, Li D R, et al. Subsidence velocity retrieval from long-term coherent targets in radar interferometric stacks. Chinese J. Geophys. (in Chinese) , 2007, 50(2): 598-604. |
| [12] | 邢学敏, 丁晓利, 朱建军, 等. CRInSAR与PS-InSAR联合探测区域线性沉降研究. 地球物理学报 , 2011, 54(5): 1193–1204. Xing X M, Ding X L, Zhu J J, et al. Detecting the regional linear subsidence based on CRInSAR and PSInSAR integration. Chinese J. Geophys. (in Chinese) , 2011, 54(5): 1193-1204. |
| [13] | 汤益先, 张红, 王超. 基于永久散射体雷达干涉测量的苏州地区沉降研究. 自然科学进展 , 2006, 16(8): 97–102. Tang Y X, Zhang H, Wang C. Subsidence retrieval in Suzhou city using Permanent Scatterers (PS) InSAR. Progress in Natural Science (in Chinese) , 2006, 16(8): 97-102. |
| [14] | 陈强, 丁晓利, 刘国祥, 等. 雷达干涉PS网络的基线识别与解算方法. 地球物理学报 , 2009, 52(9): 2229–2236. Chen Q, Ding X L, Liu G X, et al. Baseline recognition and parameter estimation of persistent-scatterer Network in Radar Interferometry. Chinese J. Geophys. (in Chinese) , 2009, 52(9): 2229-2236. |
| [15] | Huang M H, Hu J C, Hsieh C S, et al. A growing structure near the deformation front in SW Taiwan as deduced from SAR interferometry and geodetic observation. Geophysical Research Letters , 2006, 33(12): L12305. DOI:10.1029/2005GL025613 |
| [16] | Scharroo R, Visser P. Precise orbit determination and gravity field improvement for the ERS Satellites. Journal of Geophysical Research , 1998, 103(C4): 8113-8127. DOI:10.1029/97JC03179 |
| [17] | Liu G X, Buckley S M, Ding X L, et al. Estimating spatiotemporal ground deformation with improved persistent-scatterer radar interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2009, 47(9): 3209-3219. DOI:10.1109/TGRS.2009.2028797 |
| [18] | 陈强, 丁晓利, 刘国祥. PS-DInSAR公共主影像的优化选取. 测绘学报 , 2007, 36(4): 395–399. Chen Q, Ding X L, Liu G X. Method for optimum selection of common master acquisition for PS-DInSAR. Acta Geodaetical et Cartographica Sinica (in Chinese) , 2007, 36(4): 395-399. |
| [19] | Kampes B M, Hanssen R F, Perski Z. Radar Interferometry with Public Domain Tools. //Proceedings of FRINGE, 2003, December 1–5, Frascati, Italy. |
| [20] | Hu J C, Hou C S, Shen L C, et al. Fault activity and lateral extrusion inferred from velocity field revealed by GPS measurements in the Pingtung area of southwestern Taiwan. Journal of Asian Earth Sciences , 2007, 31(3): 287-302. DOI:10.1016/j.jseaes.2006.07.020 |
| [21] | Chang C P, Chang T Y, Wang C T, et al. Land-surface deformation corresponding to seasonal ground-water fluctuation determining by SAR interferometry in the SW Taiwan. Mathematics and Computers in Simulation , 2004, 67(4-5): 351-359. DOI:10.1016/j.matcom.2004.06.003 |
| [22] | Tung H. Analysis of surface deformation based on PS-InSAR technique: Case studies in coastal plainSW Taiwan. Taiwan: Department of Geosciences of National Taiwan University, 2008. |