地基GNSS 用于冻土区域地表环境多参数反演研究

李景洮，张双成,2,3，刘宁，陈雄川，王杰，冯智杰

( 1. 长安大学地质工程与测绘学院, 西安 710054；2. 地理信息工程国家重点实验室, 西安 710054；3. 西部矿产资源与地质工程教育部重点实验室, 西安 710054 )

0 引言

冻土是指温度在0℃及以下的含有孔隙冰的各种岩石和土壤，可分为季节性冻土与永久性冻土[1].从上到下主要分为活动层和永冻层，冻土冻融过程中，地表孔隙水与冰处于相变过程.冻结期孔隙水受到冻结吸力作用向活动层底部移动，遇到冰锲后形成分离冰.而随着温度周期性变化，分离冰处于周期性冻融过程，地表发生相应抬升与沉降.冻土形变受到下层分离冰状态的控制，因此出现了过渡层的概念[2].我国幅员辽阔，冻土分布区域极其广泛，尤以青藏高原为甚，周期性的监测冻土区域地表形变，对于定量评估冻土灾害规律以及采取必要性防范措施具有极其重要的意义.

冻土形变测量关键在于固定的基准点，保证其不受固体地球动力学变化影响，通常需要锚点位于永久冻土层下方.文献[3-5]制作了相对容易标记的小型原位测量仪器测量冻土形变，如升降管、铁磁探头、床架等.以上方法虽然方便简单，但受到明显的地缘范围局限，而大范围的水准测量费时费力.直至近年来以卫星遥感为代表的对地观测技术及各种数据差分处理技术的蓬勃发展，Little 等[6]首次在全球尺度使用空间差分技术，得到相对准确的美国Alaska两地点形变结果.此后大尺度地快速、精确测量地表形变成为必然趋势[7].

传统冻土区域环境监测仅局限于高精度地表形变的研究，尤以InSAR 技术的发展，为处理地表微小形变提供了简单可靠的途径，但受制于SAR 卫星重访周期与数据量大小，无法做到实时连续冻土形变监测.地基GNSS 中全球导航卫星系统干涉反射测量法(Global Navigation Satellite System Interferometric Reflectometry，GNSS-IR)技术由Larson 等[8]首次提出并验证其可行性，其基本思想是利用直射信号与反射信号产生的干涉效应来求解地表参数.其以全天时、全天候、使用穿透力较强的L 波段、可忽略云雨气候与电离层扰动影响的优势，迅速发展成为满足当前高精度、低成本、轻量化、长时序环境监测需求的一种新兴手段，对土壤湿度[8]、积雪[9]等地表参数较为敏感，研究发现其在冻土形变中具有巨大潜力.Liu 等[10]首次通过GNSS-IR 技术获得2004—2015 近十年间阿拉斯加巴罗冻土区域不同年份沉降趋势，为冻土探测提供了新的手段，Hu 等[11]通过对巴罗2004—2016 年形变分析，提出了基于stefan 方程的融合冻融季节性指数因子的复合模型.随后，对2018—2021 年巴罗地表形变进行多系统反演，提出消除GLONASS、Galileo 因重访周期所导致高程序列周期性递减的重建方法[12]，重建结果精度提升显著.文献[13-15]使用GNSS-IR 检索了加拿大北极地区与青藏高原寒区地表形变，改进解译了寒区土壤水分算法，分析并提出冻土区域GNSS-IR 与合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)协同作用.

本文通过地基GNSS 中GNSS-IR 技术与GAMIT/GLOBK 软件处理美国阿拉斯加州巴罗多年冻土区SG27 测站及附近国际GNSS 服务(International GNSS Service，IGS)站FAIR 测站2013—2021 年观测数据，分别解译了长时序地表形变、测站形变、土壤湿度及大气水汽，分析水汽与土壤湿度相关性，以期为地基GNSS 冻土区域环境监测提供新的补充.

1 研究区域及数据来源

1.1 研究区及测站概况

巴罗地区位于美国阿拉斯加州北冰洋岸最北部，为典型平原地带，地势平坦，冻土广泛分布[16].Brewer研究发现巴罗地区底部是上百米的永久冻土连续区，平均活动层厚度约30cm[17].该地区通常无雪期在7—9 月中旬，解冻时间在6—8 月中旬[18]，而2016 年属于极端异常年份，全球气温升高致使降雪从7 月初持续到10 月底，降雪迟于往年1 个月.

SG27 测站(156°36′37″W，71°19′22″N)位于美国阿拉斯加巴罗北部，隶属于美国的板块边界观测台网(plate boundary observational GNSS network，PBO)计划，该计划目前拥有1100 个高精度GPS 站点，旨在用于监测美国西部大陆板块运动.测站在WGS84坐标系下的高程为9.383m，天线高约为3.8m，接收器底部深埋在永久冻土层基岩中.测站周围环境及菲涅尔反射区如图1 所示，四周视野开阔，低矮草坪，表层为粉砂质土壤，环境适宜GNSS 信号的接收与地表参数反演，在2018 年6 月期间，测站更换POLARX5接收机，天线为TRM59800.00.

图1 SG27 测站与菲涅尔反射区

1.2 研究数据

本实验所选用数据有：1)SG27 测站2013—2021 年6—11 月GPS L1 频段观测数据、IGS 精密星历数据、广播星历数据，由于测站周围存在房屋遮挡，故选取90°~180°方位角，高度角为5°~25°的信噪比(signal to noise ratio，SNR)数据；2)为解算大气可降水量(precipitable water vapor，PWV)，选取距离SG27 约800km 的IGS 站FAIR 测站相同时段数据；3)为获取土壤湿度模型参数，选取国际土壤湿度网络(International Soil Moisture Network，ISMN)2017年有限土壤湿度实测数据.数据来源如表1 所示.

表1 研究数据来源

2 研究方法及数据处理流程

2.1 研究方法

GNSS-IR 冻土测量原理如图2 所示，该方法可有效测量出反射面至接收机天线相位中心的距离，降雪期地面反射高小于天线高，两者差值即为雪深；而无雪期测得反射高度增量Δh与地面沉降量为相反数，地表形变可由Δh计算.

图2 GNSS-IR 冻土测量原理图

GNSS-IR 利用SNR 来量化直射信号与反射信号，可表示为

式中：Ad为直射信号振幅；Am为反射信号振幅；φ为直射信号与反射信号之间的相位差，由于Ad>>Am，所以直射信号决定了混合信号SNR 的整体趋势.需利用低阶多项式拟合混合信号SNR 得到直射信号趋势项，并将其剔除获取反射信号SNR，图3 分别绘制了2018 年240 天PRN1 卫星的混合信号SNR、拟合趋势项及分离后反射信号SNR.

图3 2018年240 天的PRN1 卫星混合信号SNR、趋势项及反射信号SNR

由图3 可知，卫星高度角较低时，SNR 较低但振荡幅度明显.随着高度角增大，SNR 增大但振荡幅度逐渐减小，因此通常需选用低高度角SNR 数据.SNR 表现为正弦函数形式，反射信号SNR 与反射高度h的关系最终经过量化后可表示为

式中：RSNm为反射分量大小；A为反射信号幅值；h是接收器天线相位中心距反射面的垂直距离；λ 为载波波长；θ是高度角；φ0是相移.设t=sinθ，RSNm也可用式(3)表示：

由上式可知，反射信号SNR 振荡主频率f与天线相位中心距反射面的垂直距离h成正比关系.

L-S 谱分析是在傅里叶变化的基础上改进的针对非均匀采样间隔的数据进行分析，能够有效地从中提取弱周期信号.通过L-S 谱分析得到反射信号振荡主频率f，即可得到反射高度h.

反射高度h与天线高的差值即为雪深/地表形变，若差值为负值证明处于无雪期，地表处于形变过程.特别地，为了保证反演高度h的可靠性，本实验选取多条可靠轨道数据，筛选了多个GPS 卫星反射高，设置中值滤波0.2m，求取范围内多星平均反射高作为单天反射高.同时，使用标准差来替代反射高的不确定性.

经以上步骤，将反射信号SNR 与高度角序列迭代拟合，得到相位与振幅参数，最后通过相应土壤湿度线性模型(式(4))反演土壤湿度.

式中：CSM为土壤湿度；S为拟合斜率；∆φ为相位、振幅参数；CSMrepid为最低含水量，通常用一点土壤湿度最小值表示.

2.2 数据处理

2.2.1 GNSS-IR 地表参数解译

通过使用上述原理算法反演阿拉斯加巴罗地区2013—2021 年地表形变、土壤湿度，并反演了降雪异常年份2016 年积雪深度，将其与PBO 网站实测数据进行精度验证.具体步骤如下：1)PBO 观测数据、广播星历、精密星历数据下载；2)通过crx2rnx 工具对观测文件进行格式转换，使用gfzrnx 工具统一为采样率30s 的Rinex3 版本数据进行数据项处理；3)使用GNSS SNR 工具提取GPS SNR；4)挑选方位角90°~180°，高度角5°~25°的L1 频段的SNR 数据；5)低阶多项式分离直反射信号；6)L-S 频谱分析得到每日多星先验反射高；7)挑选多条可靠轨道求多星平均反射高，计算雪深/地表形变；8)非线性最小二乘法拟合SNR 与高度角序列，获得相位、振幅参数；9)利用土壤湿度经验模型得到土壤湿度变化.数据处理流程如图4 所示.

图4 地基GNSS 数据处理流程

2.2.2 大气水汽反演

使用GAMIT/GLOBK 软件进行高精度数据处理，具体步骤如下：1)GNSS 高精度数据处理，得到SG27 概略坐标、钟差以及对流层延迟参数等；2)计算对流层天顶总延迟(zenith total delay，ZTD)与天顶干延迟；3)由ZTD=ZHD+ZWD 计算对流层天顶湿延迟，其中ZTD 为对流层总延迟，ZHD 为天顶干延迟，ZWD 为天顶湿延迟；4)PWV=K*ZWD 计算大气可降水量，其中K为无量纲比例系数.本文为探究长时序大气水汽变化，故选取每日12:00 水汽值.数据处理流程如图4 所示.

3 异常年雪深/时序地表形变分析

3.1 异常年份雪深信息分析

2016 年为异常高温年，降雪比往年迟了近一个月，经过对该年数据处理,得到反演雪深，并将其与实测数据进行对比，如图5 所示.经过计算，反演值与实测值具有较强相关性，决定系数R2为0.8155，均方根误差(root mean squared error，RMSE)为0.0643，平均绝对误差(mean absolute error，MAE)达到0.0402，雪深偏差普遍控制在2~7cm.

图5 2016年雪深反演

在2016 年年积日45—120 天内，反演值普遍高于实测值，产生这种现象的主要原因可能有以下两点：1)与该年出现的极端高温有关，雪深每日变化差异大，在进行反射高算法计算时，由于这种差异而导致一些真实而极端的不同时段的数据被当成异常值来剔除，从而产生一些与实测数值偏差较大的误差；2)极端高温导致PBO 大部分实测数据被视为异常值剔除，导致偏差的产生与积累.

3.2 2013—2021 地表形变信息分析

测站锚点处于永久冻土层，可确保高程基准不受固体地球动力影响.在无雪期，地表经过季节性变化，活动层底部分离冰周期性冻融而导致地表周期性沉降与抬升.通常认为在有雪期，地表活动层处于饱和状态，不易发生形变.针对无雪期，绘制2013—2021年6—11 月冻土活动层形变，如图6 所示.同时，为证明其反演结果为冻土活动层形变，联测IGS 站FAIR测站解算同时段SG27 测站形变，结果如图7 所示，测站形变相对稳定，存在部分数据缺失，天顶(up，U)方向最大形变在5mm 之内，可证明反演结果为冻土活动层形变.

图6 2013—2021 年冻土活动层形变

图7 2013—2021 年SG27 测站形变

本实验以反射高度等于天线高时为无雪期起点，融化期通常开始在6 月中旬至7 月初期间，而冻结期一般发生在8 月中下旬至11 月初期间.在2013—2021 年，基本沉降趋势为大幅度沉降后抬升，部分趋势中最初抬升是由于该年雪水当量较多，导致低温时段过渡层分离冰冻结，地表短期性抬升.由于GNSS-IR无法精确监测降雪期垂直形变，因此部分抬升未完全显示趋势.2013 年为长时期沉降，年最大沉降速率达到-0.101m/a.2014 年形变趋势呈现为沉降-抬升-沉降-抬升，主要与7 月初持续低温有关，最大沉降速率为-0.075m/a，抬升速率为0.069m/a；2015 年与2019 年趋势相似，呈现线性沉降趋势后缓慢抬升，最大沉降量分别为-0.087m，-0.075m；2017 年由于设备老化(2018 年8 月更换设备)、高度角方位角条件、算法中值等原因，导致少量数据情况下，滤值提取结果趋势项较为离散，大体表现为全年具有沉降抬升趋势，沉降速率为-0.054m/a，抬升速率为0.033m/a；2016、2018、2020 年能够清晰看出垂直形变趋势，最大沉降点均发生在8 月中下旬至9 月中上旬，沉降速率分别为-0.077m/a、-0.027m/a、-0.067m/a，抬升速率分别为0.05m/a、0.025m/a、0.06m/a.其中，2016年为异常降雪年，无雪期持续至11 月，因此形变量最大，比其余两年持续跨度长约一个月.而2018 年整体形变较小，在3cm 之内变化，地表活动平缓.2021 年7 月时，该地区温度极端寒冷，致使地表活动层分离冰增大，地表抬升，出现抬升-沉降-抬升-沉降趋势变化，最大沉降量达到-0.038m.

4 大气水汽与土壤湿度时序分析

针对土壤湿度反演，由于缺少站点实测数据，故通过ISMN 获取该地区少量有限土壤湿度实测数据，使用土壤湿度经验模型反演土壤湿度，大气水汽与土壤湿度变化趋势如图8 所示.

图8 2013—2021 年无雪期大气水汽与土壤湿度变化

土壤湿度变化与降雨存在直接关系，而复杂水汽分布与降雨存在相关关系[19]，因此土壤湿度在一定程度上能显示大气水汽的相对含量变化，但只能大致表示在变化趋势而非幅度上.该区域无雪期水汽变化应符合先上升后降低的趋势，水汽在解冻期与冻结初期蒸腾作用明显并且冻结期大气干延迟较大，水汽逐渐降低.

由图8 可知，该地区2013—2021 年间水汽化符合常规趋势，水汽逐渐降低.，土壤湿度与大气水汽保持相对滞滞后性变化趋势.水汽在积累后形成降雨并滞后作用于土壤湿度，土壤湿度在0.1~0.32cm3/cm3变动，在解冻初期，由于积雪前期融化致使水汽和土壤湿度处于较高水平，随后呈现动态变化.在反演过程中，缺少实测数据致使无法验证多星反演相位精度，从而出现部分相对平缓集中的数据.2013 年土壤湿度能够较好的表示出水汽变化趋势，对于水汽急剧增大或者降低的点，土壤湿度变化均能够体现.2014 年解冻期水汽逐渐增大，土壤水分持续积累在0.24cm3/cm3之上.2015 年土壤湿度变化比水汽变化更加敏感，在8 月初期水汽达到最大，随后土壤湿度积累至最大0.242cm3/cm3后开始在冻结期下降.2016 年为全球气候异常年，水汽出现高蒸腾现象，最高可达36mm，高温使其部分数据缺失以及土壤湿度整体处于较低水平.2017 年和2018 年解冻初期雪水当量较多，土壤湿度从高逐渐降低，最终趋于0.16~0.23cm3/cm3之间，水汽大幅度增高，解冻期极为活跃，随后在冻结期降低.2019 年水汽与土壤湿度响应关系活跃，冻结期结束时水汽骤增导致土壤湿度相对2020 与2021 两年部分数据反演相位偏差较小，土壤湿度结果波动性较小，维持约在0.23cm3/cm3.

5 结束语

本文利用地基GNSS 技术对2013—2021 年冻土区阿拉斯加巴罗地区的无雪期地表参数进行了解译，得到如下结论：反演该地区异常年份2016 年积雪深度，雪深反演数据与实测数据绝对系数R2为0.8155，RMSE 为0.0643，MAE 为0.0402.雪深偏差控制在2~7cm，具有较高的反演精度.

处理了2013—2021 年期间无雪期地表形变与土壤湿度，形变结果能够较为明显的表示出地表沉降/抬升趋势，得到了每年沉降/抬升速率.同时解算同时段SG27 测站形变，其U 方向形变在5mm 之内，有效证明反演形变结果为冻土活动层形变，而非测站本身变化导致.

选取距离SG27 测站约800km 的IGS 站FAIR测站，使用GAMIT/GLOBKE 软件解算了该地区2013—2021 年无雪期大气水汽分布，分析其与土壤湿度相关性.结果显示：大气水汽与土壤湿度存在较弱滞后性对应关系，在水汽达到阈值时，土壤湿度在后一天发生增大/减小变化，大气水汽在前一天存在相应增大/减小趋势，但仅表示在趋势上而非幅度值.

地基GNSS 技术以其高精度等优点，在冻土环境监测领域存在巨大潜力，但目前仍旧存在以下现状：

1)缺少站点土壤水分实测数据，无法保证反演结果的可靠性.虽然目前GNSS-R 关于土壤水分产品较多，但受分辨率、重访周期等因素影响，缺失数据较多，同时精度可能满足不了GNSS 高精度测量结果.

2)高度与方位角条件限制致使数据量过少时，地表形变反演中少量异常值会被当成中值进行保存，从而产生误差.

针对上述问题，存在一些工作，如构建融入气象数据的机器学习插值方法，对实测土壤湿度缺失数据进行合理插值；多系统多频数据融合，构建合理定权方法，从而求解合理可靠的反射值等.

致谢：感谢美国PBO 网站提供实验数据.