0.   引言

软土地区公路易发生过大沉降，影响交通安全与畅通，对其开展工后长期监测具有重要意义^[1]。传统公路沉降监测手段，如水准测量、GNSS测量等方法多基于离散点监测，且需要监测人员定期进入现场，不但成本昂贵，且无法生成公路面状沉降结果，这严重制约了其在实际工程中的应用^[2-3]。合成孔径雷达干涉测量（interferometric synthetic aperture radar, InSAR）技术作为一种新兴的雷达遥感对地观测技术，可用于提取地表固定目标在观测时段内与雷达卫星之间距离变化信息，以其监测范围大、空间分辨率高、无接触式测量等特点，大大弥补了传统测量手段的不足，为公路沉降监测提供了新的契机^[4]。然而，InSAR技术易受时空失相关影响，且无法获取高精度时序形变结果，限制了其在公路，尤其是软土地区公路长期时序沉降监测中的应用。

小基线集干涉测量技术（small baseline subset InSAR, SBAS-InSAR）是InSAR技术的延伸，是近年来发展起来的一种高级InSAR技术^[5]，可弥补InSAR技术时空相干弱的劣势，提升沉降监测的精度，已被证实在大型交通设施（铁路、公路和桥梁等）的工后沉降监测中极具潜力^[6-8]。在SBAS-InSAR技术处理中，形变建模是生成形变的重要环节，是建立形变与形变参数之间的时序函数关系的过程。然而，传统SBAS形变模型大多采用线性模型表征形变随时间的变化特征^[7-9]，模型中的系数也多为数学系数，并无明确的物理含义，忽略了观测对象的物理性质。在软土区域修建的公路，其地表沉降随时间的变化过程是较复杂的非线性过程，用假定的线性函数去拟合这一过程，显然不能合理反映其真实沉降规律，会严重影响沉降监测的精度^[10]。文献[11，12]将表征外界环境影响因素的物理参数引入冻土时序InSAR模型中，论证了在形变模型中考虑到监测对象形变规律可极大地提高监测精度。根据大量岩土力学领域学者的研究表明，软土的沉降受土的流变性影响尤其严重。流变变形表征的是软土在外届荷载恒定不变时，土体内部应变随时间发生缓慢蠕变的现象^[13-14]，是公路工后运行期主要的沉降特征。如果可以建立一种反映软土沉降蠕变特征的InSAR物理模型，可以显著提高软土地区公路沉降建模精度和后续沉降监测的精度，增强InSAR模型的物理意义，并辅助后续沉降解译和预测。

本文提出了一种基于雷达遥感对地观测技术（InSAR）的软土地区公路沉降自动化、大范围监测方法，将非线性黏弹塑性四元件组合流变模型——Maxwell-NVPB^[15]引入InSAR技术中，构建时序InSAR相位与流变参数之间的函数关系，以替换传统SBAS技术中线性模型，克服其物理意义不明确和建模精度差的局限，实现软土地区公路条带时序沉降的自动获取，为软土地区公路沉降早期识别和养护管理提供依据。通过模拟数据和广东伦桂路水准实测数据验证本文方法的可行性和可靠性，并将其应用于湖南岳阳湖区公路大范围沉降监测。

1.   方法原理

1.1   基于非线性流变模型的软土地区公路沉降模型

Maxwell-NVPB模型是一种四元件非线性黏弹塑性流变模型，可完整描述软土衰减、稳态和加速3个蠕变阶段，由于考虑了力学参数随时间的弱化规律，更接近材料蠕变真实状况。当应力不变时，其流变状态方程为^[15]

式中：$\varepsilon$为材料产生的应变；${\sigma _{\text{c}}}$为恒定的外界荷载，为一常数；${\sigma _{\text{S}}}$为应力阈值，通过开展材料长期强度试验分析获取；n为加速蠕变阶段的蠕变影响指数；E为流变材料的弹性模量；${\eta _1}$为黏度，又称黏滞系数；${\eta _2}$为塑性系数；E，${\eta _1}$和${\eta _2}$是流变学中重要的蠕变参数，也是模型中的待估参数；t为应变发生过程的总时间。

软土沉降${S_{\text{v}}}$与应变$\varepsilon$之间的函数关系可表示为^[16]

式中：$H$为软土层平均厚度；${t_1}$，${t_2}$分别为形变发生的起始日期；$h$，$t$均为积分变量。由式（1），（2）推导可得${S_{\text{v}}}$如下式所示：

式（3）为软土在垂直方向的沉降与蠕变参数$E$，${\eta _1}$和${\eta _2}$的时序函数关系，可以描述软土在工后期的流变变形特征。

1.2   基于Maxwell-NVPB沉降模型的时序InSAR建模

对干涉图中任一高相干点，其InSAR相位可以表示为^[5]

式中：$\Delta {\varphi _j}$为干涉对相位，可通过雷达卫星影像差分干涉处理直接获取；j∈(1, 2, 3, …, M)，为干涉对序号；$(x,r)$表示高相干点的SAR坐标，${S_M}$表示高相干点上的低频形变分量，即可用函数模型建模的形变分量；$\Delta \varphi _{{\text{topo}}}^j(x,r) = (4{\text{π }}{B_ \bot }/\lambda R\sin \theta )\Delta h(x,r)$为残余地形相位分量，其中${B_ \bot }$为垂直基线长度，$R$和$\theta$分别为传感器到地面固定目标的距离和雷达入射角，均可从雷达影像头文件中获取；$\Delta h$为高程改正值，$\lambda$为波长；$\Delta \varphi _{{\text{res}}}^j(x,r)$为残余相位，由噪声相位、大气延迟相位和高频形变相位组成。在传统SBAS技术中，低频形变分量${S_L}$用假定的线性形式表示：

式中：${v_j}$为干涉影像对应时段内的线性形变速率；${t_{\text{A}}}$，${t_{\text{B}}}$分别为用于干涉组合的两景影像获取时间。InSAR获取的形变${S_{{\text{LOS}}}}$为雷达视线向（LOS：line of sight）的形变，需要依据${S_{\text{v}}} = {S_{{\text{LOS}}}}/\cos \theta$（其中$\theta$为雷达入射角）转换为垂直向沉降。在此将低频形变分量${S_M}$改写为线性分量与流变沉降分量的组合，即${S_M} = {S_L} +$ ${S_{\text{v}}}$。由式（3）~（5），基于Maxwell-NVPB沉降模型的InSAR相位方程组可表示为

式中：未知参数为蠕变参数$E$，${\eta _1}$，${\eta _2}$，线性速率$v$，高程改正$\Delta h$。当观测方程个数j大于或者等于5时，即可解算出影像上所有高相干点的未知参数估值。上述未知参数的求解，是一个非线性参数估计的问题。本文采用遗传算法（genetic algorithm，GA）和单纯形搜索法（simplex method, SM）组合，对未知参数进行求解^[17]。

1.3   时序InSAR沉降监测方法技术流程

基于Maxwell-NVPB沉降模型的时序InSAR沉降监测方法共包括以下几个重要步骤：①前期SAR影像差分干涉处理，生成用于形变建模和参数估计的差分干涉相位序列；②基于Maxwell-NVPB的时序InSAR形变建模，构建如式（6）的视线向雷达相位方程组；③基于遗传算法和单纯形搜索法组合的模型未知参数估计，生成用于计算时序形变的参数估值$E$，${\eta _1}$，${\eta _2}$，$v$以及高程改正值$\Delta h$；④软土地区公路时序沉降生成；⑤与SBAS方法进行对比验证。方法流程图如图 1所示。

图  1  时序InSAR沉降监测方法技术流程图

Figure  1.  Flow chart of time-series InSAR subsidence monitoring method

下载: 全尺寸图片 幻灯片

1.4   基于模拟数据的算法验证

在InSAR数据处理中，形变参数的估计是准确求解模型未知参数的过程，是保障后续形变结果生成的必要环节。为保证在研究区域真实雷达影像处理中形变参数估计算法的顺利运行和参数的准确求解，需要预先利用模拟数据对参数估计算法（1.3节中步骤③，④）的可行性和可靠性进行验证。设置参数$E$的取值范围[0，50] MPa，${\eta _1}$[0，10]×10⁶ MPa·s，${\eta _2}$[100，300]×10⁶ MPa·s，$v$[-200，100] mm/a，各参数真实值利用二维高斯函数模型模拟生成，在模拟的真值矩阵中随机提取出100个高相干点的参数值作为真值，与1.3节中步骤③计算获取的参数估计值进行比较以验证本文算法可行性和可靠性。图 2为噪声水平0.5 rad时流变参数估值与真实值的对比结果，图 2中显示红色与蓝色折线基本重合，表明模拟试验获取的参数估计值与真实值吻合较好。未知参数均方根误差计算结果为：E为±0.61 MPa，${\eta _1}$为±0.27×10⁶ MPa·s，${\eta _2}$为±15.74×10⁶ MPa·s，v为±0.01 mm/年和$\Delta h$为±1.90 m。

图  2  参数估计值偏差与真实值对比(噪声水平为0.5 rad)

Figure  2.  Comparison between deviation of parameter estimation and the real data (noise level 0.5 rad)

下载: 全尺寸图片 幻灯片

2.   真实数据验证

2.1   研究区域

试验选取广东佛山伦桂路的部分路段对前述所提方法开展验证。研究区域位于珠江三角洲平原，是典型的软土地区。道路沿线的土壤由泥黏土、粉质黏土和泥砂组成，土质不均匀，易发生滑坡、不均匀沉降等地质灾害。研究区域地理位置见图 3所示，图中矩形框为测区范围。

图  3  伦桂路测区位置

Figure  3.  Location of Lungui Road

下载: 全尺寸图片 幻灯片

2.2   试验结果及验证

试验共采集23景覆盖研究区域的TerraSAR-X雷达影像，影像的获取时间为2015年1月1日—2017年1月18日，空间分辨率为3 m，时间分辨率为22 d，入射角为26.39°，采用HH极化方式，2016年1月10日影像为主影像。通过查阅公路设计资料，试验中软土层平均厚度H取10 m，${\sigma _{\text{c}}}$取0.25 MPa，${\sigma _{\text{s}}}$取0.24 MPa；n通过对材料进行敏感性分析获取，取值范围为（0，3.7）。图 4为本文获取的部分伦桂路测区时序形变结果，从空间分布上可发现A，B两个区域出现明显的沉降漏斗。截至2017年1月18日，研究路段最大沉降量达到113 mm。

图  4  伦桂路测区时序沉降

Figure  4.  Time-series subsidences of Lungui Road (reference date: 1 January—December 2015)

下载: 全尺寸图片 幻灯片

为验证获取的沉降结果精度，收集了测区内两个水准点的地面水准测量结果（水准点位置见图 4中K1、K2，K1位于安利特大桥附近；K2位于容桂特大桥附近)。水准测量结果与Maxwell-NVPB沉降模型获取的结果和传统线性模型获取的结果对比表示在图 5中。3组沉降序列数据总体趋势保持一致，证实了本文方法应用于软土地区的沉降监测是可行的。

图  5  与水准点形变结果进行对比

Figure  5.  Comparison of results by proposed method and levelling measurements on benchmark

下载: 全尺寸图片 幻灯片

本文方法获取的沉降与水准结果的均方根误差（RMSE：root mean square error）分别为：K1点±8.4 mm，传统线性模型均方根误差为±11.2 mm，精度提升28.5%；K2点±2.7 mm，传统线性模型为±3.2 mm，精度提升16.9%。本文方法总RMSE为±5.6 mm，相比传统线性模型提升了22.7%。K1点处与水准结果的相关系数为0.97，而传统线性模型为0.67，相比提升了33.7%。

3.   洞庭湖区域软土地基公路形变监测应用实例

3.1   研究区域概况

本文选取洞庭湖区域岳阳楼附近基础设施及居民密集区作为研究对象，开展本文方法的应用实例研究，以为当地公路管理部门的公路沉降早期识别和养护提供参考。洞庭湖位于湖南省北部长江南岸，面积辽阔，是中国五大淡水湖之一。其周边水系发达，广泛分布着冲、沉积软土层。软土层主要由淤泥、淤泥质黏土、塑性状粉质黏土、粉质砂土组成，具有含水量高、孔隙比大、压缩性强、黏聚力低、塑性系数大、抗剪强度低等特点。正是由于这些特点，洞庭湖沿岸区域易发生地表沉降，需对其开展长期沉降监测。图 6为洞庭湖测区在Google Earth上的地理位置放大图。如图可见，沿湖路是岳阳市区通往城陵矶港的交通要道之一，同时连接着洞庭湖大桥、杭瑞高速公路等重要工程基础设施。其行车宽度24 m，双向六车道，路面由沥青混凝土组成。岳阳市交通干道密集分布，每逢汛期，洞庭湖周边区域居民区都面临着堆积变形、地表坍塌等安全隐患。因此，对洞庭湖区（尤其是沿湖沿岸区域）公路进行长期的沉降监测具有重要意义。

图  6  洞庭湖测区位置

Figure  6.  Location of Dongting Lake

下载: 全尺寸图片 幻灯片

试验共收集了24景覆盖研究区域的TerraSAR-X波段Stripmap降轨影像。影像获取时间为2012年12月—2013年4月。时间分辨率为11 d，入射角为39.23°，极化方式为HH极化，主影像选取2012年11月11日影像。图 7为本文获取的洞庭湖测区2011年12月28日—2013年4月3日的部分时间序列沉降结果。如图所示，从空间分布上来看，研究区域沉降量级较大的沉降区域集中分布在洞庭湖沿湖东岸，而靠近内陆城区相对较为稳定。截至2013年4月3日，沿湖大道的最大沉降量达到46 mm，而内陆城区累计沉降量仅12 mm。沿湖大道沿线软土层主要由淤泥和淤泥质黏土组成，软土层较厚，地表水系发达且水流速快，水位和水量受到降水影响较大，因此沉降相较于内陆区域更为明显。而内陆城区土质多含粉土、粉砂，道路路基淤泥质黏土含量少，软土层厚度较小，且地表渗入水流量小，其水量、水位虽受降水影响，但城区排水设施较发达，因此沉降相对较弱^[18]。

图  7  洞庭湖测区时序沉降（起始时间：2011年12月28日）

Figure  7.  Time-series subsidences in Dongting Lake area (reference date: 28 December 2011)

下载: 全尺寸图片 幻灯片

3.2   试验结果及分析

为了定量分析这一时序沉降演化过程，将提取的两个特征点（图 7中标记分别为Y1，Y2）上的时序沉降结果与传统线性模型获取的沉降结果进行对比分析，如图 8所示。Y1点位于沿湖路与沐河园路十字路口交汇处，处于沉降量级变化较大区域，截至2013年4月3日获取的最大沉降量为29 mm，而线性模型获取的最大沉降量为23 mm。Y2点位于洞庭大道和冷水铺路十字路口交汇处，处于沉降量级变化较小区域，截至2013年4月3日获取的最大沉降量为12 mm，而线性模型获取的最大沉降量为9 mm。通过对Y1点沉降量进行分析，得出该点总沉降中，流变分量占总沉降73%，线性分量仅占总沉降25%，高通形变分量仅占3%。这表明总沉降中起主导作用的非线性沉降分量主要由流变变形引起。相比之下，传统线性模型结果中，Y1点线性沉降分量占总形变81%，残余相位中提取的非线性沉降仅占19%。这一结果表明当软土沉降呈现强烈非线性特征时，线性模型获取的结果将无法提取出那些掩盖在残余相位中真实非线性沉降信息，导致获取的最终沉降结果精度不高。Maxwell-NVPB模型能更好地反映研究区域内软土沉降趋势与蠕变参数和形变参数之间的时序演化规律。

图  8  特征点时序沉降结果对比

Figure  8.  Comparison of subsidence results at feature points

下载: 全尺寸图片 幻灯片

通过图 8中特征点的时序沉降曲线，在2012年8月26—2012年10月9日，2012年12月14日—2013年1月16日和2013年3月12日—3月23日3个时段内3点均出现了不同程度的隆起。结合图 8中湖南省岳阳市统计局提供的温度与降水曲线，2012年9月岳阳市降水量较8月增加70 mm，气温下降5℃。在降水量突增和温度持续下降的共同作用条件下，由于外部降水补给，地下含水层的水位上升，以Y1点为例，地表沉降出现3 mm回升。2012年12月底—2013年2月，降水量增加72 mm。虽然温度几乎保持恒定，但较大降水增量仍导致Y1点地表抬升3 mm。2013年3月—2013年4月，降水量增加38 mm，导致地表抬升2 mm。根据本文获取的洞庭湖区域时序形变结果，岳阳当地公路管理部门重点对沉降明显路段开展了现场监测。

3.3   精度评估

由于洞庭湖区域没有可用的外部水准测量数据，为了验证提出模型的改进效果，文献[11]中指出，残余相位$\Delta \varphi _{{\text{res}}}^j(x,r)$的均方差可衡量模型的建模精度，在此将本文方法（图 9中Maxwell-NVPB模型）和传统SBAS中线性模型对应的残余相位进行对比。图 9为获取的干涉图中对应的残余相位值。干涉图中残余相位越小，表明模型建模精度越高。结果表明，Maxwell- NVPB模型对应的残余相位明显小于传统线性模型。通过计算各高相干点上对应的两种模型残余相位均方根值（RMS：root mean square）得出传统线性模型的RMS为0.42 rad，而Maxwell-NVPB为0.24 rad，提升了41.9%。

图  9  模型残余相位对比结果

Figure  9.  Comparison of results of residual phase of models

下载: 全尺寸图片 幻灯片

4.   结论

（1）本文提出了一种基于雷达遥感对地观测（InSAR）技术的软土地区公路沉降监测方法。利用模拟数据验证本文的沉降参数估计算法时，各未知参数估计值的误差在高噪声情况下仍小于参数估值的4%，表明本文所提方法是可行且可靠的。

（2）通过广东伦桂路的真实测区试验验证可知，各流变参数值与对应区域沉降量呈反相关关系，整体呈现下沉速度由快到慢的非线性变化趋势；与外部水准实测结果对比，本文方法获取形变的精度提升了24%，证实了Maxwell-NVPB流变模型引入传统InSAR技术中，可显著提升软土地区公路沉降监测的精度。

（3）将本文方法应用于湖南岳阳湖区公路大范围沉降监测，获取了该区1.5 a的时序沉降结果；该区域公路沉降呈现先快后慢的趋势，沿湖区沉降明显大于内陆区，最大累计沉降量达46 mm；内陆区相对沉降较小，最大沉降量为12 mm；与传统SBAS线性模型获取的沉降结果对比可知，本文方法建模精度提升了41.9%。

（4）本文将非线性流变模型引入高级雷达遥感对地观测InSAR技术的建模环节，提高了沉降序列获取精度，实现了公路大范围、高精度、高效率的沉降自动获取，可为路基路面工程监测提供良好借鉴。下一步研究中，将会综合考虑公路的交通荷载作用机制，将描述车辆循环动应力的模型及参数融入现有模型，深入探讨在交通荷载作用下的软土变形规律，以进一步提高沉降监测精度。