土壤湿度是全球水循环过程的重要状态参数，在生态研究、农业灌溉及灾害预警等方面都起着非常重要的作用，因此研究土壤湿度的时空变化特征及大范围监测是十分必要的^[1]。基于全球卫星导航系统(Global Navigation Satellite System, GNSS)单天线技术的土壤湿度反演，是近年来全球导航卫星反射信号的一项新兴的应用，具有非接触、大面积、实时性、连续性的优点，弥补了传统的烘干称重法及土壤湿度传感器在空间和时间分辨率方面的不足。
GNSS单天线技术又称为全球卫星导航系统干涉测量法(GNSS-IR)，是利用单个GNSS接收机天线同时接收导航卫星直、反射信号，来反演信号反射面物理特征的一种新兴手段，也称为干涉图技术，最早由Rodriguez-Alvarez等^[2]于2009年提出，目前该方法已广泛应用在土壤湿度^[3]、海风海浪^[4-5]、积雪深度^[6]反演方面。Larson等^[7]提出利用传统的GPS接收机可以实现单天线的干涉测量，实测数据表明土壤湿度与卫星在低仰角时的信噪比(Signal to Noise Ratio，SNR)的振幅和相位有关。2016年，Roussel等^[8]利用测绘级接收机，对GPS和GLONASS高度角为2°~70°的SNR数据进行了处理，表明单测绘级接收机可以很好地探测表层土壤的湿度变化。国内方面，2016年，汉牟田等^[9]根据干涉效应和GNSS接收机信噪比估计方法，推导了利用GNSS干涉信号幅度进行土壤湿度反演的模型，并进行了仿真验证。2018年，吴继忠等^[10]针对GPS-IR获取土壤含水量的参数估计问题，提出了一种改进的反射信号参数估计方法，并研究了土壤含水量反演模型的建立过程。2019年，笔者团队提出GA-SVM辅助的土壤湿度反演方法，并通过实验证明该方法能有效提高土壤湿度的反演精度^[11]。
以上基于GNSS-IR技术的土壤湿度反演^[7-11]主要采用GPS系统单颗星的L1或L2波段。但是，不同频率电磁波对土壤湿度的敏感性不一样，仅采用单一频段单颗卫星数据忽略了不同频点卫星数据的差异性，而且数据来源也受到了较大的观测条件限制。因此，本文提出了基于GPS导航卫星系统的L1、L2和L5频段的多星数据融合方法进行土壤湿度反演。首先，对GPS多星的L1、L2和L5频段的测量数据划分训练集和测试集，利用训练集数据建立融合值与土壤湿度的反演模型；然后，通过测试集数据进行验证，同时与单星单频及其他反演方法进行比较。证明了该数据融合的反演方法在同等条件下效果更优，可以弥补单频单星因频率和轨道限制，造成的反演精度低及卫星信号差、多径效应强等问题。
1 单天线GNSS-IR土壤湿度反演原理 GNSS-IR模式采用一个右旋圆极化(Right Hand Circular Polarization)天线接收GNSS直射信号的同时，也接收地表的反射信号。直射信号为右旋圆极化信号，在低卫星高度角的情况下，反射信号极性仍以右旋圆极化为主^[12]，其场景如图 1所示。反射信号的相位和振幅特性的变化与地表反射面(土壤、海面、积雪等)的介质特性相关，通过测量直反射信号特性即可反演出地表的介质参数如土壤湿度。

图 1 干涉场景 Fig. 1 Scenario of interference

图选项

地基情况下由于平台固定，天线架设高度较低，直射跟反射信号频率近似相等，在接收天线处叠加会产生较为稳定的干涉信号，且在卫星低仰角情况下干涉现象更加明显。接收天线接收到的干涉信号SNR可以表示为^[13]

(1)

式中：A_d、A_m分别为直射、反射信号的幅度；ψ为直射信号和反射信号的相位差。
则有

(2)

式中：φ₁为直反射路径差导致的相位差；φ₀为干涉初始相位；H为接收机天线的等效高度，会随着反射面的介电特性的变化而变化；θ为卫星的高度角；λ为GNSS卫星信号的波长。
通过式(2)可以得到多径振荡的频率为

(3)

式中：分别为卫星高度角和天线等效高度相对于时间的变化率。
考虑到天线等效高度的变化率，在每天SNR观测的几小时周期内可以忽略不计，并进一步令x=sin θ，式(3)可简化为

(4)

式(4)说明, 振荡频率f与天线等效高度H呈线性关系。结合式(1)~式(4)，将直射信号剔除，只保留与反射面参数有关的多径信号，该多径信号SNR_m可以表示为

(5)

式中：A为振荡幅度; φ为初始相位。
通过最小二乘法对SNR_m进行拟合^[14]，可以得到多径信号SNR_m的振荡频率f、振荡幅度A和初始相位φ，在以往的文献中，通常从以上观测量中选取一种来建立与土壤湿度的经验模型进行反演，统称为Larson方法。本文选取实验与土壤湿度实测值相关性较好的SNR_m的振荡幅度作为反演观测量。
2 基于自适应融合算法的GPS多星三频数据融合 数据融合是20世纪80年代提出并得到快速发展的一种信息数据综合处理的方法^[15-17]，其充分利用多元数据的互补性来提高测量信息的质量。综合利用多星多频的测量数据，可以降低反演对象的不确定性，从而提高监测及反演的准确度和可靠性。
GPS目前可观测到共有32颗卫星，分布在6个不同的轨道平面。不同轨道、不同频率、不同功率的卫星，在同一反射面同一时间的测量数据存在较大差别，因此不能直接进行线性组合。本文研究了基于最小方差的自适应融合算法，不需要卫星测量数据的任何先验知识，对GPS多星的L1、L2和L5频段的训练集测量数据进行融合，建立融合值与土壤湿度的反演模型。
基于GPS多星三频数据融合的基本思想是：从接收机中提取整个观测周期内GPS系统L1、L2和L5频段所有卫星的SNR数据，利用单天线GNSS-IR土壤湿度反演原理，分别计算所有卫星每天的SNR_m振荡幅度观测量，并与土壤湿度传感器实测的土壤湿度真值进行相关性分析，选取符合同一观测时段(2 h内)，3个频段均相关性较好(与土壤湿度相关系数大于0.5)的多颗卫星数据利用自适应融合算法得到融合后的振荡幅度观测量。对于选取的第k颗卫星第i天的振荡幅度观测量可以表示为X_i^k(k=1, 2, …, m; i=1, 2, …, n)，m为符合条件的卫星数量，n为观测周期的天数，在实际计算时取第k颗卫星3个频段观测量的均值。其在整个观测周期n天内的振荡幅度序列可以表示为X^k=(X₁^k，X₂^k，…, X_n^k)，待估计的土壤湿度每天的测量真值为X_i(i=1, 2, …, n)。根据L1、L2和L5的观测值，可以计算出X_i^k在整个观测周期n天内每天对应的方差分别为(σ^k₁)², (σ^k₂)², …, (σ_n^k)²，每颗卫星的振荡幅度观测量X₁^k，X₂^k，…, X_n^k彼此相互独立，且为X_i的无偏估计。每颗卫星的加权因子定义为W₁^k，W₂^k，…，W_n^k，且满足。
则第i天的振荡幅度融合值可以表示为。自适应融合算法的总体思想就是在保证总方差最小的条件下，根据各卫星的振荡幅度观测量，通过自适应的方式找到最优的加权因子，使融合后的值最接近测量真值X_i。第i天的融合观测量的总方差为

(6)

由于各卫星的振荡幅度观测量相互独立，且都为X_i的无偏估计，考虑到：

(7)

则有总方差为

(8)

式中：(σ_i^k)²为第i天的振荡幅度观测量X_i^k的方差。
显然，总方差σ_i²为加权因子W_i^k的二次函数，存在最小值，且当σ_i²取最小值时，可以得到最优的加权因子W_i^k。利用拉格朗日乘数法构造函数：

(9)

对式(9)中的W_i^k求偏导数得到

(10)

使用拉格朗日乘数法求解极值，令，利用的条件，可以求得总方差最小时所对应的最优加权因子为

(11)

对应的总方差σ_i²的最小值为

(12)

由式(12)可以看出，自适应融合算法的特点是方差越小，测量值越接近真值，在算法中该测量值所占的比重越大，对应的加权因子越大；反之，方差大的所占比重越小，加权因子越小。从而可以得到最优的融合值为

(13)

式中：为GPS多星的振荡幅度最优融合观测量，利用该融合值与土壤湿度真值之间建立经验模型用于土壤湿度的反演。
数据处理流程如图 2所示。

图 2 数据处理流程 Fig. 2 Flowchart of data processing

图选项

3 实验与结果 3.1 实验数据采集 实验场地位于北京市通州区国家蔬菜工程技术研究中心的一处试验田(实验位置如图 3(a)所示)。2018-09-10—2018-11-09在这里进行了为期61 d的土壤湿度探测实验，实验场地东西长约200 m，南北长约50 m，范围如图 3(b)所示。

图 3 实验场地示意图 Fig. 3 Schematic diagram of experimental site

图选项

实验期间正逢玉米收割完毕，地表植被稀少，可以视为裸土，且周边较为空旷，有利于卫星信号接收，如图 3(c)所示。实验方面使用华测N72接收机和Antcom G5Ant-52AT1天线，天线架设高度为3 m，同时接收北斗的直反射信号，采样频率为1 Hz。由于卫星轨道大部分分布在南向，天线安装在试验田北侧边沿中心位置，如图 3(b)中的红点所示。在接收天线附近安装了3个土壤湿度传感器并连接到自动气象站，其中2个传感器水平埋设，深度分别为2 cm和4 cm，第3个传感器垂直放置，测量0~6 cm的土壤平均湿度，气象站采样间隔为1 min。实验期间有3次明显的降雨过程，有利于实验验证。
3.2 实验结果 根据第1节的原理对接收到的GPS SNR数据进行处理，得到L1、L2和L5的振荡幅度观测量，与3个土壤湿度传感器的测量均值进行相关性运算，根据相关性和时间段进一步筛选，由于GPS发射L5频段的卫星数量较少，最后得到符合融合条件的卫星为PRN1、PRN6和PRN8。将数据集按2:1比例划分成训练集和测试集，利用Larson方法可以分别得到3颗星L1、L2和L5频段的单星单频训练集反演模型和测试集反演结果(见图 4)。从图 4中可以看出，振荡幅度观测量与土壤湿度具有明显的线性关系，反演结果反映了土壤湿度变化的趋势。但由于实验场地玉米收割后，土壤粗糙度较高、地势不平，影响了单星的反射信号接收，因此单频单星的反演结果差异较大且反演结果较差，以反演相关性最好的PRN8的L5频段为例(见图 4(c)、(f))，相关系数R仅为0.646 3，拟合度较低，反演精度较差。

图 4 单星反演模型及结果 Fig. 4 Inversion model and results of single satellite

图选项

为了充分利用单星单频数据的差异性和互补性，根据第2节的处理过程，对PRN1、PRN6和PRN8的3个频段的训练集振荡幅度观测量进行基于最小方差的自适应融合，得到多星三频的融合观测量与土壤湿度的训练集反演模型及测试集反演结果(见图 5(a)、(b))。由图 5(b)中可以看出，本文自适应融合算法得到的测试集反演值与土壤湿度实测值具有较好的相关性，相关系数R达到了0.805 9，并且从反演结果上可以看出反演值与实测值拟合度较高。在图 5(c)、(d)中，还进行了自适应融合算法与传统的利用单星单频的Larson方法，以及对所有频段观测量取均值的均值融合法进行了比较，结果表明，本文提出的自适应融合算法反演值更接近于实测值，相比于Larson方法(PRN8，L5，R=0.646 3), 相关系数R提高了24.69%，相比于均值融合法(R=0.635 7)，相关系数提高了26.77%，说明本文提出的自适应融合算法反演精度高，更能克服土壤粗糙度及地形起伏的影响。

图 5 融合反演模型及结果对比 Fig. 5 Fusion inversion model and results comparison

图选项

为了进一步验证自适应融合反演模型的优越性，本文还通过平均绝对误差(MAE)、最大相对误差(MRE)、均方根误差(RMSE)、平均绝对百分比误差(MAPE) 4项评价指标对上述3种反演模型的反演结果进行评价对比, 对比结果如表 1所示。表 1中, x_i为实测值; x_i为反演值，i=1, 2, …, n。从表 1中可以看出, 自适应融合算法各项误差均较小，均方根误差RMSE为2.075%，相比于Larson方法下降了22.28%，相比于均值融合法下降了23.26%。
表 1 土壤湿度反演结果评价比较 Table 1 Evaluation and comparison of soil moisture inversion results %

评价指标	计算公式	自适应融合算法	Larson方法	均值融合法
MAE		1.46	1.919	2.088
MRE		19.427	24.995	22.938
RMSE		2.075	2.67	2.704
MAPE		4.836	6.372	6.887

表选项






4 结论 本文在基于GNSS-IR的Larson方法的基础上提出了GPS多星三频数据的自适应融合土壤湿度反演方法，给出了数据处理的一般流程，并开展了实验对该方法进行了验证测试。实验结果表明：
1) 在低高度角条件下(2°~30°)，利用该方法得到的测试集反演值与土壤湿度实测值相关度较好，相关系数达到了0.805 9，可以较好地实现对固定区域土壤湿度的连续监测。
2) 与传统的Larson方法相比，本文提出的自适应融合算法相关系数提高了24.69%，比均值融合法提高了26.77%，均方根误差RMSE较后2种方法均有比较明显下降，进一步验证了本文模型的有效性。
综上所述，基于GPS多星三频数据融合的GNSS-IR土壤湿度反演方法充分利用了不同频率、不同轨道卫星测量数据的差异性，利用多星数据的互补性来提高融合数据的质量，从而提高了GNSS-IR技术在土壤湿度监测及反演方面的准确性和可靠性。
下一步随着北斗卫星系统的完善，基于北斗系统的土壤湿度反演是后续研究的主要内容之一。
致谢 感谢国家蔬菜工程技术研究中心提供的实验场地支持，同时感谢韦孝海先生在实验中提供的帮助。

参考文献

[1]	袁苇, 李宗谦, 刘宁, 等. 不同微波遥感模态和不同数据组合的等湿度区域分布与土壤湿度的反演[J]. 中国工程科学, 2004, 6(6): 50-56. YUAN W, LI Z Q, LIU N, et al. Analysis of data sets with different microwave remote sensing model in soil moisture retrieval[J]. Engineering Science, 2004, 6(6): 50-56. DOI:10.3969/j.issn.1009-1742.2004.06.008 (in Chinese)
[2]	RODRIGUEZ-ALVAREZ N, BOSCH-LLUIS X, CAMPS A, et al. Soil moisture retrieval using GNSS-R techniques:Experimental results over a bare soil field[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(11): 3616-3624. DOI:10.1109/TGRS.2009.2030672
[3]	LARSON K M, SMALL E E, GUTMANN E, et al. Using GPS multipath to measure soil moisture fluctuations:Initial results[J]. GPS Solutions, 2008, 12(3): 173-177. DOI:10.1007/s10291-007-0076-6
[4]	SANTAMARíA-GóMEZ A, WATSON C, GRAVELLE M, et al. Levelling co-located GNSS and tide gauge stations using GNSS reflectometry[J]. Journal of Geodesy, 2015, 89(3): 241-258.
[5]	CLARIZIA M P, RUF C S, JALES P, et al. Spaceborne GNSS-R minimum variance wind speed estimator[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(11): 6829-6843. DOI:10.1109/TGRS.2014.2303831
[6]	JIN S G, NAJIBI N. Sensing snow height and surface temperature variations in greenland from GPS reflected signals[J]. Advances in Space Research, 2014, 53(11): 1623-1633. DOI:10.1016/j.asr.2014.03.005
[7]	LARSON K M, BRAUN J J, SMALL E E, et al. GPS multipath and its relation to near-surface soil moisture content[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2010, 3(1): 91-99. DOI:10.1109/JSTARS.2009.2033612
[8]	ROUSSEL N, FRAPPART F, RAMILLIEN G, et al. Detection of soil moisture variations using GPS and GLONASS SNR data for elevation angles ranging from 2° to 70°[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(10): 4781-4794. DOI:10.1109/JSTARS.2016.2537847
[9]	汉牟田, 张波, 杨东凯, 等. 利用GNSS干涉信号振荡幅度反演土壤湿度[J]. 测绘学报, 2016, 45(11): 1293-1300. HAN M T, ZHANG B, YANG D K, et al. Soil moisture retrieval utilizing GNSS interference signal amplitude[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(11): 1293-1300. DOI:10.11947/j.AGCS.2016.20160145 (in Chinese)
[10]	吴继忠, 王天, 吴玮, 等. 利用GPS-IR监测土壤含水量的反演模型[J]. 武汉大学学报信息科学版, 2018, 43(6): 887-892. WU J Z, WANG T, WU W, et al. Retrieval model for soil moisture content using GPS-interferometric reflectometry[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 887-892. (in Chinese)
[11]	孙波, 梁勇, 汉牟田, 等. 基于GA-SVM的GNSS-IR土壤湿度反演方法[J]. 北京航空航天大学学报, 2019, 45(3): 486-492. SUN B, LIANG Y, HAN M T, et al. GNSS-IR soil moisture inversion method based on GA-SVM[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(3): 486-492. (in Chinese)
[12]	ARROYO A A, CAMPS A, AGUASCA A, et al. Dual-polarization GNSS-R interference pattern technique for soil moisture mapping[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(5): 1533-1544. DOI:10.1109/JSTARS.2014.2320792
[13]	BILICH A, LARSON K M. Correction published 29 March 2008:Mapping the GPS multipath environment using the signal-to-noise ratio(SNR)[J]. Radio Science, 2016, 42(6): 1-16.
[14]	CHEW C C, SMALL E E, LARSON K M, et al. An algorithm for soil moisture estimation using GPS-interferometric reflectometry for bare and vegetated soil[J]. GPS Solutions, 2016, 20(3): 1-13.
[15]	CHAIR Z, VARSHNEY P K. Optimal data fusion in multiple sensor detection systems[J]. IEEE Transactions on Aerospace & Electronic Systems, 2007, 22(1): 98-101.
[16]	CULP M, MICHAILIDIS G. A co-training algorithm for multi-view data with applications in data fusion[J]. Journal of Chemometrics, 2010, 23(6): 294-303.
[17]	GHANBARI M, YAZDANPANAH M J. Delay compensation of tilt sensors based on mEMS accelerometer using data fusion technique[J]. IEEE Sensors Journal, 2015, 15(3): 1959-1966.