洞庭湖软土区域时序InSAR形变与环境物理参数联合估计方法

朱 珺,朱凌杰,2,邢学敏,3,张 锐,鲍 亮,张腾飞,鲍皓丹,3

1. 长沙理工大学交通运输工程学院,湖南 长沙 410114; 2. 中南大学地球科学与信息物理学院,湖南 长沙 410083; 3. 洞庭湖生态环境遥感监测湖南省重点实验室,湖南 长沙 410114

洞庭湖位于长江中游荆江南岸,跨湘、鄂两省,是我国第二大淡水湖,湖区所处的洞庭湖生态经济区作为国家级发展片区,是中部地区崛起的重要战略支撑点。洞庭湖周边水系发达,且有湘江贯穿湖泊南北,区域内土质以软土为主,主要由淤泥、淤泥质黏性土、软塑状亚黏土、亚砂土及松散状粉细砂组成[1]。由于软土含水量大、可压缩性高、强度较低、结构松散等一系列不良土质特点,在软土区域修建的基础设施更易发生危险性沉降[2-3],累计形变甚至会导致各种构筑物及基础设施严重变形和损毁,是社会生产生活的严重安全隐患。因此对该区分布密集的基础设施进行长期持续稳定性监测,并剖析环境物理因素对时序形变特征的影响具有十分重要的意义。

针对基础设施的形变监测,传统的水准测量方法费时费力,成本高昂且难以实现大尺度监测[4-5]。永久散射体差分干涉测量技术(permanent scatterer interferometric synthetic aperture radar,PSInSAR)通过提取测区内长时间保持稳定散射特性的点(permanent scatterers,PS)进行建模分析,解算形变参数,能够高时效大尺度地获取测区时序形变结果,同时其在理论上具备亚毫米级的形变反演精度,在城市大型基础设施形变监测中极具潜力[6-9]。

利用PSInSAR技术进行时序形变反演的关键是对高相干点形变分量与形变参数之间的时序函数关系进行建模。准确且符合物理规律的形变模型不仅可以合理地描述形变随时间演化的规律,提高形变估计的精度,更可反演相应的环境物理参数,为后续形变结果的解译及灾害防治提供参考[10]。目前在PSInSAR技术中,针对城市基础设施变形监测的形变建模普遍采用简单的线性模型,即假定干涉组合内形变速率随时间呈线性变化趋势。但是,软土区域内基础设施形变在时序上具有显著非线性特征。已有研究表明,温度变化对基础设施变形的影响比运行荷载或结构损伤的影响更显著[11-12]。X波段的TerraSAR数据对热膨胀有高灵敏度,特别是当建筑设施以金属成分为主时[11,13],由于温度变化导致的材料膨胀或收缩会导致InSAR时间序列位移的强烈季节性变化[14]。因此,简单的线性假设无法揭示这一实际变形规律,不但会影响变形序列的反演精度,还难以获取可供灾害防治分析的环境物理参数。如果将民用基础设施的热膨胀特性引入InSAR形变建模和解译中,将会获取到更可靠的城市基础设施变形趋势信息,为本文研究提供了启示。

针对上述问题,本文提出一种基于热膨胀双曲线模型的时序InSAR形变与环境物理参数联合估计方法。该方法在PS基线网络的相位建模环节以软土区域形变的先验模型——双曲线模型为基础模型,将热膨胀系数和降水参数融入其中,即分别建立受热膨胀和外界环境降水影响的相位分量模型。用上述模型取代原线性速率模型,解算时不但可以获取每个PS点的变形值,还可同时获取模型中相应的环境物理参数,实现洞庭湖区域基础设施总形变序列与环境物理参数的联合反演。获得的形变序列与环境物理参数可为基础设施的安全稳定分析和灾害防治提供可靠的数据支撑。

1 时序InSAR模型构建

1.1 PSInSAR传统线性速率模型

PS点提取时,主要是通过对覆盖同一地区的多幅SAR影像的相位和幅度信息进行统计分析,筛选出具有较高后向散射特性的目标点。为有效减弱大气延迟相位、噪声相位的影响,处理中对筛选出来的PS点构建基线网络,即将任意相邻两PS点连接为一条PS基线,然后对每一条PS基线进行时序建模。

对于PS网络中任意一条基线(分别连接PS点i和j),对这两点的相位值作差,可建立如下模型[15]

(1)

(2)

1.2 热膨胀双曲线模型

(3)

双曲线模型是一种反映时序非线性关系的函数模型,可以表征沉降速率随时间表现为先快后慢的变化特征。符合软土形变随时间演化的非线性规律,因此被广泛用于路基路面工程中的软土次固结沉降预测[17]。双曲线模型的基本形式为

(4)

(5)

式中,Tm为干涉图m的从影像获取时刻。对等号右边进行泰勒级数展开,取前两项后,式(5)可表示为[18]

(6)

(7)

(8)

(9)

式中,Tempm为研究区域在SAR影像过境时刻的环境温度;Precm为研究区域当月降水量;Th和Pr分别为热膨胀系数和降水参数,为待求模型未知参数。

(10)

式中,a、b为变形参数;Th和Pr为环境物理参数。基于式(10)可对变形和环境物理参数联合求解。

(11)

式(11)顾及了湖区软土形变随时间呈非线性变化的特征,且考虑了基础设施受热膨胀和外界环境降水的影响,具有更为明确的物理意义。

2 模型的未知参数估计

2.1 基于LAMDBA-SVD的时间维参数估计

将式(11)构建的热膨胀双曲线模型代入式(1)中,则任一基线边上两PS点相位差可写为

(12)

完成对式(12)的时间维参数估计后,为衡量基线参数估值质量,采用PS基线边的时序相关系数作为质量指标[4]。基线边的时序相关系数可表示为

(13)

2.2 基于雅可比迭代的空间维参数估计

在获取了形变参数增量的估值之后,需要以此为观测值,恢复出每个PS点的绝对参数值,这一过程即为空间维参数估计。在此,利用间接平差的方法对整个PS基线网络来实现空间维解缠[22]。在求解过程中,需要在研究区域中选取一稳定的参考点以解决观测方程的法方程系数阵秩亏问题。由于对PS基线网络构建的观测方程,其法方程的系数阵B是一个大型稀疏矩阵,如果直接利用间接平差法进行求逆操作需要耗费大量计算时间,且求解结果不可靠。因此,本文采用雅可比迭代法实现法方程系数阵的迭代求逆。雅可比迭代是求解线性方程组中一种常用的迭代方法[23],主要用于求解带有大型稀疏矩阵的线性方程组,主要包括以下几个步骤。

(1) 对观测方程的系数矩阵进行LDU分解:B=D+L+U,其中D为对角阵,L为下三角阵,而U为上三角阵。

(2) 分解后原方程可转化为

X=B0X+f

(14)

式中,B0=-D-1(L+U);f=D-1b。由于D为对角阵,对其对角线元素求倒数即可得出D-1,因此极大地降低了运算的压力。

(3) 给X初值后可采用式(15)进行迭代

(15)

在时间维和空间维参数估计完成之后,即获取了研究区所有PS点的形变参数值,对这些形变参数在时间域上积分,即得到PS点上的低通形变分量。在PS基线网络布设时,设置基线边距离均不超过800 m,以尽最大可能减弱大气延迟相位的影响。通过对残余相位进行时空滤波以分离出残余相位中残留的高通形变分量部分,再与低通形变叠加,得到整个研究区域的时序形变场。算法流程如图1所示。

图1 基于热膨胀双曲线模型的洞庭湖软土区域时序InSAR形变估计处理流程Fig.1 Processing flow of InSAR soft soil area deformation estimation algorithm based on thermal expansion hyperbolic model

3 试验研究与分析

3.1 测区概况

洞庭湖周边区域是典型的软土地区。其中洞庭湖大桥位于岳阳市北门渡口下游,是湘北干线上跨越洞庭湖口的一座特大型桥梁,于1996年12月19日动工兴建,2000年12月26日通车,运营总长为5 784.5 m,是连接君山区和岳阳楼区的城市主干快速路[24]。试验中采用的TerraSAR-X影像空间覆盖范围和研究区域范围如图2所示。岳阳是一个拥有580万人口的城市,居住区众多,交通网络发达。每逢汛期,洞庭湖周边居民区都存在着堆积变形、地表塌陷甚至边坡滑坡的安全隐患。因此,对该地区进行长期的时空变形监测,对预防交通安全隐患具有重要意义。

图2 研究区概况Fig.2 Study area overview

3.2 试验处理流程

试验收集了2011年12月至2013年4月的24幅TerraSAR-X雷达卫星遥感影像。影像的前期干涉处理利用ENVI SARScape 5.2的雷达干涉处理模块,外部高程数据为30 m分辨率的SRTM-DEM(shuttle radar topography mission digital elevation model)数据。为了保证影像原始分辨率,将距离向和方位向的多视比设置为1∶1。干涉组合采用传统PSInSAR技术的单一主影像干涉模式,主影像选择2012年9月17日的影像,参数影像具体见表1。前期干涉处理的时空基线的阈值分别设置为130 m和300 d,通过平地相位的去除、滤波、去地形相位、轨道精化等处理后,共生成23幅差分干涉图。然后利用Matlab编程实现PS点选取、PS基线网络布设、时序InSAR建模及时空参数估计。在PS点识别时,采用三重指数阈值串行的方法,即振幅离差阈值、相干性阈值和强度值三重阈值均满足要求的点为候选PS点。在实际处理时,采用基于振幅离散指数、相干系数的平均值和强度值的三重指数法来挑选高相干点,分别设置三重指数阈值为0.46、0.80和0.68[25-27],同时手动去掉了河流里少数由采砂船等造成的高强度值的PS候选点。在PS基线网络布设时采用Delaunay三角网方法。图3为筛选前后的对比图。对每一条PS基线开展时间维参数估计,利用式(13)进行基线边质量评估,设置时间相关系数阈值为0.8,并迭代删除质量差的基线和孤立PS候选点。如图3所示,经过筛选后,测区共有4420个PS点和9224条基线。对整个PS网络进行空间维参数估计后,可获取PS点上的各参数估值,反演得到低通形变分量,最后对残余相位进行时间维高通滤波和空间维低通滤波,以获取残余相位中的高通形变分量,与低通形变分量累加后即可获取高相干点上时序形变结果。

表1 所用SAR影像的基本参数

图3 PS点与PS基线迭代筛选前后对比Fig.3 Comparison of PS points and PS baselines before and after the iterative screening

3.3 试验结果与分析

图4(a)显示了测区内热膨胀系数Th的分布结果,图4(b)为PS点上热膨胀系数的定量分布统计结果。可以发现热膨胀效应对于此测区地表形变的影响非常明显。通过统计分析,58%的PS点其热膨胀系数值分布在[-0.1, 0.1] mm/℃内,86%的PS点在[-0.2, 0.2] mm/℃内,95%的PS点在[-0.3, 0.3] mm/℃内。从颜色空间分布可以看出,热膨胀系数值在不同区域存在显著差异,从北至南呈现数值降低的趋势,且数值有正有负。这一现象的原因是在热膨胀特性参数恒定的情况下(当不同的基础设施由类似的材料制成时,其热膨胀特性参数视为恒定),热膨胀系数值与建筑物的高度成正比[28-29]。通过研究区域历史谷歌影像可知,在2012—2013年期间,研究区域由北向南建筑物高度整体呈由高向低的趋势,与本文获取的热膨胀系数分布整体保持一致。数值有正有负的原因如图5所示,热膨胀导致的监测对象上某个点的总形变是水平形变和垂直形变的矢量之和,当这个点沿视线向远离SAR传感器时,该点的水平形变表现为负值。因此当建筑比较高时,其水平方向膨胀在雷达视线向分量较小,而竖向的分量则较大,使得热膨胀引起的总形变为正[28-29]。为详细分析热膨胀系数分布的特征,本文分别在研究区域内的北部、中部和南部3个典型区域开展局部分析。以图4中A、B、C几个典型区域为例,由局部放大图(图4(c)、(d))可以看出,A、B区域都以高楼为主,因此其热膨胀系数值均表现为正。其中A区内最大热膨胀系数达0.31 mm/℃,B区的热膨胀系数为0.1 mm/℃。而C区的情况则相反(4(e))。C区主要以学校、老住宅区、长条型道路为主,建筑物高度相对较低,水平方向膨胀在雷达视线向负向的分量,远远大于其竖直向的热膨胀抬升,使得到的热膨胀参值以负值为主,热膨胀系数最小为-0.23 mm/℃。

图4 热膨胀参数分布Fig.4 Distribution of thermal expansion parameters

图5 单点热膨胀形变的几何图示Fig.5 Geometric illustration of single point thermal expansion deformation

图6(a)、图6(b)分别为双曲线形状参数va和vb的结果。由图6(a)—(b)可知,双曲线参数va在整体空间分布上并无显著趋势,其原因主要是va为双曲线形状参数,与各PS点空间分布并不相关;而vb为直线的斜率,与各PS点的沉降量相关,因此在空间上有明显的趋势性。经统计发现,va、vb为负值的占比分别为78%和84%,表明整个测区的形变以沉降趋势为主。图6(c)为降水参数Pr分布图。由图6(c)可知,图中靠近水域的1、2区域降水参数偏小,都在1.6×10-3左右,而处在内陆的3区域则达到了2.8×10-3,呈现出从水域附近到内陆逐渐增大的趋势。图6(d)为获取的高程改正结果,主要集中分布在-10～10 m。

图6 模型参数估计结果Fig.6 Results of model parameter estimation

需要特别指出的是,由于洞庭湖大桥并非修建在软土地基之上,不能利用本文的热膨胀双曲线模型。因此,对大桥区域的点进行掩膜处理,将这一部分点采用周期模型作为基础模型,替换式(7)中的双曲线模型分量部分进行形变建模[30],如式(16)所示

(16)

图7为获取的研究区域2011年12月28日至2013年4月3日的时间序列形变结果。从空间分布上来看,沉降较为明显的区域一般分布在洞庭湖东岸(D区域),由湖岸向内部城区沉降区域逐渐减少。洞庭湖东岸沿线一带沉降较为明显,最大累积沉降达13.3 mm,而城区中部区域沉降量较小,最大累积沉降仅为4.7 mm。从时序变化上看,研究区域沉降速率在时序上呈现出明显的先快后慢的特点,符合软土沉降的非线性规律。2011年12月28日至2012年11月11日,洞庭湖东岸沿线一带沉降趋于稳定,而从2012年11月11日至2013年4月3日,沉降趋势逐渐增大,其中最大沉降发生在2013年4月3日,局部(D区内)累积最大沉降超过15 mm。除此之外,测区右下角(E区域)从2013年1月5日开始逐渐下沉,直到2013年4月3日,最大沉降量达到12 mm。通过查找历史资料,发现测区右下角在2013年前后开展了厂房、建筑物拆除等基建工作,工程沉降或许是导致此处沉降的主要原因。试验中还探测到了少部分抬升区域,如东风湖沿湖区域(F区域)分布了部分抬升点,抬升量主要集中在2～6 mm。由于该区位于元江凹陷与阜山隆起之间的地质构造边缘地带[1],因此推断该抬升与岳阳市所处的地理位置,江汉洞庭盆地的地质构造、阜山隆起的影响有关。

图7 研究区的时序形变Fig.7 Time-series deformation of the study area

为了更清晰地显示变形的时序演化过程,分别提取了两个特征点开展时序形变分析。如图8所示,PS1位于岳阳市巴陵广场瞻岳门城楼楼顶东南角处,PS2位于岳阳市岳阳楼区建设北路与外滩一路交叉口以南50 m处。图9为两个特征点上各形变分量的时间序列演化情况。PS点上各形变分量的占比如表2所示。可以看出,无论是PS1还是PS2,总形变均以双曲线分量为主,占比分别为81.5%和78.6%;热膨胀及降水相关形变分量均小于2 mm。PS1的累积沉降量达9.5 mm,发生在2013年1月16日,而PS2累积沉降达7.4 mm。

表2 PS点上各形变分量占比

图8 PS点位置及实地照片Fig.8 PS points location and field photos

图9 PS1、PS2各形变分量时序形变结果Fig.9 Time-series deformation of each component on PS1 and PS2

分别提取了洞庭湖大桥的3个特征点(见图8中PS3、PS4和PS5),其时间序列形变如图10所示。由于PS3与PS5位于河两岸的桥墩附近,形变结果接近,相比于位于桥梁中间的PS4形变整体浮动偏大。由图10中可以看出,在3个特征点上的总时序形变都呈现出明显的周期性变化,其中周期形变分量都达到了80%以上,最大可达-13 mm,而热膨胀分量形变都在±2.5 mm之内。

图10 PS3、PS4和PS5各形变分量时序结果Fig.10 Time-series deformation of each component on PS3, PS4 and PS5

根据文献[13],热膨胀特性参数能够反映观测对象的物理性质,用简单的函数估计桥梁材料的热膨胀特性参数,其表达式为

K=Th/H

(17)

式中,K为热膨胀特性参数;H表示桥梁主梁相对于桥墩底部的高度;热膨胀系数Th可由本文算法获取。根据本文的估算,桥梁材料的平均热膨胀特性参数为12.1×10-6/℃。根据收集的现场施工设计资料,洞庭湖大桥主梁采用C50预应力混凝土连续肋板结构。该混凝土的热膨胀特性参数约为6×10-6～13×10-6/℃。因此,本文的估算结果与桥梁材料的物理性质有很好的一致性[31]。对该桥3个特征点的观测结果显示,在观测周期内,该桥累积最大沉降值为8 mm,发生在2012年8月26日。表明在整个观测期间,整个桥梁基本处于稳定状态。

为了揭示热膨胀分量与环境温度之间的相关性,分别对3个特征点上的热膨胀形变分量开展了相关性分析。图11为各特征点热膨胀形变分量与环境气温的相关性分析结果。由图11可以看出热膨胀形变分量与环境气温之间具有明显的相关性。通过计算可得,3个特征点与温度的平均相关系数分别为0.92、0.93和0.91。

图11 热膨胀分量与环境气温Fig.11 The thermal expansion-related deformation component and the air temperature

表3 热膨胀双曲线模型参数的不确定度

图12 模型参数敏感性分析结果Fig.12 Sensitivity analysis results of model parameters

4 精度分析与讨论

由于研究区域没有可用的外部形变测量数据,为了验证算法的精度,分别利用残余相位和时间相关系数评估本文方法的可靠性。根据文献[33]所述,残余相位的大小反映了模型与真实形变的拟合程度,即可以用来表征模型的建模精度。干涉图中残余相位越小,表征建模精度越高。因式(7)相位模型是基于每一条PS基线建立,需将每一基线边的残余相位恢复至每一个PS点上,再分别对23个干涉对中所有PS点的残余相位进行标准差(standard deviation,STD)的计算,再与传统线性速率模型估计的结果进行对比。对比结果如图13所示。由图13可以明显看出,热膨胀双曲线模型的残余相位明显小于线性模型,说明热膨胀双曲线模型在建模本测区的软土时序形变中更加适用。热膨胀双曲线模型在所有干涉对上的残余相位均小于0.60 rad,STD为0.40 rad,而线性模型残余相位的STD为0.63 rad,提升了36.5%。

图13 各干涉图的残余相位对比Fig.13 Residual phase comparison of each interferogram

由文献[34—35]可知,时间相关系数也可以作为模型建模精度的又一指标。对应PS点上的时间相关系数越高,则表明其建模精度越高。因此,本文计算了所有PS点在不同模型下的平均时间相关系数,结果如图14所示,可以明显看出热膨胀双曲线模型的平均时间相关系数(整个图都是深黄色)整体显著高于线性模型结果。热膨胀双曲线模型在所有PS点上的平均时间相干系数STD为0.92,而线性模型为0.63。

图14 模型的平均时间相关系数对比Fig.14 Comparison of the average temporal coherence of the models

根据前述试验结果及精度分析,现针对本文方法进行讨论如下。

(1) 本文提出将双曲线模型作为PSInSAR处理基础模型,并融入热膨胀和降水参数分量,取代原线性速率模型,在理论上顾及了软土形变随时间呈非线性变化的物理特征、湖区基础设施的热膨胀效应和外界环境降水的影响因素。相比于传统线性速率模型,本文方法残余相位的显著提升也证实了方法的可靠性。本文方法可利用InSAR相位方程组估计出各PS点的热膨胀系数值,由此可以计算出基础设施的热膨胀特性参数,提供了一种利用InSAR相位观测值获取材料热膨胀参数的手段。

(2) 由于研究区域内存在着时空异质性,不同区域、不同时期软黏土区分布面积,土质成分等特征并不完全相同。例如,洞庭湖大桥并不是全部修建在软土地基之上,利用本文的模型建模仍存在局限性。因此在试验处理时,将反映周期变化特征的周期模型作为基础模型,用于洞庭湖大桥上的点的形变建模,以更好地建模桥梁随热膨胀变化的周期特性。试验结果表明周期模型能更好地反映大桥这一特殊设施的实际形变随时间演化的规律。而本文的双曲线模型对建设在软黏土区域、高程相对较低的基础设施(如路基路面)则更有优势。

(3) 试验中获取的数据集为2011年12月至2013年4月,由于过境时间的限制未跨越两年的完整周期,且在数据集期间没有相匹配的地面水准测量数据作为外部验证数据。对于基础设施的时序形变特征的研究,如果可以用长达两年以上的数据进行分析则更具说服力。除利用残余相位来验证本文模型建模效果的提升之外,为弥补外部水准数据的缺失,搜集了洞庭湖区域早期沉降历史资料以进行验证。通过文献[36]得知,在研究时段内,洞庭湖区沉降量在1 cm以内,这与本文研究得到的沉降量一致。

5 结 论

本文提出一种软土区域时序InSAR形变与环境物理参数联合估计方法。根据软土区的变形特性,选择双曲线模型作为基础形变模型,并引入热膨胀和降水参数分量,改进了传统的PSInSAR线性速率形变模型,并可同时获取区域特征点的热膨胀系数,提供了一种获取材料热膨胀参数的遥感手段。将洞庭湖区域基础设施分布密集区作为试验区域,设计了包含PS点选点、基线网络构建、InSAR时序建模、时空维参数估计、时序总形变的生成这一完整算法流程。试验利用2011年12月至2013年4月的24幅TerraSAR-X雷达卫星遥感影像,获取了测区热膨胀双曲线模型的各参数值,并反演得到了测区最终的时序形变场。试验结果表明,截至2013年4月3日,洞庭湖东岸沿线一带沉降较为明显,最大累积沉降达13.3 mm,而城区中部区域沉降量较小,最大累积沉降仅为4.7 mm。根据生成的热膨胀系数图,估计了洞庭湖大桥的热膨胀特性参数值,与桥梁混凝土的物理性质一致。为了验证结果的可靠性,分析了23个干涉对的残余相位,得出本文方法的残余相位相比传统线性速率模型提升了36.5%。通过搜集洞庭湖区域早期沉降历史资料得知,在研究时段内,洞庭湖城区沉降量总体分布在1 cm以内,与本文研究得到的沉降量级一致。