2. 中南大学地球科学与信息物理学院, 湖南 长沙 410083;
3. 洞庭湖生态环境遥感监测湖南省重点实验室, 湖南 长沙 410114
2. School of Geosciences and Info-physics, Central South University, Changsha 410083, China;
3. Hunan Key Laboratory of Remote Sensing of Ecological Environment in Dongting Lake Area, Changsha 410114, China
洞庭湖位于长江中游荆江南岸，跨湘、鄂两省，是我国第二大淡水湖，湖区所处的洞庭湖生态经济区作为国家级发展片区，是中部地区崛起的重要战略支撑点。洞庭湖周边水系发达，且有湘江贯穿湖泊南北，区域内土质以软土为主，主要由淤泥、淤泥质黏性土、软塑状亚黏土、亚砂土及松散状粉细砂组成[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=
(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为超级主影像获取时刻。此时，若令 
(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点上的各参数估值, 反演得到低通形变分量，最后对残余相位进行时间维高通滤波和空间维低通滤波，以获取残余相位中的高通形变分量，与低通形变分量累加后即可获取高相干点上时序形变结果。
| 编号 | 获取日期 | 时间基线/d | 空间基线/m |
| 1 | 2011-12-28 | -264 | -95 |
| 2 | 2012-02-10 | -220 | -26 |
| 3 | 2012-03-03 | -198 | -100 |
| 4 | 2012-03-25 | -176 | -17 |
| 5 | 2012-04-16 | -154 | 35 |
| 6 | 2012-05-08 | -132 | 58 |
| 7 | 2012-05-30 | -110 | 58 |
| 8 | 2012-08-26 | -22 | 61 |
| 9 | 2012-09-17 | 0 | 0 |
| 10 | 2012-10-09 | 22 | -145 |
| 11 | 2012-10-31 | 44 | 16 |
| 12 | 2012-11-11 | 55 | 52 |
| 13 | 2012-11-22 | 66 | -32 |
| 14 | 2012-12-14 | 88 | 144 |
| 15 | 2012-12-25 | 99 | 3 |
| 16 | 2013-01-05 | 110 | -71 |
| 17 | 2013-01-16 | 121 | -38 |
| 18 | 2013-01-27 | 132 | -94 |
| 19 | 2013-02-07 | 143 | 70 |
| 20 | 2013-02-18 | 154 | -47 |
| 21 | 2013-03-01 | 165 | -206 |
| 22 | 2013-03-12 | 176 | -120 |
| 23 | 2013-03-23 | 187 | 133 |
| 24 | 2013-04-03 | 198 | -125 |
|
| 图 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 |
| 点号 | 双曲线分量 | 热膨胀分量 | 降水分量 | 高通分量 |
| PS1 | 81.5 | 8.4 | 4.6 | 2.9 |
| PS2 | 78.6 | 8.2 | 2.7 | 3.8 |
分别提取了洞庭湖大桥的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范围。模型的不确定度是根据函数
| 参数 | 参数模拟值 | 参数误差范围 | 不确定度/(%) |
| a | -1.5 | -3~0 | -100~100 |
| b | -0.3 | -0.44~0.16 | -47~47 |
| Th | 0.1 | 0~0.2 | -100~100 |
| Pr | 2×10-3 | 0.2×10-3~3.8×10-3 | -90~90 |
|
| 图 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以内，与本文研究得到的沉降量级一致。
| [1] |
孙立群, 张鑫, 梁杏, 等. 洞庭盆地资水三角洲地区第四纪沉积环境判别及其特征[J]. 地球科学, 2021, 46(9): 3245-3257. SUN Liqun, ZHANG Xin, LIANG Xing, et al. Identification and characteristics of the sedimentary environment since the quaternary in Zi River delta, Dongting Basin[J]. Earth Science, 2021, 46(9): 3245-3257. |
| [2] |
薛祥, 宋连亮, 贾亮, 等. 高速公路软土路基工后沉降预测的新方法[J]. 岩土工程学报, 2011, 33(S1): 132-137. XUE Xiang, SONG Lianliang, JIA Liang, et al. New prediction method for post-construction settlement of soft-soil roadbed of expressway[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(S1): 132-137. |
| [3] |
SHI Xuguo, LIAO Mingsheng, WANG Teng, et al. Expressway deformation mapping using high-resolution TerraSAR-X images[J]. Remote Sensing Letters, 2014, 5(2): 194-203. DOI:10.1080/2150704X.2014.891774 |
| [4] |
ZHOU LÜ, ZHAO Yizhan, ZHU Zilin, et al. Spatial and temporal evolution of surface subsidence in Tianjin from 2015 to 2020 based on SBAS-InSAR technology[J]. Journal of Geodesy and Geoinformation Science, 2022, 5(1): 60-72. |
| [5] |
贺跃光, 蒋裕, 肖亮, 等. 基于SBAS-InSAR的钻井水溶矿山地表时序沉降场反演[J]. 长沙理工大学学报(自然科学版), 2019, 16(4): 42-47, 56. HE Yueguang, JIANG Yu, XIAO Liang, et al. Time series of subsidence inversion on drilling water-soluble mining area using SBAS-InSAR[J]. Journal of Changsha University of Science and Technology (Natural Science), 2019, 16(4): 42-47, 56. |
| [6] |
FERRETTI A, PRATI C, ROCCA F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2202-2212. DOI:10.1109/36.868878 |
| [7] |
CROSETTO M, MONSERRAT O, CUEVAS-GONZÁLEZ M, et al. Persistent scatterer interferometry: a review[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 115: 78-89. DOI:10.1016/j.isprsjprs.2015.10.011 |
| [8] |
秦晓琼, 杨梦诗, 王寒梅, 等. 高分辨率PS-InSAR在轨道交通形变特征探测中的应用[J]. 测绘学报, 2016, 45(6): 713-721. QIN Xiaoqiong, YANG Mengshi, WANG Hanmei, et al. Application of high-resolution PS-InSAR in deformation characteristics probe of urban rail transit[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(6): 713-721. DOI:10.11947/j.AGCS.2016.20150440 |
| [9] |
YANG Qin, KE Yinghai, ZHANG Dongyi, et al. Multi-scale analysis of the relationship between land subsidence and buildings: a case study in an eastern Beijing urban area using the PS-InSAR technique[J]. Remote Sensing, 2018, 10(7): 1006. DOI:10.3390/rs10071006 |
| [10] |
张永红, 吴宏安, 孙广通. 时间序列InSAR技术中的形变模型研究[J]. 测绘学报, 2012, 41(6): 864-869, 876. ZHANG Yonghong, WU Hongan, SUN Guangtong. Deformation model of time series interferometric SAR techniques[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6): 864-869, 876. |
| [11] |
FORNARO G, REALE D, VERDE S. Bridge thermal dilation monitoring with millimeter sensitivity via multidimensional SAR imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(4): 677-681. DOI:10.1109/LGRS.2012.2218214 |
| [12] |
ZHU Jinsong, MENG Qingling. Effective and fine analysis for temperature effect of bridges in natural environments[J]. Journal of Bridge Engineering, 2017, 22(6): 04017017. DOI:10.1061/(ASCE)BE.1943-5592.0001039 |
| [13] |
CROSETTO M, MONSERRAT O, CUEVAS-GONZÁLEZ M, et al. Measuring thermal expansion using X-band persistent scatterer interferometry[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 100: 84-91. DOI:10.1016/j.isprsjprs.2014.05.006 |
| [14] |
QIN Xiaoqiong, ZHANG Lu, DING Xiaoli, et al. Mapping and characterizing thermal dilation of civil infrastructures with multi-temporal X-band synthetic aperture radar interferometry[J]. Remote Sensing, 2018, 10(6): 941. DOI:10.3390/rs10060941 |
| [15] |
ZHU X, WANG Y, MONTAZERI S, et al. A review of ten-year advances of multi-baseline SAR interferometry using TerraSAR-X data[J]. Remote Sensing, 2018, 10(9): 1374. DOI:10.3390/rs10091374 |
| [16] |
HOOPER A, ZEBKER H, SEGALL P, et al. A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[J]. Geophysical Research Letters, 2004, 31(23): 1-5. |
| [17] |
杨涛, 李国维, 杨伟清. 基于双曲线法的分级填筑路堤沉降预测[J]. 岩土力学, 2004, 25(10): 1551-1554. YANG Tao, LI Guowei, YANG Weiqing. Settlement prediction of stage constructed embankment on soft groundbased on the hyperbolic method[J]. Rock and Soil Mechanics, 2004, 25(10): 1551-1554. |
| [18] |
KIM S W, WDOWINSKI S, DIXON T H, et al. Measurements and predictions of subsidence induced by soil consolidation using persistent scatterer InSAR and a hyperbolic model[J]. Geophysical Research Letters, 2010, 37(5): 87-105. |
| [19] |
MONSERRAT O, CROSETTO M, CUEVAS M, et al. The thermal expansion component of persistent scatterer interferometry observations[J]. IEEE Geoscience and Remote Sensing Letters, 2011, 8(5): 864-868. DOI:10.1109/LGRS.2011.2119463 |
| [20] |
KAMPES B M, HANSSEN R F. Ambiguity resolution for permanent scatterer interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(11): 2446-2453. DOI:10.1109/TGRS.2004.835222 |
| [21] |
邢学敏. CRInSAR与PSInSAR联合监测矿区时序地表形变研究[D]. 长沙: 中南大学, 2011. XING Xuemin. Study on monitoring the time series ground deformation in mining area based on CRInSAR and PSInSAR integration[D]. Changsha: Central South University, 2011. |
| [22] |
陈强, 刘国祥, 丁晓利, 等. 永久散射体雷达差分干涉应用于区域地表沉降探测[J]. 地球物理学报, 2007, 50(3): 737-743. CHEN Qiang, LIU Guoxiang, DING Xiaoli, et al. Radar differential interferometry based on permanent scatterers and its application to detecting regional ground subsidence[J]. Chinese Journal of Geophysics, 2007, 50(3): 737-743. |
| [23] |
张后苏. 误差理论与测量平差基础[M]. 长沙: 中南工业大学出版社, 1995: 41-45. ZHANG Housu. Error theory and surveying adjustment basis[M]. Changsha: Central South University Press, 1995: 41-45. |
| [24] |
廖建宏, 李迪清, 胡建华, 等. 岳阳洞庭湖大桥三塔斜拉桥设计[C]//中国公路学会2001年学术交流论文集. 北京: 中国公路学会, 2001: 231-237. LIAO Jianhong, LI Diqing, HU Jianhua, et al. Three-tower cable-stayed bridge design of Yueyang Dongting Lake bridge[C]//Proceedings of 2001 Academic Exchange of China Road Society. Beijing: China Road Society, 2001: 231-237. |
| [25] |
陈强, 丁晓利, 刘国祥, 等. 雷达干涉PS网络的基线识别与解算方法[J]. 地球物理学报, 2009, 52(9): 2229-2236. CHEN Qiang, DING Xiaoli, LIU Guoxiang, et al. Baseline recognition and parameter estimation of persistent-scatterer network in radar interferometry[J]. Chinese Journal of Geophysics, 2009, 52(9): 2229-2236. |
| [26] |
邢学敏, 丁晓利, 朱建军, 等. CRInSAR与PSInSAR联合探测区域线性沉降研究[J]. 地球物理学报, 2011, 54(5): 1193-1204. XING Xuemin, DING Xiaoli, ZHU Jianjun, et al. Detecting the regional linear subsidence based on CRInSAR and PSInSAR integration[J]. Chinese Journal of Geophysics, 2011, 54(5): 1193-1204. |
| [27] |
ZHU Lingjie, XING Xuemin, ZHU Yikai, et al. An advanced time-series InSAR approach based on poisson curve for soft clay highway deformation monitoring[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 7682-7698. DOI:10.1109/JSTARS.2021.3100086 |
| [28] |
MA Peifeng, LIN Hui, LAN Hengxing, et al. Multi-dimensional SAR tomography for monitoring the deformation of newly built concrete buildings[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 106: 118-128. DOI:10.1016/j.isprsjprs.2015.04.012 |
| [29] |
QIN X, LIAO M, ZHANG L, et al. Structural health and stability assessment of high-speed railways via thermal dilation mapping with time-series InSAR analysis[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2017, 10(6): 2999-3010. DOI:10.1109/JSTARS.2017.2719025 |
| [30] |
COLESANTI C, FERRETTI A, NOVALI F, et al. SAR monitoring of progressive and seasonal ground deformation using the permanent scatterers technique[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(7): 1685-1701. DOI:10.1109/TGRS.2003.813278 |
| [31] |
黄杰, 吴胜兴, 沈德建. 不同粗骨料混凝土早期热膨胀系数试验研究[J]. 结构工程师, 2010, 26(3): 154-158. HUANG Jie, WU Shengxing, SHEN Dejian. Experimental study on thermal expansion coefficient of different type coarse aggregates concrete at early ages[J]. Structural Engineers, 2010, 26(3): 154-158. |
| [32] |
YANG Zefa, LI Zhiwei, ZHU Jianjun, et al. An InSAR-based temporal probability integral method and its application for predicting mining-induced dynamic deformations and assessing progressive damage to surface buildings[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2018, 11(2): 472-484. DOI:10.1109/JSTARS.2018.2789341 |
| [33] |
ZHAO Rong, LI Zhiwei, FENG Guangcai, et al. Monitoring surface deformation over permafrost with an improved SBAS-InSAR algorithm: with emphasis on climatic factors modeling[J]. Remote Sensing of Environment, 2016, 184: 276-287. DOI:10.1016/j.rse.2016.07.019 |
| [34] |
WANG Chao, ZHANG Zhengjia, ZHANG Hong, et al. Seasonal deformation features on Qinghai-Tibet railway observed using time-series InSAR technique with high-resolution TerraSAR-X images[J]. Remote Sensing Letters, 2017, 8(1): 1-10. |
| [35] |
XING Xuemin, CHANG H C, CHEN Lifu, et al. Radar interferometry time series to investigate deformation of soft clay subgrade settlement—a case study of Lungui highway, China[J]. Remote Sensing, 2019, 11(4): 429. |
| [36] |
姜端午, 黄树春, 张苑平, 等. 基于地质环境遥感调查与监测数据探讨洞庭湖演变规律[J]. 国土资源遥感, 2010(S1): 124-129. JIANG Duanwu, HUANG Shuchun, ZHANG Yuanping, et al. A discussion on the evolution of the Dongting Lake based on geo-environmental remote sensing investigation and monitoring data[J]. Remote Sensing for Land & Resources, 2010(S1): 124-129. |