本文内容来源于《测绘学报》2023年第12期（审图号GS京(2023)2441号）

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

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

1. 长沙理工大学交通运输工程学院, 湖南 长沙 410114;

2. 中南大学地球科学与信息物理学院, 湖南 长沙 410083;

3. 洞庭湖生态环境遥感监测湖南省重点实验室, 湖南 长沙 410114

基金项目：国家自然科学基金(42374046;42330717; 42074033;51878078; 41701536; 41904003; 61701047);湖南省自然科学基金(2022JJ30589; 2019JJ50639; 2020JJ5571); 湖南省教育厅重点项目(18A148);洞庭湖生态环境遥感监测湖南省重点实验室开放课题(2021-011);湖南省交通运输厅科技创新项目(202211);长沙市杰出创新青年培养计划项目(kq2209011);中国水利水电第八工程局有限公司科研项目(2023060);湖南省自然资源厅科技计划项目(20230118CH); 长沙理工大学大学生研究性学习与创新实验计划项目

摘要：洞庭湖生态经济区是国家级重要发展片区, 是国家战略发展的重要部分, 但洞庭湖区域土质以软土为主, 必须对区域内密集分布的基础设施开展长期持续的安全稳定性监测。目前使用PSInSAR进行监测时普遍是采用线性速率模型, 但洞庭湖软土区域基础设施形变随时间演化具有明显的非线性特征, 并且还受热膨胀、降水等环境物理因素影响。针对这一问题, 本文提出一种用PSInSAR同时估计变形和环境物理参数的方法。该方法利用洞庭湖软土形变预测的双曲线模型和热膨胀效应先验模型代替原方法中的线性速率模型, 同时融入热膨胀参数和环境降水参数, 在PSInSAR解算过程中将形变与环境物理参数一并求解。在此基础上本文还提出了基于LAMDBA-SVD的时间维参数估计和基于雅可比迭代的空间维参数估计的解算策略。试验选取了覆盖岳阳市洞庭湖区的24景TerraSAR-X雷达卫星遥感影像, 设计了包含PS点选点、基线网络构建、InSAR时序建模、时空维参数估计、时序总形变的生成这一完整算法流程, 获取了此区域典型基础设施的热膨胀系数, 以及2011年12月至2013年4月的时间序列形变场。利用模型的残余相位来评估建模精度, 结果显示本文方法的残余相位为0.4 rad, 相比传统线性速率模型提升了36.5%。

关键词：合成孔径雷达干涉测量 软土 热膨胀 基础设施 形变监测

朱珺, 朱凌杰, 邢学敏, 等. 洞庭湖软土区域时序InSAR形变与环境物理参数联合估计方法[J]. 测绘学报，2023，52(12)：2127-2140. DOI: 10.11947/j.AGCS.2023.20220287

ZHU Jun, ZHU Lingjie, XING Xuemin, et al. Joint estimation method of time-series InSAR deformation and environmental physical parameters for soft clay area over Dongting lake[J]. Acta Geodaetica et Cartographica Sinica, 2023, 52(12): 2127-2140. DOI: 10.11947/j.AGCS.2023.20220287

阅读全文：引 言 洞庭湖位于长江中游荆江南岸，跨湘、鄂两省，是我国第二大淡水湖，湖区所处的洞庭湖生态经济区作为国家级发展片区，是中部地区崛起的重要战略支撑点。洞庭湖周边水系发达，且有湘江贯穿湖泊南北，区域内土质以软土为主，主要由淤泥、淤泥质黏性土、软塑状亚黏土、亚砂土及松散状粉细砂组成[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)式中，m为干涉对序号；Δϕi, jm、Δki, jm分别为相邻两PS点在第m幅干涉图中对应的相位增量和整周模糊数增量；βi, jm为高程相位改正系数，βi, jm=，Bj为空间基线，θ为影像的入射角，Ri表示雷达传感器与PS点之间的距离；ΔδHi, j为相邻两PS点高程改正差值；Δresi, jm表示残余相位，主要由大气延迟相位，非线性形变相位和传感器几何噪声等组成[16]；Δϕdefo_i, jm为相邻两PS点上的形变分量的增量。传统PSInSAR技术获取地表形变过程的关键环节就是对Δϕdefo_i, jm建模，目前在城市基础设施变形监测中，较为普遍的是采用线性模型，即认为干涉组合内，形变速率随时间呈线性变化趋势，即可表示为 (2)式中，Δvi, j为相邻两PS点的线性沉降速率的增量；Ti, j为时间基线。式(2)中未知参数为PS基线上的参数增量Δvi, j、ΔδHi, j及Δki, jm。显然，式(2)仅适用于描述那些随时间呈线性变化的形变特征。由于软土变形在时序演化上具有显著的非线性特征，因此式(2)不再适用。

1.2 热膨胀双曲线模型

软土变形在时序演化上可分为3个阶段：瞬时变形、主固结变形和次固结变形。其中，次固结变形非线性特征明显，线性形变速率假设不再适用。此外，像桥梁、建筑物这样的大型基础设施变形除与桩基土体沉降有关之外，更与钢结构表面、拉索等部件随温度变化而产生的物理形变(即热膨胀变形)有直接关系[14]，这一热膨胀效应对变形速率图有明显的影响。考虑这3个方面的因素，Δϕdefo_i, jm可表示为 (3)式中，Sdefo_i, jm为软土形变分量；Stemp_i, jm和Sprec_i, jm分别为热膨胀效应和外界环境降水引起的形变分量。双曲线模型是一种反映时序非线性关系的函数模型，可以表征沉降速率随时间表现为先快后慢的变化特征。符合软土形变随时间演化的非线性规律，因此被广泛用于路基路面工程中的软土次固结沉降预测[17]。双曲线模型的基本形式为 (4)式中，S0、t0分别为起始参考点的沉降量与观测时间；St、t分别为任意点上的沉降量及对应的时间；va、vb为待求的双曲线形状参数。由于式(3)中Sdefo_i, jm为相邻两PS点的形变增量，即Sdefo_i, jm=Sdefo_jm-Sdefo_im，则任意PS点i的形变分量Sdefo_im可以用双曲线模型进行表示 (5)式中，Tm为干涉图m的从影像获取时刻。对等号右边进行泰勒级数展开，取前两项后，式(5)可表示为[18] (6)式中，Tc为超级主影像获取时刻。此时，若令 ，则式(6)可化简为 (7)则待求的双曲线形状参数va、vb可简化为求a、b。根据文献[13, 19]所述，PS点i的热膨胀效应Stemp_im和外界环境降水引起的形变分量Smprec_i可以表示为 (8) (9)式中，Tempm为研究区域在SAR影像过境时刻的环境温度；Precm为研究区域当月降水量；Th和Pr分别为热膨胀系数和降水参数，为待求模型未知参数。以双曲线模型作为软土形变基础模型，融入热膨胀系数和降水参数后，PS点i的形变相位分量ϕdefo_im可以表示为 (10)式中，a、b为变形参数；Th和Pr为环境物理参数。基于式(10)可对变形和环境物理参数联合求解。对任意一条PS点基线，Δϕdefo_i, jm可写为 (11)式(11)顾及了湖区软土形变随时间呈非线性变化的特征，且考虑了基础设施受热膨胀和外界环境降水的影响，具有更为明确的物理意义。

2

模型的未知参数估计2.1 基于LAMDBA-SVD的时间维参数估计将式(11)构建的热膨胀双曲线模型代入式(1)中，则任一基线边上两PS点相位差可写为 (12)式中，Δki, jm、ΔδHi, j、Δai, j、Δbi, j、ΔThi, j、ΔPri, j为待求参数。求解这些未知参数增量的过程为时间维参数估计的过程。而将求解出的未知参数增量估值结果则作为已知数据，再恢复至每一个PS点的绝对未知参数的过程则为空间维参数估计。Kampes首次将整数最小二乘搜索算法(least squares ambiguity correlation adjustment，LAMDBA)成功应用于PSInSAR技术的时间维相位解缠中，可估计出每一条PS基线上的整周模糊数Δki, jm[20]。本文将这一思想用于求解式(12)中的未知参数。由于利用LAMDBA算法求解式(12)时，式(12)的法方程系数为秩亏，需要采用补充虚拟观测方程的方法来求解[21]。补充虚拟观测方程后，式(12)的方程组法方程系数矩阵仍呈现明显的病态，这会对参数估计结果有较大的影响。因此，本文对观测方程的系数矩阵采用奇异值分解(singular value decomposition，SVD)的方法，可使由补充虚拟观测导致的病态条件数减半，并且回避了对大型病态矩阵的求逆，提高了运算效率。完成对式(12)的时间维参数估计后，为衡量基线参数估值质量，采用PS基线边的时序相关系数作为质量指标[4]。基线边的时序相关系数可表示为 (13)式中，Δresi, jm可利用估计出的参数增量值代入后求得。在实际处理中，对时序相关系数设定一个阈值，去掉低于阈值的基线边及孤立PS点后，再重新进行布网、建模和时间维参数估计。如此迭代处理，直到基线边的时序相关系数阈值满足要求为止。

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

在获取了形变参数增量的估值之后，需要以此为观测值，恢复出每个PS点的绝对参数值，这一过程即为空间维参数估计。在此，利用间接平差的方法对整个PS基线网络来实现空间维解缠[22]。在求解过程中，需要在研究区域中选取一稳定的参考点以解决观测方程的法方程系数阵秩亏问题。由于对PS基线网络构建的观测方程，其法方程的系数阵B是一个大型稀疏矩阵，如果直接利用间接平差法进行求逆操作需要耗费大量计算时间，且求解结果不可靠。因此，本文采用雅可比迭代法实现法方程系数阵的迭代求逆。雅可比迭代是求解线性方程组中一种常用的迭代方法[23]，主要用于求解带有大型稀疏矩阵的线性方程组，主要包括以下几个步骤。(1) 对观测方程的系数矩阵进行LDU分解：B=D+L+U，其中D为对角阵，L为下三角阵，而U为上三角阵。(2) 分解后原方程可转化为 (14)式中，B0=-D-1(L+U)；f=D-1b。由于D为对角阵，对其对角线元素求倒数即可得出D-1，因此极大地降低了运算的压力。(3) 给X初值后可采用式(15)进行迭代 (15)(4) 式(15)的迭代终止条件为|Xi(k+1)-Xi(k)| < ε，ε的大小根据实际运算要求设置。此时求得的X(n+1)则为未知参数的最优解。在时间维和空间维参数估计完成之后，即获取了研究区所有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影像的基本参数Tab. 1 Basic parameters of the used SAR acquisitions

图 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)式中，Sdefo_im为周期模型形变分量；T1为周期性形变的周期(取为一年)；ΔAi, j和ΔBi, j为待求的周期形变速率参数。图 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。

图 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图选项

表 2 PS点上各形变分量占比Tab. 2 The proportion of each deformation component at the feature points

分别提取了洞庭湖大桥的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]，热膨胀特性参数能够反映观测对象的物理性质，用简单的函数估计桥梁材料的热膨胀特性参数，其表达式为 (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图选项

根据文献[32]所述，InSAR形变模型获得的形变精度主要取决于模型参数估计的精度。如果估计形变对某个参数高度敏感，说明这一参数的微小误差则会导致最终估计形变的较大变化。而对于那些不敏感的参数，即使引入了一定的误差，对最终估计的时序变形结果也不会造成较大影响。因此，对模型开展参数敏感性分析，不但可以揭示模型参数对形变的影响规律，还可以有效地指导InSAR数据处理，具有一定的理论和实践意义。在此，利用模拟数据对热膨胀双曲模型的参数和估计形变开展敏感性分析试验。利用8组不同的模型参数来生成最终的形变，表 3显示了4个参数的模拟值Ms和每个模型参数的误差Mr范围。模型的不确定度是根据函数×100%来定义的。敏感性分析结果如图 12所示，各参数敏感性系数可由图 12中对应直线的斜率表示。由图 12不难看出，总变形对双曲线形状参数a、热膨胀系数Th高度敏感。这表明对于本文的试验区域，这两个参数的估计误差将会对总形变的误差有较大影响。因此，准确估计a和Th值，对于后续准确预测区域基础设施的变形，预防基建工程安全事故的发生至关重要。

表 3 热膨胀双曲线模型参数的不确定度Tab. 3 Inaccuracy of parameters of thermal expansion hyperbolic model

图 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以内，与本文研究得到的沉降量级一致。作者简介第一作者简介：朱珺(1986-), 男, 博士, 讲师, 研究方向为时序InSAR形变监测、PolInSAR沙漠地区穿透深度反演。E-mail: jzhu@csust.edu.cn通信作者：邢学敏, E-mail：xuemin.xing@csu.edu.cn

初审：张 琳复审：宋启凡
终审：金 君

资讯