沿海沉降变化GNSS定位及GNSS-IR组合监测

王笑蕾，牛紫瑾，何秀凤，李润川

河海大学地球科学与工程学院，江苏 南京 211100

在地面沉降与全球海平面上升的共同作用下，部分沿海地区相对海平面上升速率达到平均海平面上升速率的2倍，导致沿海地区洪水频发、湿地丧失、盐水入侵、海岸线侵蚀、人员伤亡、经济损失等，因此对沿海进行沉降监测十分必要[1]。目前最常使用的沉降监测方法是全球导航卫星系统(global navigation satellite system，GNSS)定位技术。GNSS定位获取的数据包括天线相位中心的三维坐标序列等，其中垂直分量可以反映较为基底下方的沉积物厚度的变化，包括冰川均衡前凹凸塌陷、生长断层运动和地下水抽取、矿区开采等人为运动[2-4]。GNSS测站都安装在建筑物顶部或桅杆上，基地一般深入地下2～30 m，因此GNSS定位结果中的垂直坐标分量反映的是测站基底下方的沉降变化。而由于海洋冲积效应，沿海地区基底上方的沉积物容易快速堆积，使得沿岸的沉积厚度出现较大变化。这类沉积物形成于地质时代的最新阶段—全新世(距今约11 700 a)。如果只根据GNSS定位结果的垂直坐标分量推算整体沉降，则会忽略基底上方沉降变化信息，导致低估了整体沉降和相对海平面上升值。因此，需要一种技术能够实现对基底上方沉降变化的监测，与GNSS定位技术相互配合，以实现整体沉降监测，对沿海安全及分析提供更好的保障。传统的基底上方沉降监测主要依靠杆表面高程表-标记层(rod surface elevation table marker horizon，RSET-MH)方法测量，该方法是采集沿海地表高程变化和垂直增生数据计算得到基底上方沉降速率[5-6]。由于该技术监测站点造价高，测站较少，使得基底上方沉降监测的空间分辨受限。随着GNSS的发展，GNSS-interferometric reflectometry (GNSS-IR)技术被证明可以进行反射面高度变化监测。GNSS-IR技术只基于传统的大地测量型接收机即可完成高度监测，根据其信噪比(signal-to-noise ratio,SNR)中的频率特性，便可完成反射面高度变化监测，具有全自动、小成本、长期连续、框架固定等优势[7-8]。目前，GNSS-IR技术已经用于海面[9-13]、雪面[14-19]、冻土[20]、地面[21]的高度监测中。本文利用GNSS-IR技术进行基底上方沉降变化监测，同时结合GNSS定位技术进行基底下方沉降变化监测，实现对沿海整体沉降变化监测。

利用GNSS定位进行沉降监测的技术已经相对成熟。文献[22]利用连续10 a的GNSS大地高与精密水准测量数据的对比研究，精度可达厘米级。文献[4]利用GNSS精密定位获取西安市地面沉降，说明地下水抽取和地面建设是影响地下变形沉降的主要因素。文献[23]对比分析了GNSS和GRACE得到的地表垂直形变，表明GNSS和GRACE垂直形变具有较好的一致性。文献[24]利用包含GNSS在内的导航卫星数据对青岛验潮站沉降变化进行了分析，并结合沉降信息得出平均海平面的绝对上升速率。GNSS-IR高度变化监测技术则是始于1993年。文献[25]首次提出利用GNSS反射信号测量海面高度的设想。文献[14]提出了较为完善的GNSS-IR测高原理和方法，并将其应用于雪深监测中。文献[9]进一步改进了GNSS-IR测高原理和方法，并将其用于海面高度监测中。文献[20]利用GNSS-IR技术监测了冻土高度变化。文献[21]则利用GNSS观测数据初步验证了GNSS-IR具有反演地形起伏的潜力。相关研究证实了GNSS-IR技术可以测量反射面高度变化，其具备对基底上方沉降变化的监测能力。

目前，国内外已经有一些综合利用GNSS定位和GNSS-IR技术进行绝对高度变化监测的相关研究。文献[26]采用GNSS定位方法测量陶波湖岸GNSS站点的沉降速率，并使用GNSS-IR技术测量湖面水位变化，将GNSS测量的两种方式相结合，得到陶波湖绝对湖面高度的长期序列。文献[15]利用GNSS和GNSS-IR对格陵兰岛雪深进行计算，同时得到冰川移动速度场和垂直方向变化的信息。本文综合利用GNSS和GNSS-IR技术，进行沿岸沉降变化监测研究，并采用基于最小二乘的线性拟合方法对沉降速率进行拟合，获得基底下方、基底上方和整体的沉降速率；相关原理详见第2节。算例使用密西西比河三角洲的测站数据进行沉降监测，同时使用路易斯安那州建立的沿海参考监测系统(coastwide reference monitoring system，CRMS)提供的RSET-MH基底上方沉降实测数据进行对比；相关算例分析详见第2节。

1 沿海沉降监测原理

世界范围内存在成千上万的GNSS连续跟踪站，这些跟踪站的基底一般较深，以起到维持框架的作用。以GNSS跟踪站基底为界，以上为基底上方沉积厚度，以下为基底下方沉积厚度，如图1所示。

图1 沿海平原地层和GNSS测站Fig.1 The stratigraphic and the GNSS station of coastal plain

沿海地区全新世层较易产生形变，基地一般需要打到全新世层下部或更新世层；因此，一般基底上方为全新世层，而基底下方为更新世层(有时也包括全新世层下部)。下面详细介绍利用GNSS定位和GNSS-IR技术进行基底下方沉降监测和基底上方沉降监测的原理。

1.1 基底下方沉降监测原理

传统的GNSS定位利用直射信号来确定天线的三维坐标，其中的垂直分量可以反映基底下方沉积厚度变化Dbelow

Dbelow=UGNSS(t)-UGNSS(t0)

(1)

式中，UGNSS为GNSS定位结果的垂直方向分量；t0为初始时刻；t为待求结果的试验时刻，对应的基底下方沉降速率为vbelow=∂Dbelow/∂t，该速率以抬升为正，以下降为负，与RSET-MH实测结果保持方向统一。根据式(1)，可得

UGNSS(t)=UGNSS(t0)+vbelow×(t-t0)

(2)

本文使用GAMIT/GLOBK软件采用组合解算平差获得GNSS定位垂直分量UGNSS。GAMIT/GLOBK软件采用组合解算平差，GAMIT对GNSS站的观测文件、精密星历等数据进行基线解算，解算过程将部分GNSS站作为已知控制站点，其他站点作为待解算测站点，利用站间差分和星间差分消除误差的同时进行对流层、电离层等改正模型以获取测站高精度的三维坐标。GLOBK是一个卡尔曼滤波器，对处理后的近似坐标进行平差和统一约束得到更高精度的坐标序列。GAMIT/GLOBK软件在处理数据过程中参数设置和解算策略见表1。

表1 基线解算参数和解算策略Tab.1 Baseline solving parameters and strategies

将GAMIT/GLOBK解算得到的垂直分量数据UGNSS，剔除绝对值大于2倍中误差的粗差，然后对式(2)进行最小二乘参数拟合，估计获得式(2)中的vbelow参数，即可获得基底下方沉降速率。

1.2 基底上方沉降监测原理

多路径效应导致GNSS天线接收到两种信号：一种是来自卫星的直射信号，另一种是反射信号，如图1所示。因此，接收机接收的是卫星的直射信号和反射信号的干涉信号，这种干涉信息体现在观测值文件中SNR数据中。卫星信号在一次反射情况下，反射能量最强，SNR数据中的主频一般对应一次反射。

在一次多路径反射条件下(图1)，从卫星发出信号被天线捕获时，有直射信号，也有周围地表反射来的反射信号；与相应的直射信号相比，多路径信号的传播路径较长，路程增加了D=D1-D2。利用其中的几何关系，可以得到如下公式

D1=h/sine

(3)

D2=D1cos 2e

(4)

D=(D1-D2)=2hsine

(5)

式中，e为卫星高度角；h为反射面与天线中心之间的垂直距离，一般称为垂直反射高度(reflector height,RH)。因此，直射信号与反射信号之间的相位差Δφ为

Δφ=2πD/λ=4πhsine/λ

(6)

式中，λ为信号波长。相位差中隐藏了一个频率信息f

(7)

式中，ωφ为角频率。由式(9)可得

(8)

SNR由直射信号振幅Ad、反射信号振幅Am及信号相位差Δφ决定

(9)

式中，Ad表示直射信号振幅；Am表示反射信号振幅；Ac表示合成信号振幅。由于天线增益Ad值远大于Am，可通过对输入SNR数据进行趋势项拟合，剔除直射信号，利用LSP(lomb-scargle periodogram)方法对残余项进行频谱分析可以获得频率f，从而获得h[27]。为了减小粗差，本文工作剔除了频谱峰值能量小于15 W/Hz、峰值噪声比小于3.5的RH结果。GNSS基座由水泥浇灌制成，其基座高度不变，根据图1，基底上方沉积厚度变化Dabove为

Dabove=-[h(t)-h(t0)]

(10)

对应的浅层沉降速率为vabove=∂Dabove/∂t，该速率以抬升为正，以下降为负，与RSET-MH实测结果保持方向统一。则根据式(10)，可得

-h(t)=-h(t0)+vabove×(t-t0)

(11)

根据GNSS-IR技术获得-h序列，剔除绝对值大于2倍中误差的粗差，然后对式(11)进行最小二乘参数拟合，估计获得式(11)中的vabove参数，即可获得基底上方沉降速率。

1.3 整体沉降监测原理

如上文所示，利用GNSS定位和GNSS-IR技术可分别获得基底下方沉降速率和基底上方沉降速率，将二者相加，即可获得整体沉降速率。基底上方沉降速率、基底下方沉降速率以及整体沉降均以抬升为正，以下降为负，与RSET-MH实测结果保持方向统一。本文的技术路线如图2所示。

图2 沉降监测流程Fig.2 The flowchart of settlement monitoring

具体步骤如下：

(1) 获得测站的观测文件和星历文件。

(2) 利用GAMIT/GLOBK软件进行基线解算，获取GNSS测站垂直方向坐标；剔除粗差，根据垂直坐标和式(2)进行最小二乘拟合，得到基底下方沉降速率。

(3) 根据站点坐标和星历文件，解算获得高度角和方位角，对SNR序列进行LSP分析，获得频谱图中峰值对应的RH值；剔除粗差，根据-h值和式(11)进行最小二乘拟合，得到基底上方沉降速率。

(4) 将基底下方沉降速率和基底上方沉降速率相加，获得整体沉降速率。

2 算例分析

2.1 站点及数据介绍

试验区域选取密西西比河三角洲。密西西比河三角洲是世界十大三角洲之一，是由密西西比河冲入墨西哥湾时沉积产生的。河流持续冲积使得密西西比河三角洲沉积厚度增加，导致该地区洪灾频繁发生。本文在密西西比河三角洲靠近墨西哥湾的区域选取5个测站(FSHS、GRIS、BVHS、ENG5、MSIN，如图3红色三角形所示)，下载其2015—2019年期间的观测值文件和星历文件，5个测站目前接收全球定位系统(GPS)发射的GPS L1和L2信号；其中，BVHS和ENG5站观测文件有大段缺失，而FSHS、GRIS和MSIN站观测文件几乎无缺失。5个测站基线距离在54～210 km之间(最短基线为BVHS-GRIS：54.102 km，最长基线为BVHS-FSHS：209.605 km)，距离较远，且测站深入地下深度不同、全新世厚度不同，理论上，各站间的沉降变化相似性不强。利用GAMIT/GLOBK对这5个站点进行基线解算和平差，计算坐标序列，具体的解算策略详见1.1节。GNSS-IR基底上方沉降变化反演中，GPS L2信号由于信号加密的原因，经常表现不佳，因此仅使用GPS L1信号对应的SNR数据进行RH计算。计算过程中，由于BVHS和ENG5观测文件缺失，其RH序列不连续，无法获得准确的基底上方沉降速率。因此，在沉降分析中，仅使用了FSHS、GRIS和MSIN站点；3个站点的经纬度、接收机类型、天线类型、采样间隔信息等在表2中给出。

图3 站点分布Fig.3 The distribution of stations

表2 站点信息Tab.2 Stations information

沿海参考监测系统(coastwide reference monitoring system，CRMS)建设了许多RSET-MH测量站点(如图3黄色圆点所示)。RSET-MH是将垂直的不锈钢杆深入基底10～25 m，量化相对于基准底部的表面高程变化，考虑基底垂直增生的影响，从而测量基底以上基底上方垂直沉积物的堆积状况。CRMS0513、CRMS0292和CRMS0086 3个站点是分别距离FSHS、GRIS和MSIN最近的CRMS站点，收集CRMS0086、CRMS0513和CRMS0292实测基底上方沉降数据，用于对比分析。同时，文献[28]中绘制了路易斯安那州沿海平原全新世-更新世地表的地图，公开了相对于平均海平面高度的地表深度数据，根据站点所处位置插值获得对应的地表深度数据(表2)，以方便后续分析。

2.2 基底下方沉降分析

利用GAMIT/GLOBK对测站基线解算后的数据进行平差，获得FSHS、GRIS、MSIN站点的垂直方向坐标分量。FSHS、GRIS、MSIN站点垂直方向分量每日平均精度分别为6、9、11 mm，表明GAMIT/GLOBK解算结果良好。垂直坐标分量序列及最小二乘线性拟合结果如图4所示。

图4 GPS坐标垂直分量序列Fig.4 The sequences of GPS up-coordinates

图4中，红线所示的线性趋势既是垂直坐标分量的变化趋势，也是基底下方沉积厚度的变化趋势。其对应的斜率即为基底下方沉降速率。由图4(a)、(b)可以看出，FSHS、GRIS两个测站坐标明显呈下降趋势，基底下方沉降速率分别为-3.4±0.1 mm/a、-4.5±0.3 mm/a，图4(c)中MSIN测站沉降趋势不明显，基底下方沉降速率为1.8±0.4 mm/a。文献[2,4]的研究表明沉降速率取决于沉积层序的厚度、年龄和特征。FSHS、GRIS测站坐标下降趋势比MSIN大，说明前两者基底下方沉积物厚度比后者大，且沉积物的地质年龄要小，测站基底下方为易压实的、地质年代较近的沉积物。

2.3 基底上方沉降分析

在对选定的3个测站进行GPS-IR(GPS-interferometric reflectometry)反演之前要选择高度角、方位角范围，避免植被、水域、建筑物等物体的影响。根据测站周围实际情况选取，FSHS站高度角范围是5°～25°，方位角范围是32°～62°及280°～330°；GRIS站高度角范围是15°～25°，方位角范围是118°～173°及210°～270°；MSIN站高度角范围是7°～23°，方位角范围是180°～240°。通过高度角和方位角限制后绘制站点GPS信号的菲涅尔反射区域，如图5所示。距离FSHS站点由远到近的反射区域分别对应高度角15°、20°、25°;距离GRIS站点由远到近的反射区域分别对应高度角5°、10°、15°、20°、25°;距离MSIN站点由远到近的反射区域分别对应高度角7°、12°、17°、23°，图5中以不同颜色标注。

图5 有效卫星轨迹Fig.5 The trajectories of effective satellite

对选定高度角和方位角区间内的GPS L1的SNR序列进行LSP分析，获得相关频率，并根据式(4)转换到h。根据式(6)，基底上方沉积厚度变化与h呈负相关，因此-h序列的变化即基底上方沉积厚度变化，如图6所示。图6中的反演序列经过了与2.1节同样的数据预处理和线性拟合工作。

图6中，红线所示的线性趋势既是垂直坐标分量的变化趋势，也是基底上方沉积厚度的变化趋势。其对应的斜率即为基底上方沉降速率：GRIS站为-10.3±1.6 mm/a，FSHS站为-4.9±0.7 mm/a，MSIN站为-1.8±0.7 mm/a。同时选取与GPS站点距离较近的CRMS站点(图3)，将其RSET-MH实测基底上方沉降速率与GPS-IR基底上方沉降速率反演结果进行对比，相关结果如图7所示。

图6 GPS-IR反演结果Fig.6 The GPS-IR retrieval results

图7 RSET-MH与GPS-IR基底上方沉降速率结果对比Fig.7 The comparison of RSET-MH and GPS-IR

图7中，黑色代表GPS-IR结果，蓝色代表RSET-MH结果，GPS-IR和RSET-MH是通过基于最小二乘法的多元线性回归估计的趋势和误差，中间值为最优拟合值，上下误差在置信区间为95%之内。从二者对比可以看出来，在误差范围内GPS-IR和RSET-MH计算得出的基底上方沉降速率保持很好一致性，证实了GPS-IR技术具有监测基底上方沉降的能力。

如表2所示，路易斯安那州海岸的格兰德岛(GRIS)的GPS测站建在大陆架边缘近60 m厚度的全新世沉积物上，基底深度大于20 m，基底上方沉降速率为-10.3±1.6 mm/a；富兰克林(FSHS)的站点建在20 m左右的厚度上，基底深度却在1 m深的地方，基底上方沉降速率为-4.9±0.7 mm/a，与GRIS站点基底上方沉降速率相比较小；皮尔林顿(MSIN)的站点建在1.2 m的全新世厚度之上，但是基底深度在2.4 m的更新世沉积物中，基底上方沉降速率为-1.8±0.7 mm/a，小于FSHS基底上方沉降速率。根据2.1节中所述，沉降速率与沉积物年龄和厚度有关。GRIS测站基底以上是40 m厚度的全新世沉积物，FSHS站点基地以上只有1 m左右的全新世沉积厚度，因此GRIS沉降速率较大。MSIN测站基底建在地质年代相对久远的更新世沉积物中，被压实的是更新世沉积厚度，使得基底上方沉降速率小于其他两个测站。说明基底上方沉降速率与基底深度以及深入的地质年代有关，随着基底以上全新世沉积物厚度的增大而加快。

2.4 整体沉降分析

大多数沿海平原都有河流三角洲，形成于全新世中期至晚期，位于更新世之上，这里的沉积层更容易产生沉积物压实作用。利用公开数据[28]，本文通过Surfer12.0软件中克里金插值法计算出路易斯安那州沿海地区各地全新世-更新世地表深度，如图8所示，蓝色越浅的部分表示全新世沉积厚度越大，红色越深越接近200 m厚度，因此自30.5°N，93°W向东南沿海地区，全新世厚度从0～10 m增加到超过200 m，说明越靠近墨西哥海湾，全新世沉积物堆积得越厚。统计2.2节和2.3节获得的基底下方沉降和基底上方沉降速率，同时根据第1节所述原理，获得整体沉降速率，见表3。

由表3可知，对于该区域来说，基底上方沉降速率大，不可忽略。基底上方沉降速率很大程度上影响了对整体沉降速率的监测，而传统的GPS定位方法无法监测到基底上方全新世沉积物压实导致的沉降，可利用GPS-IR技术对基底上方沉降速率补充监测。同时，根据表3和图8可以发现：全新世沉积厚度、测站基底深度与站点沉降速率有关，全新世沉积厚度越大、测站基底越深，沉降速率越大。

表3 各站点沉降速率统计Tab.3 The station sedimentation rate results mm/a

3 结论与讨论

本文综合利用GNSS和GNSS-IR技术进行沿岸的沉降监测变化，获得了基底下方、基底上方和整体的沉降速率。其中，GNSS接收机测得基底下方沉降速率，反映基底以下的地质年代相对久远的全新世和更新世沉积物的压实；GNSS-IR测量的基底上方沉降速率反映基底以上较新全新世的压实作用，在误差范围内与RSET-MH结果较为一致，证实了GNSS-IR技术可以进行基底上方沉降监测。沉降主要是由于基底上方沉积物压实造成的，因此基底较深的站点基底上方沉降速率比较浅的站点的沉降速率要大，而建立在基底较浅站点的基底上方沉降不明显。整体沉降速率随着全新世厚度的增大而加快，进一步证实地质年代越新沉降速率越大。此外，整体沉降速率与海平面上升速率相当，可用来估计陆地损失和洪水易发性。因此，GNSS接收机同时具备了监测基底下方沉降和基底上方沉降的能力，无须额外投入，只需要基于现有沉降监测跟踪站，即可实现更好的沉降监测性能。

考虑到在GNSS站点建设过程中的人工压实问题，由于GNSS-IR监测的区域较大(图5)，因此监测站点附近的人工压实并不会影响该区域基底上方监测的结果；但是，由于GNSS定位是对站点基座底部基底下方沉降的监测，因此人工压实对于基底下方沉降会有影响。同时，InSAR技术也可以进行地面沉降监测，但是会遇到失相干问题；本文的相关原理和方法在一定程度上也可以与InSAR技术互补，来提高沿海地面沉降监测的精度和分辨率。