地面沉降灾害是一个全球性问题,仅在美国,已经有遍及45个州的土地受到地面沉降的影响.传统基于水准测量等手段的地面沉降监测技术需要大量的外业工作,强度大、作业周期长、费用高,覆盖范围和采样密度均受到很大的限制,难以适用于大范围的地表形变调查和监测.进入21世纪以来,随着雷达卫星观测能力的增强和合成孔径雷达干涉测量技术(InSAR)的发展,星载InSAR技术已被广泛应用于城市地表沉降监测(王艳等,2007;张勤等,2009;Hu et al., 2018).
本文研究区域美国南加州地区地壳运动活跃,地表形变模式复杂多样,局部地区地表沉降、抬升和季节性形变等现象并存.众多学者利用InSAR技术对该地区地表形变开展了大量研究工作,Ferretti等(2000)将永久散射体InSAR技术(PS-InSAR)用于监测Pomona-Ontario盆地形变;Bawden等(2001)监测到Santa Ana盆地季节性形变;Zhang等(2012)提出临时相干点InSAR技术(TCP-InSAR),利用1995年10月至2000年12月的32景ERS-1/2数据得到的监测结果与GPS结果较为一致;Hu等(2013)提出加权最小二乘方法,利用18景ENVISAT ASAR数据得到2003年至2007年间洛杉矶地区的地表形变.上述应用均采用单一平台卫星数据,但受其重访周期和一维形变测量能力的限制,很难反映真实的地表形变特征及其演化规律.随着越来越多的SAR卫星系统发射运行,使得相同时间段内覆盖同一区域的多源多轨道InSAR数据的融合成为可能,不仅可以提高形变时间序列的连续性,还可以重建三维形变场(Wang et al., 2018).然而,要融合多源InSAR数据,除了要建立不同平台观测数据准确的函数模型之外,还需要确定观测数据的随机模型,进而根据观测数据质量确定各数据源对应的权重.但是,由于时空失相干、地形误差、轨道误差、解缠误差、大气延迟误差等因素的影响,导致不同平台通常具有不同误差统计特性,所以很难甚至不可能准确获取先验方差.邓琳等(2016)建立约束模型并结合多源数据提取南加州地区长时间序列形变结果,但该方法没有考虑到不同平台之间误差分布的差异性,而且忽略水平方向形变,结果仅提取出了垂直方向的形变.
本文针对时间和空间上均重叠的多源多轨道序列InSAR数据,根据先验形变趋势信息改进小基线集技术(SBAS)形变信号数学模型,并采用方差分量估计(Variance Component Estimation,VCE)方法(张勤等,2011)对多源InSAR观测值进行分组,通过解算出来的单位权中误差对观测值的方差(或权重)进行迭代估计更新,最终得到形变参数的最优估值.实验结果表明,本文方法得到的自适应融合形变测量结果在垂直方向的形变时间序列和水平方向的形变速率与位于研究区域内的9个GPS观测站点的监测资料均较为一致.1融合多源InSAR数据获取时间序列和三维形变场
1.1 改进的时间序列形变信号数学模型首先回顾单一平台获取时间序列形变信号的方法,并提出改进的形变信号数学模型.
假设研究区域在时间序列t0到tN上以重复轨道观测方式获取了N+1幅SAR影像,按照SBAS方法可组合生成M幅差分干涉图,所有干涉图均已经过解缠处理,对于某一像素,其解缠干涉相位δφT=[δφ1, …, δφM]与对应的SAR影像相位φT=[φ1, …, φN](定义t0时刻的相位φ0=0)之间的关系可以表示为(Berardino et al., 2002):
(1) |
其中,A为M×N维系数矩阵,每行对应一幅干涉图,其主影像和辅影像对应系数分别为1和-1,其他均为0.如果M≥N,A的秩等于N,由最小二乘原则:φ=(AT A)-1ATδφ.
为了得到一个物理意义上平滑连续的解,将公式(1)中的未知参数用平均形变相位速率来替代(Berardino et al., 2002),即:
(2) |
公式(1)经过变换可以表示为:
(3) |
其中,B是M×N维系数矩阵,其位置(j, k)的数值可以表示为B(j, k)=tk+1-tk,ISj+1≤k≤IEj,j=1, …, M,IE表示主影像,IS表示辅影像;否则B(j, k)=0.通常,采用奇异值分解(Singular Value Decomposition, SVD)方法计算B的伪逆矩阵(Berardino et al., 2002).
公式(1)还可以表示为
(4) |
其中,d(n)表示第n次观测对应的视线方向累计形变量,n=1, …, N,R为传感器到目标的斜距,θ为入射角,B⊥k为垂直基线,Δh为高程误差,Δφk为噪声相位.
该式与公式(3)联合可以写成:
(5) |
其中,
(6) |
(7) |
上述形变速率为视线向(Line-of-Sight,LOS)的一维形变反演结果,通过分析南加州的GPS监测资料可发现这些站点的位移在垂直向上波动较为明显,且具有一定周期性,而在水平方向(东西向和南北向)上仅呈现线性趋势.为此建立相应的形变信号数学模型,同时引入东西向和南北向形变速率参数,以获取三维形变测量结果.假设[Se Sn Su]为LOS向在东西向、南北向和垂直向上的投影矢量,则公式(5)可以表示成:
(8) |
其中,
(9) |
(10) |
(11) |
(12) |
其中,v″是维数为(N+3)×1的未知参数矢量,vUi,vE,vN分别为以相位记录的垂直向、东西向和南北向形变速率,i=1, …, N,D是维数为M×(N+3)的系数矩阵,δφ是维数为M×1的观测值矢量,tbm为第m幅干涉图的时间基线,m=1, …, M.
根据上述模型,假设观测值的方差(权重P)已知,采用最小二乘估计可得:
(13) |
直接利用公式(13)进行单一卫星平台数据求解三维形变参数时会导致严重的秩亏问题,因此需要融合多源InSAR数据进行多维时序形变估计.
1.2 验后方差分量估计要融合多源InSAR数据,不仅需要建立不同平台观测数据的函数模型,还需要确定观测数据的随机模型.然而如前所述,InSAR观测值的方差是难以精确获取的.因此本文利用方差分量估计对观测值的方差或权阵进行验后估计.方差分量估计的基本理论是:首先对各观测量初始定权,并采用最小二乘方法进行预平差;然后依据一定原则,利用平差得到的观测值改正数来迭代估计观测量的方差,直至各观测量的单位权中误差相等(张勤等,2011).多类观测值要根据它们的统计特性分成不同的组,由于不同的InSAR平台的观测值相互独立,我们将一类InSAR数据观测量分为一组.本文采用的实验数据为三种InSAR平台获取的观测数据,所以将观测值分为三组.设每一组为Li,其权重为Pi,其中i=1, 2, 3.首先利用最小二乘方法得到第一次估值:
(14) |
其中,Di(i=1, 2, 3)为三组观测方程系数矩阵,Vi(i=1, 2, 3)为三组观测方程改正数矢量,N=
(15) |
其中,
S可以表示为:
(16) |
其中,ni(i=1, 2, 3)为各组观测值Li(i=1, 2, 3)中观测量的个数,tr(*)表示矩阵求迹.
然后,利用单位权中误差计算新的权阵:
(17) |
将新的权阵估计值
(18) |
当迭代终止时,得到的
由于不同的SAR卫星系统具有不同的成像几何,导致其获取形变测量结果的坐标系和参考基准不一致,因此需要对雷达坐标系下的解缠相位进行地理编码,使得多源平台的观测值位于统一的地理坐标系下.为了实现参考基准的统一,不同平台的参考点均选在同一位置.在重叠区域内,选取坐标相近的测量点作为同名点,提取其解缠相位,建立1.1节所述形变信号数学模型,按照1.2节提出的验后方差分量估计方法进行解算融合处理.数据处理流程如图 1所示.
美国南加州地区由于受到太平洋板块和北美板块的南北向挤压,成为世界上构造运动最为活跃的区域之一.该区域地震多发,有文献研究记录的强震包括:1933年MW6.4的Long Beach地震、1971年MW6.7的San Fernando地震、1987年MW6.0的Whittier Narrows地震、1994年MW6.7的Northridg地震和1999年MW7.1的Hector Mine地震(Bawden et al., 2001; Davis et al., 1989; Hauksson, 1987; Hauksson and Jones, 1989; Hauksson et al., 1995; Mellors et al., 2004; Shaw et al., 1999).同时,该地区地质构造复杂,分布了许多断层,如Newport-Inglewood断层、San Andreas断层、Chino断层、Raymond断层等(Wason et al., 2002).
该区域人口密集,地下水开采强度大,如Pomona-Ontari盆地(Ferretti et al., 2000)、San Bernardino盆地(Lu and Danskin, 2001)为典型的地下水开采区,同时矿产资源丰富,油田众多,如Inglewood油田、Santa Fe Springs油田和Wilmington油田(许文斌等,2012).由于地下水抽取与回灌和石油的开采,同时受到断层影响,导致局部地区同时出现地表沉降、抬升以及季节周期性形变现象.位于该区域内的南加州综合GPS观测网(SCIGN)是世界上站点分布最密集的GPS观测网之一,但其较低的空间采样率(≥10 km)使得无法获取广域范围高空间分辨率的连续地表形变场特征信息.而InSAR技术因具有监测范围大、测量精度高、不受天气影响等特点,使其在地面沉降监测领域具有独特的优势.
2.2 数据集本实验研究共收集到22景L波段ALOS PALSAR升轨数据(时间跨度为2006年12月至2010年10月),47景C波段ENVISAT ASAR升轨数据(时间跨度2005年1月至2010年9月),以及51景ASAR降轨数据(时间跨度为2005年5月至2010年9月),三组数据集的参数如表 1所示,覆盖范围如图 2所示.从表 1和图 2中我们可以看出,三个数据集在时间和空间上均有一定的重叠度,这为提高时间序列连续性和反演三维形变场提供了数据保证.
我们首先采用StaMPS/SBAS方法对三组数据分别进行处理(时间基线和空间基线分布如图 3所示),30 m分辨率的SRTM DEM用于去除地形相位和平地效应,参考点均选择位于GPS网的参考基准站ELSC站处,获取其各自的形变速率估计和干涉图解缠结果并进行地理编码,然后在重叠区域内,在地理坐标系下识别同名点,并提取同名点在不同数据集得到所有干涉图中的解缠相位值,建立形变观测模型.针对该模型,利用方差分量估计方法得到验后方差,进行迭代处理,直到获取最优形变参数.位于三个数据集重叠区域内的9个GPS监测站数据从南加州综合GPS网(SCIGN, http://www.scign.org)下载,用以和多源InSAR数据解算结果进行比对验证.
由三组时间序列InSAR数据分别解算得到各自观测时间段内的雷达视线方向(LOS)年平均形变速率如图 4所示.
从三个数据集的解算结果可以看出,三种数据可识别到的形变区域基本一致,形变量级大致相同:在位于油气田和城市区域的主要沉降区,由石油开采和地下水抽取和回灌引起的年平均形变量最大可达20 mm;Newport-Inglewood断层两侧表现为明显不均匀形变(如图 4中标注所示).其中ASAR升降轨数据表现出较好的一致性,因为其具有相同的波长和相近的入射角;而PALSAR数据与ASAR数据存在较小的差异,由于ASAR数据相对于PALSAR数据波长较短,对于地表植被穿透能力较弱,因此ASAR数据对时间失相干更为敏感,导致在山区等植被茂密的地区均无法识别到足够数量的测量点目标,另外ASAR升降轨数据集获取时间较为接近,个别数据时间间隔仅为一天,而PALSAR数据观测时间间隔较大,且起始自2006年,也是导致其形变速率存在微小差异的原因.
3.2 时间序列和三维形变反演结果上述年形变速率测量结果均为雷达视线方向,我们可以根据此结果获知形变区的位置和形变量级,但受到单一平台重访周期等因素的限制,很难得到时间分辨率较高的时间序列结果,而且仅仅从一维形变也很难确定该研究区域的形变特征,针对以上难点,我们采用验后方差分量估计的方法将三个序列的InSAR形变测量结果进行融合.为了验证本文验后方差分量估计方法的优势,采用等权估计方法作为对照组.等权估计方法不考虑多源观测数据之间的不同误差统计特性,对各类观测值赋以相同的权重,解算过程不需要迭代.方差分量估计方法和等权估计方法得到的垂直向上InSAR时间序列形变解算结果与位于重叠区域内的9个GPS监测站观测的时间序列形变进行对比,结果如图 5所示.
从图 5中可以看出,等权估计方法明显受到不同平台数据的误差影响,解算的形变时间序列结果波动较大,且与方差分量估计方法存在偏差,说明本文采用的方差分量估计方法可以有效自适应调节不同平台观测数据的权重,经过迭代处理,最终获取形变参数的最优估值.多源InSAR数据融合后,时间节点增加至120个,提高了形变时间序列的连续性,与GPS监测数据对比来看,较为准确地反映了地表形变波动特征,且周期性也比较一致,说明本文采用的形变信号数学模型对于重建垂直向形变时间序列具有较好效果.
图 6a—b为东西向和南北向形变速率结果,为了便于比较,GPS水平方向上监测数据经过拟合得到的形变速率在图中用带有颜色的三角形表示,从多源InSAR数据解算结果以及GPS监测结果可以看出,该区域存在水平形变,而且不可忽略,东西向主要表现为向西形变(向东为正值,向西为负值),年平均速率为20~40 mm·a-1,南北向主要表现为向北形变(向北为正值,向南为负值),年平均速率为10~30 mm·a-1,但在Santa Fe Springs油田区域,存在方向相反的水平形变,这与该油田的开采导致形成漏斗状的形变有关(图 6中剖线A-A′提取形变如图 7a—b所示).另外在Newport-Inglewood断层两侧,也存在明显不均匀水平形变(图 6中B-B′剖线提取形变如图 7c—d所示).水平方向上形变速率方向和量级分布如图 8所示,可以看到该区域水平形变以西北方向为主,说明该区域受到太平洋板块和北美洲板块的相互作用,产生了向西北方向形变的趋势.
多源InSAR数据融合解算结果与GPS监测资料在水平方向上的差异统计结果列于表 2,东西向和南北向的均方根误差分别为7.567 mm·a-1和10.294 mm·a-1,可以看到两个方向的形变测量精度大致相当.虽然InSAR对南北向形变不敏感,但由于本文中使用的数据包含升降轨数据,而且水平形变具有明显线性特征,因此建立的形变信号数学模型对南北向形变也起到约束作用,有助于改善南北向的形变测量精度.
由于受到SAR卫星成像几何的影响,InSAR测量对南北方向形变不敏感,而要得到三维形变结果,至少需要包括升降轨在内的三个独立的InSAR数据集,然而目前普遍采用的多源InSAR数据融合方法均为针对大尺度形变监测设计,不适用于缓慢形变情况,或者忽略南北向形变甚至水平形变,只解算垂直向上的形变.因此在建立准确的函数模型和随机模型的基础上发展更为有效的多源InSAR数据自适应融合方法是本文的研究重点.
首先我们改进小基线集技术时序InSAR分析采用的形变反演模型,加入东西向和南北向的水平形变速率参数,从而得到垂直向上形变时间序列结果和水平向上形变速率结果,虽然该形变信号数学模型不具有普遍性,但可以利用研究区域的先验信息,发掘研究目标的形变特征,建立模型时做适当取舍,这样就会避免忽略有用形变信息,造成误判.另外,多源数据融合要考虑到不同平台数据之间误差分布的差异性,采用方差分量估计方法将多源InSAR数据进行分组,对验后方差进行迭代计算,用以精确确定权函数,最终得到形变参数的最优估值,实现多源InSAR数据融合.
5 总结本文获取覆盖南加州地区的22景升轨ALOS PALSAR数据,47景升轨和51景降轨ENVISAT ASAR数据,结合南加州综合GPS网观测数据,发现该区域在垂直向上呈现明显形变波动,且带有一定周期性,而在东西向和南北向则表现为线性形变.因此,本文借鉴小基线集技术获取时间序列的方法,根据GPS观测数据提供的先验信息,建立相应的形变信号数学模型,以得到垂直向上的形变时间序列和水平向上的形变速率结果.另外,要融合多源InSAR数据,除了要建立不同平台之间准确的函数模型,还需要确定观测数据的随机模型,进而用以获取精确的权函数,但是考虑到多源InSAR数据通常具有不同的误差统计特性,很难确定其先验方差,因此本文采用验后方差分量估计方法,对观测值的方差(或权重)进行迭代估计,最终获得形变参数最优估值,从而实现提高时间序列的连续性,重建三维形变场.与位于重叠区域的9个GPS站的监测数据对比来看:垂直向上,经过多源InSAR数据融合之后,时间节点增加至120个,较为准确地反映了地表形变波动,且周期性也比较一致,说明本文采用的形变信号数学模型对于恢复垂直向形变时间序列具有较好效果;水平方向上的形变速率结果显示,该区域存在水平形变,而且不可忽略,东西向表现为向西形变,年平均速率为20~40 mm·a-1,南北向表现为向北形变,年平均速率为10~30 mm·a-1,东西向和南北向的形变速率与GPS数据拟合得到的形变速率均方根误差分别为7.567 mm·a-1和10.294 mm·a-1,也证明本文使用的包含升降轨在内的多源InSAR数据以及建立的形变信号数学模型对南北向形变起到一定约束作用,有助于改善南北向的形变测量精度.
致谢 感谢美国地质调查局提供的数字高程模型(SRTM DEM),日本宇航局(JAXA)通过ALOS RA6项目提供的ALOS PALSAR数据(PI 3248),欧洲航天局(ESA)通过龙计划项目提供的ENVISAT ASAR数据(id 32278),美国南加州综合GPS网(SCIGN)提供的GPS观测数据.
Bawden G W, Thatcher W, Stein R S, et al. 2001. Tectonic contraction across Los Angeles after removal of groundwater pumping effects. Nature, 412(6849): 812-815. DOI:10.1038/35090558 |
Berardino P, Fornaro G, Lanari R, et al. 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing, 40(11): 2375-2383. DOI:10.1109/TGRS.2002.803792 |
Davis T L, Namson J, Yerkes R F. 1989. A cross section of the Los Angeles area:Seismically active fold and thrust belt, the 1987 Whittier Narrows Earthquake, and earthquake hazard. Journal of Geophysical Research:Solid Earth, 94(B7): 9644-9664. DOI:10.1029/JB094iB07p09644 |
Deng L, Liu G X, Zhang R, et al. 2016. A multi-platform MC-SBAS method for extracting long-term ground deformation. Acta Geodaetica et Cartographica (in Chinese), 45(2): 213-223. |
Ferretti A, Prati C, Roeea F. 2000. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 38(5): 2202-2212. DOI:10.1109/36.868878 |
Hauksson E. 1987. Seismotectonics of the Newport-Inglewood fault zone in the Los Angeles Basin, southern California. Bulletin of the Seismological Society of America, 77(2): 539-561. |
Hauksson E, Jones L M. 1989. The 1987 Whittier Narrows earthquake sequence in Los Angeles, southern California:Seismological and tectonic analysis. Journal of Geophysical Research:Solid Earth, 94(B7): 9569-9589. DOI:10.1029/JB094iB07p09569 |
Hauksson E, Jones L M, Hutton K. 1995. The 1994 Northridge earthquake sequence in California:Seismological and tectonic aspects. Journal of Geophysical Research:Solid Earth, 100(B7): 12335-12355. DOI:10.1029/95JB00865 |
Hu J, Li Z W, Ding X L, et al. 2013. Spatial-temporal surface deformation of Los Angeles over 2003-2007 from weighted least squares DInSAR. International Journal of Applied Earth Observation and Geoinformation, 21: 484-492. DOI:10.1016/j.jag.2012.07.007 |
Hu Q, Ma R H, Chen F L, et al. 2018. Urban landscape monitoring based on high-resolution spaceborne TerraSAR-X data:a case study of Nanjing City, China. Remote Sensing Letters, 9(5): 447-456. |
Lu Z, Danskin W R. 2001. InSAR analysis of natural recharge to define structure of a ground-water basin, San Bernardino, California. Geophysical Research Letters, 28(13): 2661-2664. DOI:10.1029/2000GL012753 |
Mellors R J, Magistrale H, Earle P, et al. 2004. Comparison of four moderate-size earthquakes in southern California using seismology and InSAR. Bulletin of the Seismological Society of America, 94(6): 2004-2014. DOI:10.1785/0120020219 |
Shaw J H, Shearer P M. 1999. An elusive blind thrust fault beneath metropolitan Los Angeles. Science, 283(5407): 1516-1518. DOI:10.1126/science.283.5407.1516 |
Wang Y, Liao M S, Li D R, et al. 2007. Subsidence velocity retrieval from long-term coherent targets in radar interferometric stacks. Chinese Journal of Geophysics (in Chinese), 50(2): 598-604. |
Wang Z W, Yu S W, Tao Q X, et al. 2018. A method of monitoring three-dimensional ground displacement in mining areas by integrating multiple InSAR methods. International Journal of Remote Sensing, 39(4): 1199-1219. DOI:10.1080/01431161.2017.1399473 |
Wason K M, Bock Y, Sandwell D T. 2002. Satellite interferometric observations of displacements associated with seasonal groundwater in the Los Angeles basin. Journal of Geophysical Research, 107(B4): 2074. |
Xu W B, Li Z W, Ding X L, et al. 2012. Application of small baseline subsets D-InSAR technology to estimate the time series land deformation and aquifer storage coefficients of Los Angeles area. Chinese Journal of Geophysics (in Chinese), 55(2): 452-461. DOI:10.6038/j.issn.0001-5733.2012.02.009 |
Zhang L, Lu Z, Ding X L, et al. 2012. Mapping ground surface deformation using temporarily coherent point SAR interferometry:Application to Los Angeles Basin. Remote Sensing of Environment, 117: 429-439. DOI:10.1016/j.rse.2011.10.020 |
Zhang Q, Zhao C Y, Ding X L, et al. 2009. 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), 52(5): 1214-1222. DOI:10.3969/j.issn.0001-5733.2009.05.010 |
Zhang Q, Zhang J Q, Yue D J, et al. 2010. Advanced Theory and Application of Surveying Data (in Chinese). Beijing: Surveying and Mapping Press.
|
邓琳, 刘国祥, 张瑞, 等. 2016. 多平台MC-SBAS长时序建模与形变提取方法. 测绘学报, 45(2): 213-223. |
王艳, 廖明生, 李德仁, 等. 2007. 利用长时间序列相干目标获取地面沉降场. 地球物理学报, 50(2): 598-604. DOI:10.3321/j.issn:0001-5733.2007.02.034 |
许文斌, 李志伟, 丁晓利, 等. 2012. 利用InSAR短基线技术估计洛杉矶地区的地表时序形变和含水层参数. 地球物理学报, 55(2): 452-461. DOI:10.6038/j.issn.0001-5733.2012.02.009 |
张勤, 赵超英, 丁晓利, 等. 2009. 利用GPS与InSAR研究西安现今地面沉降与地裂缝时空演化特征. 地球物理学报, 52(5): 1214-1222. DOI:10.3969/j.issn.0001-5733.2009.05.010 |
张勤, 张菊清, 岳东杰, 等. 2011. 近代测量数据处理与应用. 北京: 测绘出版社.
|