GNSS-IR双频数据融合的土壤湿度反演方法

荆丽丽，杨磊，*，汉牟田，洪学宝，孙波，梁勇

（1.山东农业大学 信息科学与工程学院，泰安271019； 2.北京航空航天大学 电子信息工程学院，北京100083）

土壤水分又称土壤湿度或土壤含水量，是联系陆-气相互作用的关键物理量之一[1]。准确测定土壤水分可以有效管理灌溉、保证作物稳产高产等[2-3]，因此研究快速、实时、准确、高效、廉价的大面积农田土壤湿度监测方法是精准农业领域一个至关重要的研究方向。全球导航卫星系统反射信号干涉测量（Global Navigation Satellite Signal-Interferometer and Reflectometry，GNSS-IR）是一种新兴的地基遥感技术[4]，其利用 GNSS直射与地表反射信号的干涉效应，通过获取干涉信号所携带的反射面的相关特性，进行地物参数的遥感，如土壤湿度[5-7]、积雪厚度[8-11]等。

基于GNSS-IR的土壤湿度测量于2008年由Larson等[12]首次提出，该技术使用接收机记录的信噪比（Signal-to-Noise Ratio，SNR）数据进行土壤湿度测量，其物理本质是利用直射信号和反射信号之间的干涉效应，而接收机记录的SNR数据可以看作是对干涉信号的度量，Larson等[12]通过实验证明了SNR的振荡幅度、频率和相位3个观测量与土壤湿度相关。经过近年的发展，GNSS-IR技术已得到国内外专家学者的广泛关注。国外方面，2014年，Alonso-Arroyo等[13]改变接收天线的极化方式，提出通过跟踪垂直极化（V-Pol）和水平极化（H-Pol）SNR数据间的相位差进行土壤湿度测量的方法，该方法提高布鲁斯特角度估计的准确性，从而提高土壤水分探测的准确性，并进行了实验验证。2016年，Roussel等[14]使用参考站接收机，对高度角为2°～70°的全球定位系统（Globa Positioning System，GPS）和 GLONASS的 SNR数据进行了处理分析，并对从低仰角（2°～30°）和高仰角（30°～70°）SNR数据中提取出的观测量进行了融合，结果表明融合后观测量与土壤湿度间的相关性高于融合前观测量与土壤湿度间的相关性。2017年，Yang等[15]提出干涉信号的解析模型，并证明了从SNR数据中提取介电常数进而反演土壤湿度的可行性，结果表明北斗信号在土壤湿度反演方面效果明显。

国内方面，2015年，敖敏思等[16]结合仿真和实测土壤湿度数据、GPS观测值开展对比试验，发现SNR能有效跟踪土壤湿度变化，最大有效测量范围为45 m，并指出利用指数函数能较好地描述SNR多径延迟相位与土壤湿度之间的关系。同年，徐晓悦等[17]通过 Lomb-Scargle变换分析了SNR数据中的多径反射分量，分析结果显示，幅度与土壤湿度的相关性较强。汉牟田等[18]于2016年根据GNSS的SNR和干涉效应估计方法，推导出利用GNSS干涉信号幅度反演土壤湿度的半经验模型，并对此进行了仿真验证。2016年，严颂华等[19]探讨并通过实验验证了北斗B1波段信号干涉测量土壤水分的可行性。2017年，金双根等[20]对GNSS-R干涉测量的应用和植被覆盖问题进行了综合分析。2018年，段睿等[21]针对单天线模式下土壤湿度均方根误差较差等问题，提出SVRM辅助的土壤湿度反演方法，并通过实验证明该方法能有效提高土壤湿度均方根误差。

以上研究仅利用单一频点GNSS干涉信号反演土壤湿度，忽视了不同频点信号间的差异性及潜在的更丰富的观测信息。本文通过双频数据融合，综合双频观测信息，提升土壤湿度反演性能。

1 GNSS-IR土壤湿度反演原理

GNSS-IR是一种通过土壤反射的GNSS信号与直射信号的干涉效应实现土壤湿度反演的遥感技术[22]。GNSS接收天线在接收直射信号的同时，因周围环境（海面、土壤表面等）的反射也会接收到反射信号，当天线架设高度较低时，反射信号同直射信号相比具有相同的频率，因此可在天线处产生较为稳定的干涉效应。干涉场景如图1所示。

这一干涉效应体现在接收机记录的SNR中。SNR可以用直射和反射信号表示为[12]

式中：Ad、Am分别为直射、反射信号的幅度；ψ为两信号的相位差；θ为卫星高度角。当GNSS卫星高度角θ发生变化时，Ad、Am和ψ均会发生变化，造成SNR的波动。

图1 GNSS-R干涉测量原理图Fig.1 Schematic diagram of GNSS-R interferometry

根据式（1），将二阶多项式拟合法的结果作为直射分量，并将其从SNR序列中分离得到SNR的反射分量。通过频谱分析法得到SNR序列的频率，进而求得等效天线高度，选取值最大的等效天线高度值作为等效天线高度的估计值。

去除直射分量后，SNR的反射分量可表示为

式中：h为等效天线高度；λ为波长；φ为相位。

根据式（2）所述的形式，对 SNR的反射分量序列进行最小二乘拟合，得到相位观测量，建立其与实测土壤湿度之间的映射关系，即可用来反演土壤湿度。

2 基于熵值法的双频数据融合

根据香农信息论，信息熵可以量化地衡量信息的不确定性。在GNSS-IR土壤湿度反演中，不同频点SNR序列的相位观测量所含的土壤湿度信息不同，可以根据文献[23]所述的熵值法计算各频点SNR序列相位观测量的信息熵来确定其权值，加权作差后得到融合相位观测量。

具体地，对于频点k（k=1，2），第i颗卫星的第 j天相位观测量为其在整个观测过程中的相位观测量序列为X（k）i。

首先计算相位观测量序列 X（ik）与实测土壤湿度 SM的正负相关性 r（ik）：

并根据计算结果选择式（4）求得第i颗卫星第j天相位观测量的指标：

将第i颗卫星第j天相位观测量代入式（5），共N天的数据，计算得到第 i颗卫星第 j天相位观测量的比重 z（ijk）：

将第i颗卫星相位观测量的指标熵代入式（7），求得第i颗卫星相位观测量的熵值法冗余度：

根据第i颗卫星相位观测量熵值法冗余度和所有卫星冗余度和的比值，计算得到第i颗卫星相位观测量的指标熵的权重：

将第i颗卫星相位观测量的指标熵的权重代入式（9）计算指标评价得分，得到第 i颗卫星第 j天的相位观测量的指标评价得分：

按上述方法求得第i颗卫星第 j天各个频点相位观测量的指标评价得分后，分别与相应的原始相位观测量相乘并代入式（10）进行数据融合，得到融合后的相位观测量：

式中：Xij为融合后的相位观测量。

根据式（10）进行2个频点相位观测量的融合，将数据分为训练集与测试集，用训练集建立相位观测量与土壤湿度间的经验模型，并进行土壤湿度反演。

3 实验数据处理及结果分析

3.1 实验概况

本文使用2014-02-04—2014-03-21在法国图卢兹（Toulouse）市的 Lamasquere（43°29′14.45″N，1°13′44.11″E）采集的 GPS L1和 L2 SNR数据。实验场地四周空旷无遮蔽，整个实验期间为裸土，现场如图2所示。

实验使用Leica GR25接收机和AR10基准站天线，天线指向天顶，天线架设高度为1.70 m，采样频率为1 Hz。距离天线相位中心在地表投影约2 m的位置布置2枚ML3 Theta Probe土壤湿度传感器采集原位土壤湿度数据，精度为±1%，深度分别为2 cm和5 cm，采样间隔为10 min。

图2 GNSS-IR土壤湿度地基实验场地信息Fig.2 Information of GNSS-IR soil moisture ground-based experiment site

3.2 实验数据处理

数据处理流程如图3所示。

图3 数据处理过程Fig.3 Flow of data processing

具体过程如下：①获取实验数据。从采集的数据中提取SNR、时间、卫星高度角、方位角等信息，将 L1和 L2频点的 SNR定义为 SNR（1）、SNR（2）。②数据预处理。截取卫星高度角在2°～30°范围的SNR数据，同时删除SNR为零的数据。③数据分离。将数据分为训练集和测试集。④去除趋势项。对 SNR（1）、SNR（2）的卫星高度角正弦值进行二阶多项式拟合，并以此作为直射分量的近似值，去除直射分量，得到反射分量 SNR（m1）、SNR（m2）。⑤频谱分析。利用Lomb-Scargle方法对反射分量进行频谱分析，求得频率 f（1）和 f（2），进而求得等效天线高度 h（1）和 h（2）。⑥求相位观测量。根据式（2）对SNRm进行最小二乘拟合求得相位 观 测 量 P（1）和 P（2）。⑦ 数 据 融 合。根 据式（10）进行融合处理，得到融合后的相位观测量P。⑧建模并验证。选取融合前后的相位观测量分别与对应的土壤湿度观测值进行模型建立，将测试数据代入模型进行验证。

3.3 实验结果分析

根据3.2节处理过程，将L1和L2的SNR序列相位观测量进行融合，分别建立L1、L2及融合的相位观测量与同比土壤湿度间的一元线性回归模型，结果如图4～图7所示。

图 4～图 7分别为 PRN4、PRN7、PRN12及PRN24的建模结果，其中图 4（a）、图 5（a）、图 6（a）、图 7（a）为融合后的相位观测量与同比土壤湿度的建模结果，图 4（b）、图 5（b）、图 6（b）、图 7（b）为单频点（以 L1频点为例）的相位观测量与同比土壤湿度的建模结果。从图中可以直观地看出融合后的相位观测量与土壤湿度间的相关性更高。

图4 PRN4建模结果Fig.4 Modeling results of PRN4

图5 PRN7建模结果Fig.5 Modeling results of PRN7

图6 PRN12建模结果Fig.6 Modeling results of PRN12

图7 PRN24建模结果Fig.7 Modeling results of PRN24

模型建立完成后，将测试数据代入模型得到反演结果。同时，对单频SNR序列做同样处理得到单频反演结果。图 8给出了 PRN4、PRN7、PRN12、PRN24 4颗卫星的处理结果。

图8为各颗卫星反演结果对比，其中融合后PRN4的标准差为0.68%，比L1提高42.93%，比 L2提高 43.23%；均方根误差（RMSE）为0.49%，比L1降低81.58%，比L2降低77.52%；PRN7的标准差为0.62%，比L1提高43.85%，比L2提高 47.01%；RMSE为 0.98%，比 L1降低58.82%，比L2降低58.47%；PRN12的标准差为0.81%，比L1提高51.30%，比L2提高41.03%；RMSE为1.28%，比L1降低13.51%，比 L2降低46.67%；PRN24的标准差为0.92%，比 L1提高3.14%，比 L2降低 3.84%；RMSE为 0.95%，比L1降低64.15%。

对所有卫星同一时间段（13：00—15：00）的反演结果取平均，结果如图9所示。

图8 不同卫星反演结果对比Fig.8 Comparison of different satellite inversion results

图9（a）为融合前L1频点、L2频点和融合后反演结果以及FDR采集的同比数据随时间变化的对比结果，从对比中可以清晰显示融合后的反演结果与FDR的匹配度更高。图9（b）为反演结果线性对比，图中计算了融合前两频点的反演结果、融合后的反演结果与同比数据间的决定系数及 均 方 根 误 差 RMSEL1、RMSEL2、RMSE。通过计算得到如下结果：融合后的决定系数 R2为0.97，RMSE为 0.37%，R2比 L1频点结果提高73.2%，比 L2频点提高38.5%；RMSE比L1频点降低72.8%，比L2频点降低73.4%。

图9 反演结果平均对比分析Fig.9 Inversion results average comparison analysis

通过统计所有卫星的标准差，得到平均结果如下：L1和L2双频融合反演结果平均标准差为0.6%，精度比L1单频反演结果提高64.73%，比L2单频反演结果提高32.12%。

4 结 论

本文提出了一种基于GNSS-IR双频数据融合的土壤湿度反演方法，该方法以不同频点SNR序列相位观测量所对应的指标评价得分作为融合的权重标准，通过熵值法对2个频点SNR相位观测量进行融合处理。利用地基观测获取的实验数据对该方法的性能进行了验证，结果表明该方法可显著提高土壤湿度反演精度。从处理结果中可得到如下结论：

1）不同频点间所含的土壤信息存在差异。熵值法不仅能够克服主观确定权值的随机性和臆断性，而且作为多指标评价方法能够避免不同频点间信息的重复性，综合考虑信息间的差异性，两频点中质量好的数据能得以保留，从而使反演结果得到提高。

2）综合反演结果分析，熵值法在双频数据融合处理中能够有效提高反演结果的准确性，不仅反演精度有所提高，融合后数据的变化趋势也更加趋近于同比数据。最终结果也表明融合后的反演结果对比单一频点的反演结果有显著提高，且融合后的平均标准差为0.6%。

致谢感谢北京航空航天大学杨东凯教授在学术上给予的指导和帮助，感谢 F.Baup和 K.Boniface收集气象数据，以及 Roussel.N和 F.Frappart收集GNSS观测数据，感谢J.Darrozes为本文提供的数据。