顾及多径环境差异和粗差的GNSS-MR土壤湿度反演

顾及多径环境差异和粗差的GNSS-MR土壤湿度反演 GNSS-MR soil moisture retrieval considering the multipath environments differences and gross error first-author 李婷，1993年生，女，硕士研究生，研究方向为GNSS-MR土壤湿度反演。E-mail： 1506087124@qq.com 李婷，1993年生，女，硕士研究生，研究方向为GNSS-MR土壤湿度反演。E-mail： 1506087124@qq.com 1506087124@qq.com corresp 张显云，1974年生，男，副教授，研究方向为GNSS理论及应用研究。E-mail： mec.xyzhang@gzul.edu.cn 张显云，1974年生，男，副教授，研究方向为GNSS理论及应用研究。E-mail： mec.xyzhang@gzul.edu.cn mec.xyzhang@gzul.edu.cn 贵州大学 矿业学院, 贵阳 550025 贵州大学 矿业学院, 贵阳 550025 College of Mining, Guizhou University, Guiyang 550025, China College of Mining, Guizhou University, Guiyang 550025, China 针对基于单系统单卫星GNSS-MR（GNSS Multipath Reflectometry）土壤湿度反演的可靠性不高、实际可操作性不强和最小二乘估计不具鲁棒性的缺点，为获取更优的延迟相位估值，并改善GNSS-MR土壤湿度反演的可靠性和实际可操作性，同时简化繁杂的选星过程，提出了一种基于抗差估计的多系统多卫星组合GNSS-MR土壤湿度反演算法。该算法首先顾及多径环境的差异性、多径误差的周期特性等进行信噪比SNR（Signal to Noise Ratio）观测值的筛选，然后采用基于IGGIII（Weight Function III Developed by Institute of Geodesy and Geophysics）权函数的抗差估计解求延迟相位，进而获得表征土壤湿度变化趋势的延迟相位组合。实验结果表明，相较于未采用抗差估计的多系统多卫星组合（方案1）和单卫星组合（方案3），得益于抗差估计良好的鲁棒性，基于抗差估计的多系统多卫星组合（方案2）和单卫星组合（方案4）获得了较高的建模精度，所得延迟相位与实测土壤湿度间的相关系数分别为0.97和0.95、土壤湿度拟合残差的均方根误差分别为0.010和0.012；同时，方案2和方案4还取得了较高的土壤湿度预报精度，土壤湿度预测值与土壤湿度实测值间的相关系数分别为0.92和0.91、土壤湿度预报残差的均方根误差分别为0.016和0.023；此外，相比于方案4，方案2在采用抗差估计解求延迟相位的基础上，采用多系统多卫星组合进一步提升了延迟相位的估值精度，从而不仅避免了复杂的选星过程，而且还获得了更好的建模效果和更高的土壤湿度预报精度。 In order to obtain better phase delay estimation, improve the reliability and the practical operability of GNSS-MR(GNSS Multipath Reflectometry) soil moisture inversion, and also to simplify the complex process of the satellite selection, a multi-system multi-satellite GNSS-MR soil moisture inversion algorithm based on the robust estimation was proposed in view of the poor reliability and operability of the single system single satellite GNSS-MR soil moisture inversion and the least squares estimation of no robustness. In this algorithm, the spatial difference of multipath environment and the periodic characteristics of multipath were taken into account to screen the SNR (Signal to Noise Ratio) observations. Then, the phase delay combination representing the change trend of soil moisture was obtained by using the robust estimation based on IGGIII (Weight Function III Developed by Institute of Geodesy and Geophysics) weight function. Compared with multi-system multi- satellite combination (scheme 1) and the single-satellite combination (scheme 3), the experimental results showed that the multi-system multi-satellite combination (scheme 2) and the single-satellite combination (scheme 4) based on the robust estimation achieved higher modeling accuracy, which were benefited from the positive performance of the robust estimation. The correlation coefficients between the estimated phase delays and the measured soil moisture were 0.97 and 0.95, respectively, and the root mean square error of the soil moisture fitting residual were 0.010 and 0.012, respectively. At the same time, scheme 2 and scheme 4 also achieved higher soil moisture prediction accuracy, with the correlation coefficient between the predicted soil moisture and the measured soil moisture being 0.92 and 0.91, respectively, and the root mean square error of the soil moisture forecast residuals being 0.016 and 0.023, respectively. In addition, compared with scheme 4, scheme 2 not only adopted the robust estimation, but also adopted the multi-system multi-satellite combination, which contributed to better modeling effect and higher modeling accuracy. Moreover, because it could avoid the complex process of the satellite selection, scheme 2 owned better performance in GNSS-MR soil moisture inversion. 多系统多卫星GNSS-MR 土壤湿度反演 信噪比 延迟相位 抗差估计 multi-system multi-satellite GNSS-MR soil moisture inversion signal to noise ratio delayed phase robust estimation 1　引言 1 土壤湿度亦称土壤含水率，是表征土壤干湿程度的物理量，是生态系统的重要影响因子；土壤湿度对天气预报和气候研究有着重要价值；准实时的土壤湿度动态监测对边坡稳定性预判、精准农业等均有着积极意义（ Robock等，2000 Robock等，2000 ； Schaufler等，2010 Schaufler等，2010 ； Seneviratne等，2010 Seneviratne等，2010 ； Seuffert等，2003 Seuffert等，2003 ； Brocca等，2010 Brocca等，2010 ； Wanders等，2014 Wanders等，2014 ； Vey等，2016 Vey等，2016 ； Roussel等，2015 Roussel等，2015 ； Ray和Jacobs，2007 Ray和Jacobs，2007 ）。然而，传统的土壤湿度测量手段一般是定时、定点测量，具有时空分辨率低、土壤湿度样本分布不均匀、不易于实现大范围土壤湿度监测等缺点。遥感手段虽能实现全球尺度的土壤湿度监测，但由于目前遥感卫星复现性较差，导致土壤湿度时间分辨率较低，即使采用ESA（European Space Agency）的SMOS（Soil Moisture and Ocean Salinity）卫星，其时间分辨率也仅为3 d（ Kerr 等，2010 Kerr 等，2010 ）。全球导航卫星系统GNSS除可广泛应用于定位、导航和授时外，还能无偿提供全天候的、全球覆盖的、准实时的L波段（1—2 GHz）微波信号用于遥感探测（ 金双根 等，2017 金双根 等，2017 ； 万玮 等，2016 万玮 等，2016 ）。GNSS-MR是一种直接利用一定截止高度角范围内的多路径信噪比SNRMP（Signal to Noise Ratio Multipath）（亦称为信噪比残差）反射分量之幅度、频率和相位的变化，进行近地表物理参数反演的一种全新的遥感技术，已被证实可用于潮位变化（ 张双成 等，2016 张双成 等，2016 ）、雪深（ Yu 等，2015 Yu 等，2015 ）等地表环境的实时监测，且可克服传统测量手段对被监测对象的破坏、仪器类型间数据难以同化、时空分辨率受限或高成本等缺点。 土壤湿度反演方面，国外学者 Larson等（2008） Larson等（2008） 率先采用GPS-MR技术进行实验，并得出了SNR多路径反射分量的幅度与土壤湿度间存在相关性的结论，继而揭开了采用GNSS-MR技术进行土壤湿度反演研究的序幕。然而，迄今为止国内外 GNSS-MR土壤湿度反演的相关研究采用的几乎都是单系统、单卫星的信噪比观测值（ Larson等，2010 Larson等，2010 ； Zavorotny等，2010 Zavorotny等，2010 ； Chew等，2015 Chew等，2015 ； Wan等，2015 Wan等，2015 ； 敖敏思 等，2012 敖敏思 等，2012 ， 2015 2015 ； Yang等，2017 Yang等，2017 ）。由于：（1）GNSS同一卫星单天中处于低高度角状态的持续时间较短；（2）土壤含水量及土壤构成往往存在空间差异性，而短时间内单系统、单卫星在相近方位的复现性较差；（3）卫星信号的实际反射面并非GNSS-MR假设的理想水平反射面；加之复杂的菲涅尔效应等的综合影响，导致基于单系统、单卫星GNSS-MR的土壤湿度反演可靠性及稳定性较差，甚至可能无法实现短期内土壤湿度变化的高精度、准动态监测。 此外，当前基于GNSS-MR的土壤湿度反演研究几乎都是通过构建低阶多项式去除趋势项，从而获得表征地表反射物物理信息的信噪比残差序列，并据此反演土壤湿度的变化趋势。其效果主要依赖于SNR的观测质量和SNR中直接信号的成功移除（ Yu等，2015 Yu等，2015 ）。然而，SNR观测值难免会受到异常噪声的污染，从而导致提取的多路径信噪比往往“不纯”，继而影响土壤湿度反演的准确度，而当前GNSS-MR土壤湿度反演中对SNR MP 提纯算法的研究基本未涉及。再者，要从众多的SNR观测值中选择一组最优的SNR组合，实际操作难度较大，继而将影响GNSS-MR土壤湿度反演的效率。 为此，本文在分析GNSS-MR相关理论基础上，为充分挖掘多模GNSS组合的优势，利用GPS/BDS/GLONASS/Galileo（为区别卫星系统，后续图表中分别缩写为G/C/R/E）组合进行土壤湿度的反演，以弥补基于单系统、单卫星的GNSS-MR土壤湿度反演可靠性不高、稳定性较差的缺陷；同时，引入基于IGGIII权函数的抗差估计方法解求延迟相位，以克服采用多项式去趋势项后所得多路径信噪比不纯，导致延迟相位估值相位精度不高的问题；最后，结合实测数据分析了文中算法的性能。 2　原理与方法 1 2.1　GNSS-MR土壤湿度反演原理 2 信噪比观测值SNR是常规测量型GNSS接收机的附带观测值，是表征接收机天线接收信号强弱的物理量，该值为信号功率与噪声功率的比值，其亦是反映载波信号质量的主要指标，受天线增益参数、卫星与接收机间的几何距离、多路径效应等因素的共同影响。如 图1 图1 所示，当GNSS测站位于近地表土壤时，接收机接收到的信号实际为直射信号和反射信号构成的合成信号。 图1 图1 中，E为卫星高度角。由于低卫星高度角的SNR观测值往往能够较好地承载着测站周围多路径反射物丰富的物理信息，基于GNSS-MR的相关研究一般以卫星高度角5°—30°的SNR观测值作为重点关注对象。研究表明，SNR具有如下两个特点：（1）SNR与卫星高度角呈正相关关系（ 尹子明 等，2016 尹子明 等，2016 ）；（2）多路径效应增大时，SNR明显降低（ 吴丹 等，2015 吴丹 等，2015 ）。 大多数GNSS接收机记录的SNR观测值往往为SNR的振幅，其计算式为 <math id="M1"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">B</mi><mo>-</mo><mi mathvariant="normal">H</mi><mi mathvariant="normal">z</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">20</mn><mo>⋅</mo><mi mathvariant="normal">l</mi><mi mathvariant="normal">g</mi><mfenced separators="|"><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">v</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi><mo>/</mo><mi mathvariant="normal">v</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194520&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194513&type= 46.14333344 5.24933338 、单位为对数单位dB-Hz。由于以dB-Hz为单位的SNR观测值对卫星高度角的变化响应不积极，振幅幅值波动起伏变化较小，故为便于分析和阅读，一般将SNR观测值按 式（1） 式（1） 进行线性化转换。 <math display="block" id="M2"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">v</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi><mo>/</mo><mi mathvariant="normal">v</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">s</mi></mrow></msub><mo>=</mo><msup><mrow><mn mathvariant="normal">10</mn></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">B</mi><mo>-</mo><mi mathvariant="normal">H</mi><mi mathvariant="normal">z</mi><mo>/</mo><mn mathvariant="normal">20</mn></mrow></msub></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194530&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194524&type= 34.20533371 5.07999992 在仅存单一反射的假设基础上，SNR与信号幅度存在如下关系（ Strode和Groves，2016 Strode和Groves，2016 ）： <math display="block" id="M3"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msup><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>=</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">C</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>=</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">m</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><mn mathvariant="normal">2</mn><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">m</mi></mrow></msub><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mi>φ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194548&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194541&type= 54.10200119 4.74133301 式中，SNR表示接收机记录的合成信噪比； <math id="M4"><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">C</mi></mrow><mrow/></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194561&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194555&type= 3.47133350 4.06400013 为合成信号的振幅； <math id="M5"><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow><mrow/></msubsup><mi mathvariant="normal">、</mi><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">m</mi></mrow><mrow/></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194574&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194568&type= 8.97466660 4.31799984 分别为直、反射信号的振幅； <math id="M6"><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mi>φ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194587&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194581&type= 6.94266701 3.30200005 表示延迟相位的余弦值。鉴于直反射信号数值上的差异（ <math id="M7"><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">m</mi></mrow></msub><mo>≪</mo><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">d</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194598&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194592&type= 12.36133385 3.72533321 ），加之SNR随卫星高度角整体呈二次曲线变化趋势（ Roussel等，2016 Roussel等，2016 ），构建低阶多项式便可分离直射分量，从而获得近地表土壤湿度等多路径因子引起的多路径信噪比SNR MP 。 图2 图2 为实验站点GPS卫星PRN G28（DOY：128（2017年），UTC：6.92—13.64 h）的合成信噪比（SNR）、合成信噪比拟合值（fitted SNR）及去趋势项后的信噪比残差序列（SNR MP ），横轴为信噪比采集时间。 图3 图3 为G28去除趋势项后的SNR MP 的重采样图，横轴为卫星高度角的正弦值。SNR MP 主要表征多路径噪声水平的高低。可看出，SNR MP 具有明显的周期特性。SNR MP 可表示为 <math display="block" id="M8"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mo>=</mo><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mfenced separators="|"><mrow><mfrac><mrow><mn mathvariant="normal">4</mn><mi mathvariant="normal">π</mi><mi>h</mi></mrow><mrow><mi>λ</mi></mrow></mfrac><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>E</mi><mo>+</mo><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194641&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194636&type= 52.06999969 8.80533314 式中， <math id="M9"><mi>λ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195372&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195370&type= 2.28600001 2.62466669 为载波波长、 <math id="M10"><mi>E</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195377&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195376&type= 2.28600001 2.62466669 为卫星高度角、 <math id="M11"><mi>h</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195375&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195373&type= 1.86266661 2.62466669 为垂直反射距离、 <math id="M12"><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194691&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194685&type= 4.91066647 3.72533321 和 <math id="M13"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194704&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194697&type= 4.99533367 4.40266657 分别为多路径反射分量的振幅和延迟相位。若记 <math id="M14"><mi>t</mi><mo>=</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>E</mi><mo>,</mo><mi>f</mi><mo>=</mo><mfrac><mrow><mn mathvariant="normal">2</mn><mi>h</mi></mrow><mrow><mi>λ</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194716&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194711&type= 25.06133461 7.70466709 ，则（3）式可进一步简化为 <math display="block" id="M15"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mo>=</mo><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mfenced separators="|"><mrow><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi>f</mi><mi>t</mi><mo>+</mo><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194729&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194723&type= 43.68799973 5.24933338 由于 <math id="M16"><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>E</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194743&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194734&type= 7.19666624 2.70933342 为非等间距采样，导致SNR MP 无法保证整周期截断，可首先通过Lomb-Scargle谱分析诊断出多路径反射信号的频率 <math id="M17"><mi>f</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194750&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194747&type= 1.10066664 3.47133350 ，继而借助最小二乘解求其幅度和延迟相位（ 徐斌 等，2011 徐斌 等，2011 ）。延迟相位和土壤湿度具有较强的相关性，依据延迟相位便可反演土壤湿度。 2.2　基于Multi-GNSS的土壤湿度反演 2 据上述分析知，GNSS-MR之所以能够实现对土壤湿度的反演，主要得益于低高度角卫星的信噪比较好地承载了近地表反射物丰富的物理信息。然而，GNSS卫星的运动一般具有周期性，卫星高度角在一天中的变化较大且卫星低高度角状态持续时间较短，导致单颗卫星的信噪比能较好地承载表征地表反射物物理特性的有效时间非常有限。对实验地点GPS卫星（PRN1-PRN11）5°—25°高度角区间的持续时间及方位角变化情况分析知，单颗卫星5°—25°高度角持续时间仅为60 min左右。若顾及土壤湿度采样时间、反射介质的单一性及土壤平整性等因素的影响，在限制卫星方位角情况下，大部分时段均无法挑选出适宜的GPS信噪比观测值，且每个时段所选卫星往往不同，从而导致基于单系统、单卫星GNSS-MR的土壤湿度反演可靠性不高、稳定性较差。 事实上，GPS、GLONASS、Galileo均为中高轨卫星MEO（Medium Earth Orbit），如 图4 图4 中（a）—（d）所示，低高度角区间（红色矩形框部分）MEO卫星信噪比具有明显的周期特性，其均可作为GNSS-MR土壤湿度反演的数据源。值得一提的是，BDS由地球静止轨道卫星GEO（Geosynchronous Earth Orbit）、倾斜轨道同步卫星IGSO（Inclined Geosynchronous Satellite Orbit）和MEO 3种不同轨道的卫星组成。其中，IGSO卫星为区域服务卫星，运动轨迹投影呈南北“8”字形，其低高度角状态时持续时间较长。由 图4（e） 图4（e） 中可看出，BDS IGSO卫星低高度角状态的SNR观测值具有与MEO卫星信噪比类似的周期特性，同样可作为GNSS-MR土壤湿度反演的数据源；而BDS GEO卫星由于其轨道（地球静止轨道）的特殊性，卫星高度角较高且变化幅度较小，由 图4（f） 图4（f） 中可看出，BDS GEO卫星的SNR观测值不能较好地承载近地表反射物的多路径信息，周期特性较差。因此，基于单台BDS接收机所采集的GEO卫星信噪比，尚无法实现对土壤湿度的反演。当前，GNSS已呈现出多模多频的发展格局，同一时间同一地点可视卫星数目成倍增加，使得采用多系统组合进行土壤湿度的反演成为可能。 1,2;3,4;5,6;7; 2.3　延迟相位的解求 2 为解求延迟相位，需对观测方程 式（4） 式（4） 进行线性化处理，得如下线性化方程。 <math display="block" id="M18"><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mo>+</mo><mi>v</mi><mo>=</mo><mi>F</mi><mo stretchy="false">(</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo>,</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo stretchy="false">)</mo><mo>+</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mfenced separators="|"><mrow><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi>f</mi><mi>t</mi><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></mrow></mfenced><mi mathvariant="normal">d</mi><mi>A</mi><mo>-</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mfenced separators="|"><mrow><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi>f</mi><mi>t</mi><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></mrow></mfenced><mi mathvariant="normal">d</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194920&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194914&type= 109.89733887 4.57200003 式（5） 式（5） 对应误差方程可表示为 <math display="block" id="M19"><mi mathvariant="bold-italic">V</mi><mo>=</mo><mi mathvariant="bold-italic">B</mi><mover accent="true"><mi mathvariant="bold-italic">X</mi><mo>̂</mo></mover><mo>-</mo><mi mathvariant="bold-italic">L</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194932&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194927&type= 18.37266731 3.47133350 式中， v 为单个SNR MP 观测值的残差； <math id="M20"><mi mathvariant="bold-italic">V</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194943&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194937&type= 2.28600001 2.62466669 为SNR MP 观测值的残差序列构成的列向量； <math id="M21"><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194957&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194950&type= 4.91066647 4.06400013 、 <math id="M22"><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195070&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195066&type= 4.99533367 4.65666676 分别为SNR MP 的振幅和延迟相位的初值，可在采用Lomb-Scargle谱分析诊断出各SNR MP 频率 <math id="M23"><mi>f</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194982&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194976&type= 1.10066664 3.47133350 的基础上，可由 式（4） 式（4） 采用进行周期拟合进行确定； <math id="M24"><mi>F</mi><mo stretchy="false">(</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo>,</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194992&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20194987&type= 17.44133186 4.82600021 为 式（4） 式（4） 分别对SNR MP 振幅 <math id="M25"><msub><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195007&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195000&type= 4.91066647 3.72533321 、延迟相位 <math id="M26"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195021&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195013&type= 4.99533367 4.40266657 求一阶导数对应的初值；参数估值 <math id="M27"><mover accent="true"><mi mathvariant="bold-italic">X</mi><mo>̂</mo></mover><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><mi mathvariant="normal">d</mi><mi>A</mi></mtd></mtr><mtr><mtd><mi mathvariant="normal">d</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195032&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195027&type= 16.84866714 8.89000034 ，其中 <math id="M28"><mi mathvariant="normal">d</mi><mi>A</mi><mi mathvariant="normal">、</mi><mi mathvariant="normal">d</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195046&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195040&type= 12.61533356 4.40266657 分别为SNR MP 振幅初值 <math id="M29"><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195061&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195051&type= 4.91066647 4.06400013 和延迟相位初值 <math id="M30"><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195070&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195066&type= 4.99533367 4.65666676 的改正数； <math id="M31"><mi>L</mi><mo>=</mo><mi mathvariant="normal">S</mi><mi mathvariant="normal">N</mi><msub><mrow><mi mathvariant="normal">R</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mo>-</mo><mi>F</mi><mo stretchy="false">(</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo>,</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195081&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195075&type= 38.01533508 4.82600021 为常数项；其余符号含义同上。此外，误差方程系数阵 B 可将 <math id="M32"><mi>f</mi><mi mathvariant="normal">、</mi><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mi mathvariant="normal">、</mi><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195098&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195089&type= 14.81666565 4.82600021 代入 式（5） 式（5） 得以确定。 <math display="block" id="M33"><mi mathvariant="bold-italic">B</mi><mo>=</mo><mfenced separators="|"><mrow><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi>f</mi><mi>t</mi><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo stretchy="false">)</mo><mo>-</mo><msubsup><mrow><mi>A</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mo stretchy="false">(</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi>f</mi><mi>t</mi><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo stretchy="false">)</mo></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195112&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195102&type= 66.37866211 4.57200003 继而依据最小二乘估计得 <math display="block" id="M34"><mover accent="true"><mi mathvariant="bold-italic">X</mi><mo>̂</mo></mover><mo>=</mo><msup><mrow><mfenced separators="|"><mrow><msup><mrow><mi mathvariant="bold-italic">B</mi></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msup><mi mathvariant="bold-italic">P</mi><mi mathvariant="bold-italic">B</mi></mrow></mfenced></mrow><mrow><mo>-</mo><mn mathvariant="normal">1</mn></mrow></msup><mfenced separators="|"><mrow><msup><mrow><mi mathvariant="bold-italic">B</mi></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msup><mi mathvariant="bold-italic">P</mi><mi mathvariant="bold-italic">L</mi></mrow></mfenced><mtext> </mtext></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195129&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195122&type= 33.44333267 4.91066647 <math display="block" id="M35"><mi mathvariant="bold-italic">P</mi><mo>=</mo><mi mathvariant="bold-italic">E</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195141&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195136&type= 9.39799976 2.62466669 现有的GNSS-MR土壤湿度反演研究普遍忽略了观测值本身的质量差异，默认各SNR MP 观测值对参数估值的贡献相同。然而，实际信噪比观测值难免会受到异常噪声的污染，导致由多项式拟合去趋势项后提取的SNR MP 往往“不纯”。如 图5 图5 所示，E24和R19红圈部分的SNR MP 明显受到了强噪声的污染。 1,2;3; 为此，本文在参数平差过程中引入抗差M估计，参数抗差解为 <math display="block" id="M36"><mover accent="true"><mi mathvariant="bold-italic">X</mi><mo>̂</mo></mover><mo>=</mo><msup><mrow><mfenced separators="|"><mrow><msup><mrow><mi mathvariant="bold-italic">B</mi></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msup><msub><mrow><mover accent="true"><mi mathvariant="bold-italic">P</mi><mo>¯</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub><mi mathvariant="bold-italic">B</mi></mrow></mfenced></mrow><mrow><mo>-</mo><mn mathvariant="normal">1</mn></mrow></msup><mfenced separators="|"><mrow><msup><mrow><mi mathvariant="bold-italic">B</mi></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msup><msub><mrow><mover accent="true"><mi mathvariant="bold-italic">P</mi><mo>¯</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub><mi mathvariant="bold-italic">L</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195220&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195213&type= 35.81399918 5.84200001 式中， P 为等价权；权函数采用IGGIII权函数（ 杨元喜，1993 杨元喜，1993 ），即： <math display="block" id="M37"><msub><mrow><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>P</mi></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>P</mi></mrow><mrow><mi>i</mi><mi>j</mi></mrow></msub><mfrac><mrow><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow><mrow><mfenced open="|" close="|" separators="|"><mrow><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></mrow></mfenced></mrow></mfrac><msup><mrow><mfenced separators="|"><mrow><mfrac><mrow><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>-</mo><mfenced open="|" close="|" separators="|"><mrow><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></mrow></mfenced></mrow><mrow><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>-</mo><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfrac></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn><mtext>  </mtext></mrow></msup></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195229&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195225&type= 47.92133331 25.56933212 <math display="block" id="M38"><mtable columnalign="left"><mtr><mtd><mfenced open="|" close="|" separators="|"><mrow><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></mrow></mfenced><mo><</mo><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>≤</mo><mfenced open="|" close="|" separators="|"><mrow><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></mrow></mfenced><mo><</mo><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mfenced open="|" close="|" separators="|"><mrow><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></mrow></mfenced><mo>≥</mo><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195240&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195235&type= 24.97666740 22.01333237 式中，下标 i 为SNR MP 的编号； P ij 为 i 个观测值的权； <math id="M39"><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195253&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195250&type= 2.11666679 3.72533321 为观测值残差； <math id="M40"><mi>σ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195262&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195257&type= 2.45533323 2.62466669 为观测值中误差； <math id="M41"><mrow><mrow><msub><mrow><mi>v</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mo>/</mo><mrow><mi>σ</mi></mrow></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195268&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195264&type= 7.45066643 4.23333359 为标准化残差； <math id="M42"><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195275&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195272&type= 3.47133350 3.72533321 为分位参数； <math id="M43"><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195288&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195283&type= 3.47133350 3.72533321 为淘汰点。通常 <math id="M44"><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>∈</mo><mfenced open="[" close="]" separators="|"><mrow><mn mathvariant="normal">1.0,1.5</mn></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195295&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195292&type= 22.18266678 4.14866638 ， <math id="M45"><msub><mrow><mi>K</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>∈</mo><mfenced open="[" close="]" separators="|"><mrow><mn mathvariant="normal">2.5,3</mn></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195301&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195299&type= 19.47333336 4.14866638 。 在采用（8）式或（9）式解求出延迟相位改正数 <math id="M46"><mi mathvariant="normal">d</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195311&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195306&type= 6.85799980 4.40266657 基础上，便可求得同一时段内单颗卫星的延迟相位估值（ <math id="M47"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub><mo>=</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msubsup><mo>+</mo><mi mathvariant="normal">d</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195318&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195314&type= 26.24666595 4.74133301 ）。对于多系统多卫星组合，取同一时段内所有选定卫星的延迟相位估值的算术平均值作为该时段的延迟相位，亦即： <math display="block" id="M48"><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mi>l</mi></mrow></msubsup><mo>=</mo><mfrac><mrow><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow><mrow><mi>l</mi></mrow></msubsup><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow><mrow><mi>l</mi></mrow></msubsup><mo>+</mo><mo>⋯</mo><mo>+</mo><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><msub><mrow><mi mathvariant="normal">P</mi></mrow><mrow><mi>n</mi></mrow></msub></mrow><mrow><mi>l</mi></mrow></msubsup></mrow><mrow><mi>n</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195327&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195324&type= 45.29666519 8.72066593 式中， l 为时段编号； n 为该时段SNR MP 的个数。 对于单卫星组合，每个时段仅有一个延迟相位。由各时段延迟相位估值得延迟相位组合 <math id="M49"><mfenced open="{" close="}" separators="|"><mrow><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msubsup><mi mathvariant="normal">、</mi><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mi mathvariant="normal">、</mi><mo>⋯</mo><mi mathvariant="normal">、</mi><msubsup><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">P</mi></mrow><mrow><mi>m</mi></mrow></msubsup></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195333&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195330&type= 27.51666832 4.31799984 （上标1、2、…、 m 表示时段编号），据此表征土壤湿度的变化趋势。基于抗差估计的多系统多卫星组合GNSS-MR土壤湿度反演的主要数据处理流程如 图6 图6 所示。 3　数据采集 1 3.1　信噪比观测值的采集 2 信噪比观测值采集点位于贵州大学花溪校区西区校园内。实验区呈南北狭长带状走势，东西方向宽约60 m、南北方向长约320 m，此范围内无明显遮挡物和其他多路径反射源。东、西方向60 m外两侧为公路，且有建筑物和树木覆盖。此外，GNSS站点南侧土壤平整性较好，而北侧地形起伏较大，试验区总体呈北高南低态势。GNSS站点及周围环境如 图7 图7 所示。 信噪比采集时间为2017-05-08 T 14：30到2017-05-09 T 14：30。GNSS接收机及参数设置见 表1 表1 。 参数 参量设置 接收机 Trimble R8 接收机天线 内置天线 采样频率 10 s 截止高度角 0° 天线高 1.429 m（固定不变） GNSS卫星反射信号主要来自“第一菲涅尔区”，该区域内所有点对反射均有贡献（ 马建国等，2011 马建国等，2011 ）。“第一菲涅尔区”是1组与卫星高度角、方位角和天线高相关的椭圆，其半径计算式为 <math display="block" id="M50"><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>R</mi></mrow><mrow><mi mathvariant="normal">F</mi></mrow></msub><mo>=</mo><mroot><mrow><mi>λ</mi><mfrac><mrow><mi>h</mi></mrow><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>E</mi></mrow></mfrac></mrow><mrow/></mroot></mtd></mtr><mtr><mtd><msub><mrow><mi>L</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">M</mi><mi mathvariant="normal">A</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mi>R</mi></mrow><mrow><mi mathvariant="normal">F</mi></mrow></msub></mrow><mrow><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>E</mi></mrow></mfrac></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195359&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195355&type= 24.21466637 19.89666748 式中， <math id="M51"><msub><mrow><mi>R</mi></mrow><mrow><mi mathvariant="normal">F</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195364&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195361&type= 3.55599999 3.72533321 、 <math id="M52"><msub><mrow><mi>L</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">M</mi><mi mathvariant="normal">A</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195368&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195366&type= 5.84200001 3.72533321 分别为“第一菲涅尔区”椭圆的短半径和长半径； <math id="M53"><mi>λ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195372&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195370&type= 2.28600001 2.62466669 为载波波长； <math id="M54"><mi>h</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195375&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195373&type= 1.86266661 2.62466669 为信号垂直反射高度； <math id="M55"><mi>E</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195377&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195376&type= 2.28600001 2.62466669 为卫星高度角。实验站点卫星高度角5°和25°时的“第一菲涅尔区”如 图8 图8 所示： 1,2; 由 图8 图8 可看出，卫星高度角为5°时，“第一菲涅尔区”各方向距离GNSS站点不超过40 m；卫星高度角为25°时，反射点向GNSS站点靠拢，反射区最大距离为6 m左右。结合 图7 图7 进一步分析知，实验站点南侧“第一菲涅尔区”范围内均为单一土壤，可作为GNSS-MR土壤湿度反演的理想数据源。 3.2　烘焙法测量土壤湿度 2 土壤湿度采用烘焙法进行测量。土壤每两小时进行一次采样，取样深度为0—5 cm。将采集的土壤样本以密封袋密封保存并带回实验室进行烘焙测量。 烘焙法测定土壤湿度的原理是：在105°C温度下，土壤中水分开始蒸发，将土壤样品烘烤至恒重，通过烘烤前后质量之差，计算出土壤湿度SWC（Soil Water Content）。公式如下： <math display="block" id="M56"><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">C</mi><mo>=</mo><mfrac><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">B</mi><mi mathvariant="normal">B</mi><mo>-</mo><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">A</mi><mi mathvariant="normal">B</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">A</mi><mi mathvariant="normal">B</mi><mo>-</mo><mi mathvariant="normal">C</mi><mi mathvariant="normal">B</mi><mi mathvariant="normal">W</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195393&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195391&type= 36.23733139 7.70466709 式中，SWBB（Soil Weight、Before Baking）为烘烤前土壤重量，SWAB（Soil Weight、After Baking）；CBW（Cardboard Box Weighe）为纸盒重量。可看出，SWC无量纲。 4　算法的验证 1 4.1　实验方案 2 为检验基于抗差估计的多系统多卫星组合GNSS-MR土壤湿度反演的性能，设计如下4种方案。 方案1：基于多系统多卫星组合的GNSS-MR土壤湿度反演。 方案2：基于抗差的多系统多卫星组合GNSS-MR土壤湿度反演。 方案3：基于单卫星组合的GNSS-MR土壤湿度反演。 方案4：基于抗差的单卫星组合GNSS-MR土壤湿度反演。 4.2　土壤湿度反演 2 4.2.1　数据处理 3 基于多系统多卫星组合GNSS-MR土壤湿度反演的数据处理主要包括以下3个步骤： （1）GNSS-MR土壤湿度反演备选数据源的确定。多路径信号的振幅大小体现了反射信号的强弱，主要与反射体的特性和信号的频率有关（ 叶险峰，2016 叶险峰，2016 ），采用不同载波频率上SNR观测值的反演结果存在差异（ 陈发德 等，2018 陈发德 等，2018 ； Tabibi 等，2017 Tabibi 等，2017 ； 王泽民 等，2018 王泽民 等，2018 ）。鉴于1）GPS L1C信噪比较L2P信噪比土壤湿度反演精度高（ Vey 等，2016 Vey 等，2016 ）；2）BDS B1信噪比具有与L1C信噪比相似的雪深反演结果（ 王泽民 等，2018 王泽民 等，2018 ）；3）GLONASS系统R1（C/A码和P码）和R2（C/A码和P码）4种信号均能获得与GPS L2C相似的雪深反演结果，反演精度从高到低排列依次为：R2 > L2C > R1 > L1C（ Tabibi等，2017 Tabibi等，2017 ）；4）本实验所用接收机（Trimble R8）没有记录L2C观测值（RINEX VERSION 3.02），而本文的主要目的是为了验证基于抗差估计的多系统多卫星组合GNSS-MR算法在土壤湿度反演中的性能，故采用GPS L1C、GLONASS R1C、Galileo E1及BDS B1信噪比作为土壤湿度反演的数据源。 此外，考虑到实验站点南侧土壤比较平坦，将SNR观测值的方位角范围限制在130°—195°，以削弱土壤反射面平整性和土壤湿度的空间异值性对土壤湿度反演的影响。再者，受卫星运动周期特性的影响，表征多路径现象的SNR MP 也呈现出周期性的特点。经实验分析发现，当卫星高度角在5.75°—20.48°变化时，实验站点处SNR MP 的周期特性最为明显，而较强的周期特性有助于SNR MP 主频的提取。为此，对各时段各卫星波段1上的信噪比观测值通过设定卫星方位角区间（130°—195°）和卫星高度角区间（5.75°—20.48°），即得GNSS-MR土壤湿度反演的备选数据源（以下简称“S 1 ”）。 （2）GNSS-MR土壤湿度反演数据源的确定。构建二次多项式扣除S 1 中的趋势项即得多路径信噪比SNR MP 。SNR MP 中的隐含噪声将影响其频率的正确诊断。Lomb-Scargle谱分析可甄别出信号噪声水平的高低，并能从大噪声时间序列中准确地提取出弱周期信号，且探测效率较高（ 徐斌 等，2011 徐斌 等，2011 ）。此外，由 <math id="M57"><mi>f</mi><mo>=</mo><mfrac><mrow><mn mathvariant="normal">2</mn><mi>h</mi></mrow><mrow><mi>λ</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195396&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195395&type= 10.58333302 7.70466709 可得SNR MP 的概略频率。故为保证SNR MP 数据源的质量，在识别出各SNR MP 频率的基础上，对频率明显偏离其概值的SNR MP 进行剔除处理，从而最终确定出GNSS-MR土壤湿度反演的数据源。 图9 图9 为卫星G29、R09的功率图，其最大功率谱密度PSD（Power Spectrum Density）对应的频率分别为3.996和4.993，明显偏离其各自频率的概值15.02和15.27，故数据处理时剔除了G29和R09。 1,2; 基于单卫星组合GNSS-MR土壤湿度反演的实施步骤基本同上，区别在于每个时段均需解求出落入限制方位角范围的所有SNR MP 观测值的延迟相位，然后尝试从各时段众多延迟相位中分别挑选出一个，构成延迟相位组合，据此表征土壤湿度的变化。经（1）、（2）两个步骤处理后的选星结果如 表2 表2 所示。 时段 编号 选星前 多卫星 组合 单卫星 组合 1 G23、G03、C06 G23、G03、C06 C06 2 R15、E09、G07、R05、G30 R15、E09、G07、R05 E09 3 R16、E24、R09、R07、G28 E24、R07、G28 E24 4 G13、C07 G13、C07 C07 5 C07、R21 C07、R21 R21 6 E02、R11、G29、R22 E02、R11、R22 R11 7 C08、R23、G25、E22 G25、E22 G25 8 R13、G26、G21 R13、G21 R13 9 G18、R02、G27 G18、G27 G27 10 G10、C08、R17、E08、G08 G10、C08、R17、E08、G08 G08 11 C12、R03、E26、R04 C12、R03、E26 C12 12 G11、R19 G11、R19 G11 表2 表2 中，“选星前”表示限制卫星方位角后的所有卫星；“多卫星组合”为各时段分别剔除SNR MP 频率明显偏离其概值后的卫星组合；“单卫星组合”表示每个时段从“多卫星组合”中各挑选出一颗卫星后构成的卫星组合。选星结果表明，多卫星组合中每个时段可用的卫星数目仍可维持在2—5颗，足可保证GNSS-MR土壤湿度反演的正常实施。 （3）延迟相位的解求。在确定出土壤湿度反演数据源基础上，方案1、方案3依据 式（8） 式（8） 基于单位权解求延迟相位改正数；方案2和方案4则依据 式（9） 式（9） 基于IGGIII权函数采用抗差迭代解求延迟相位改正数。4种方案延迟相位与土壤湿度间的相关性如 图10 图10 所示。 4.2.2　回归分析建模 3 为验证各方案土壤湿度反演的性能，基于延迟相位和土壤湿度样本值进行建模，继而依据所建模型由延迟相位预测土壤湿度。因此，为获得土壤湿度和延迟相位间的函数关系，对延迟相位和土壤湿度进行一元线性回归分析建模。各方案土壤湿度回归效果见 图11 图11 所示散点图。 1,2;3,4;5; 图11 图11 中，“红色直线”表示土壤湿度模型回归值，“黑色小十字”为不同时段的土壤湿度测量值。可看出，土壤湿度测量值与方案2的模型回归值吻合最好，而与方案3的模型回归值偏离程度最大。各方案土壤湿度模型回归值残差及其置信区间如 图12 图12 所示。 图12 图12 中，横轴为时段编号，纵轴为土壤湿度回归值的残差及其置信区间。4种方案相关系数及回归效果统计见 表3 表3 。 方案 相关系数 F值（显著性水平0.001） F 0.001 （1，10）= 21.04 残差均方根误差（RMSE） 方案1 0.80 17.80 0.024 方案2 0.97 156.00 0.010 方案3 0.78 15.48 0.025 方案4 0.95 97.19 0.012 4.2.3　模型的验证 3 鉴于实验站点东、西两侧较窄，且存在有其他潜在的多路径反射源，而建模时采用的SNR观测值已限制为GNSS站点南侧。因此，为验证本文算法的优劣，将SNR观测值的方位角限制为GNSS站点的北侧（方位角为280°—335°），并兼顾各SNR MP 的频率是否异常、周期特性是否明显（高度角限制为8°—22°时，模型验证所用SNR MP 周期特性最为明显），得如 表4 表4 所示的选星结果。 时段编号 选星前 多卫星组合 单卫星组合 1 R16、R09 R09 R09 2 E02、R10 R10 R10 3 R11、G12 R11、G12 G12 4 G25、R22 R22 R22 5 G29 G29 G29 6 R23、G21 R23、G21 R23 7 R24、R03 R24、R03 R03 8 R17、E26、G22 R17、G22 G22 9 G03、C14、G23 G03、G23 G03 10 R05、G09、R06 G09、R06 G09 11 E19、G07 G07 G07 12 R06、G30 G30 G30 由 表4 表4 可看出，由于该站点北侧可视卫星数目有限，导致方位角280°—335°范围内可供使用的低高度角卫星较少。 图13 图13 中（a）、（b）分别为实验第1时段和第2时段卫星天空图。 1,2; 采用上述延迟相位的解求方法，结合 表4 表4 中的选星结果，在解求出各时段延迟相位后，将其代入对应的土壤湿度回归模型，从而计算出各时段土壤湿度的预报值。为衡量土壤湿度预报值与实测值间整体变化趋势的一致性，计算得各方案土壤湿度预报值与实测值间的相关系数（ 表5 表5 ）。需指出的是，由于各方案模型训练前后的土壤湿度分别相差一常数，故各方案模型训练前后的相关系数是相同的，所以 表5 表5 中只给出了一种相关系数。可看出，由基于抗差迭代解求的延迟相位（方案2和方案4）计算所得土壤湿度预报值与实测值间的相关性明显优于方案1和方案3，但由于实验时只采集了GNSS站点南侧的土壤湿度，各方案土壤湿度预报值与实测值均呈现出一定的系统性偏差。为提高模型预报精度，对土壤湿度回归模型进行训练，即各时段土壤湿度预报值均加上第1时段土壤湿度预报值与实测值之差，得 表5 表5 所示的统计信息。 表5 表5 中，RMSE1、RMSE2分别为模型训练前后的土壤湿度预报值残差的均方根误差。 方案 相关系数 RMSE1 RMSE2 方案1 0.73 0.057 0.029 方案2 0.92 0.122 0.016 方案3 0.76 0.047 0.026 方案4 0.91 0.088 0.023 4.3　结果分析 2 （1）由 表2 表2 分析可知，“多卫星组合”和“单卫星组合”的选星结果均涵盖了4大卫星导航系统，且各时段的卫星构成不尽相同。因此，若为获得较优的延迟相位估值，同时顾及SNR观测值的质量、多径环境的空间差异性和土壤湿度的空间异值性等因素的共同影响而进行选星，将可能导致基于单系统、单卫星的GNSS-MR无法实现土壤湿度的反演；此外，从各时段“多卫星组合”中分别挑选出1个SNR MP 组成单卫星SNR MP 组合，共有103680种不同的选法。因此，基于单卫星系统的GNSS-MR可靠性不高、稳定性较差；基于多系统单卫星组合的GNSS-MR选星较繁杂，实际可操作性不强，这无疑将限制GNSS-MR在环境监测中的扩展应用。 （2）由 图10 图10 知，无论是多系统多卫星组合GNSS-MR，还是多系统单卫星组合GNSS-MR，采用IGGIII权函数进行抗差迭代后，延迟相位与土壤湿度测量值间的相关性得到了显著改善。方案2较方案1相关系数由0.80提高到了0.97、方案4较方案3相关系数由0.78提高到了0.95。其原因可能是：SNR观测值受到了强噪声的污染，导致采用二次多项式去趋势项后提取的多路径信噪比不纯，加之最小二乘估计不具备鲁棒性，从而使得基于多路径信噪比采用最小二乘估计所求延迟相位的精度不高，继而降低了延迟相位与土壤湿度间的相关性；而基于IGGIII权函数的抗差估计在延迟相位解求过程中对受异常噪声污染的SNR MP 观测值对参数估值的贡献起到了很好的抑制作用。此外，方案1与方案3、方案2与方案4相关系数相当，但方案2避免了方案4的选星过程。 （3）回归方程中因变量 <math id="M58"><mi>y</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195467&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195465&type= 1.69333339 3.47133350 与自变量 <math id="M59"><mi>x</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195472&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20195470&type= 1.69333339 2.62466669 是否确实存在线性关系，需进行回归方程的显著性检验。取显著性水平0.001，查表得F 0.001 （1，10）= 21.04。由 表3 表3 知，对于方案1—方案4，检验量F分别为17.80、160.00、15.48和97.19。因17.80 < 21.04、15.48 < 21.04，故方案1和方案3所得回归方程不显著；而160.00 > 21.04、97.19 > 21.04，表明由方案2和方案4所得延迟相位和土壤湿度间存在较强的线性关系，且方案2较方案4回归方程显著性更强。此外，由 图12 图12 和 表3 表3 进一步分析知，就土壤湿度模型回归值残差的RMSE而言，四种方案中方案2最小，仅为0.010；而方案4又小于方案1和方案3。这主要得益于：方案2首先采用基于IGGIII权函数的抗差估计解求延迟相位，从而提高了各时段单个SNR MP 观测值的延迟相位解求精度；然后取同一时段中所有延迟相位抗差解的算术平均值作为该时段的延迟相位，使得每个时段延迟相位的精度得到了进一步的改善，故而依据各时段延迟相位所建立的土壤湿度回归模型更能表征土壤湿度的变化趋势。 （4）由 表5 表5 知，经抗差迭代后，相较于方案1和方案3，方案2和方案4土壤湿度预报值与实测值间的相关系数得到了明显改善，相关系数分别由0.73、0.76提升为0.92和0.91，且方案1与方案3、方案2与方案4的相关系数相当。但可能受土壤湿度的空间异值性、GNSS站点南北侧土壤成分构成的差异性以及北侧可供使用卫星数目较少等因素的共同作用，各方案土壤湿度预报值残差的均方根误差较大，但经模型训练后，预报精度得到了显著的提高。其中方案2的土壤湿度预报值与实测值吻合最好，具有最小的均方根误差。验证结果进一步证实了基于抗差估计的多系统多卫星组合GNSS-MR（方案2）较文中其他方案在土壤湿度反演方面的性能更优。 5　结论 1 针对多径环境往往存在空间差异性和单系统、单卫星短时间弧段内的重访周期较长，导致基于单系统、单卫星的GNSS-MR的土壤湿度反演的可靠性和稳定性较差，难以实现短期内土壤湿度反演的缺陷；同时顾及信噪比观测值难免会受到异常噪声的污染，而当前GNSS-MR土壤湿度反演所普遍采用的最小二乘估计不具鲁棒性，导致延迟相位解求精度不高的问题，本文在详细介绍GNSS-MR土壤湿度反演原理的基础上，鉴于多系统多卫星组合的SNR观测值能够同时提供更广的方位角覆盖范围，从而能够在顾及多径环境的空间差异性的基础上，筛选出更多的可用卫星，提出了一种基于抗差估计的多系统多卫星组合GNSS-MR土壤湿度反演算法，并结合实验结果对其在土壤湿度反演方面的性能进行了分析。实验结果表明，本文算法不仅能简化GNSS-MR土壤湿度反演时的繁杂选星过程，且可显著提高GNSS-MR土壤湿度反演的精度，并能改善基于单系统和多系统单卫星组合GNSS-MR土壤湿度反演的可靠性和稳定性，在土壤湿度反演中更具可操作性。 此外，本文实验同样存在几点不足：首先，由于多系统多卫星SNR观测值和土壤湿度数据源的限制，文中算法的性能还有待进一步验证；其次，文中实验只是将频率明显偏离其概值的信噪比观测值进行了简单的剔除处理，而未能对异常频率进行量化。因此，在后续的研究中，将采用更长时间跨度、更多样的反射环境进行实验，并考虑对异常频率进行量化，从而充分挖掘当前多模多频GNSS大数据蕴藏的全方位价值，使其能更好地服务于地表环境监测。 参考文献（References） Ao M S , Hu Y J , Liu Y D , Zeng Y and Wu Z . 2012 . Inversion of soil moisture fluctuation based on signal-to-noise ratio of global positioning system . Journal of Geomatics Science and Technology , 29 ( 2 ): 140 - 143 Inversion of soil moisture fluctuation based on signal-to-noise ratio of global positioning system 敖敏思 , 胡友健 , 刘亚东 , 曾云 , 吴正 . 2012 . GPS信噪比观测值的土壤湿度变化趋势反演 . 测绘科学技术学报 , 29 ( 2 ): 140- 143 [ DOI: 10.3969/j.issn.1673-6338.2012.02.015 http://dx.doi.org/10.3969/j.issn.1673-6338.2012.02.015 ] GPS信噪比观测值的土壤湿度变化趋势反演 Ao M S , Zhu J J , Hu Y J , Zeng Y and Liu Y D . 2015 . Comparative experiments on soil moisture monitoring with GPS SNR observations . Geomatics and Information Science of Wuhan University , 40 ( 1 ): 117 - 120, 127 Comparative experiments on soil moisture monitoring with GPS SNR observations 敖敏思 , 朱建军 , 胡友健 , 曾云 , 刘亚东 . 2015 . 利用SNR观测值进行GPS土壤湿度监测 . 武汉大学学报(信息科学版) , 40 ( 1 ): 117 - 120 , 127 [ DOI: 10.13203/j.whugis20130170 http://dx.doi.org/10.13203/j.whugis20130170 ] 利用SNR观测值进行GPS土壤湿度监测 Brocca L , Melone F , Moramarco T , Wagner W , Naeimi V , Bartalis Z and Hasenauer S . 2010 . Improving runoff prediction through the assimilation of the ASCAT soil moisture product . Hydrology and Earth System Sciences , 14 ( 10 ): 1881 - 1893 [ DOI: 10.5194/hess-14-1881-2010 http://dx.doi.org/10.5194/hess-14-1881-2010 ] Improving runoff prediction through the assimilation of the ASCAT soil moisture product Chen F D , Liu L L , Huang L K and Li J Y . 2018 . Sea surface height monitoring based on Multi-GNSS-MR . Progress in Geophysics , 33 ( 5 ): 1767 - 1772 Sea surface height monitoring based on Multi-GNSS-MR 陈发德 , 刘立龙 , 黄良珂 , 黎峻宇 . 2018 . 基于多模GNSS-MR海平面测高研究 . 地球物理学进展 , 33 ( 5 ): 1767- 1772 [ DOI: 10.6038/pg2018BB0412 http://dx.doi.org/10.6038/pg2018BB0412 ] 基于多模GNSS-MR海平面测高研究 Chew C C , Small E E , Larson K M and Zavorotny V U . 2015 . Vegetation sensing using GPS-interferometric reflectometry: theoretical effects of Canopy parameters on Signal-to-Noise ratio data . IEEE Transactions on Geoscience and Remote Sensing , 53 ( 5 ): 2755 - 2764 [ DOI: 10.1109/tgrs.2014.2364513 http://dx.doi.org/10.1109/tgrs.2014.2364513 ] Vegetation sensing using GPS-interferometric reflectometry: theoretical effects of Canopy parameters on Signal-to-Noise ratio data Jin S G , Zhang Q Y and Qian X D . 2017 . New progress and application prospects of Global Navigation Satellite System Reflectometry (GNSS+R) . Acta Geodaetica et Cartographica Sinica , 46 ( 10 ): 1389 - 1398 New progress and application prospects of Global Navigation Satellite System Reflectometry (GNSS+R) 金双根 , 张勤耘 , 钱晓东 . 2017 . 全球导航卫星系统反射测量(GNSS+R)最新进展与应用前景 . 测绘学报 , 46 ( 10 ): 1389- 1398 [ DOI: 10.11947/j.AGCS.2017.20170282 http://dx.doi.org/10.11947/j.AGCS.2017.20170282 ] 全球导航卫星系统反射测量(GNSS+R)最新进展与应用前景 Kerr Y H , Waldteufel P , Wigneron J P , Delwart S , Cabot F , Boutin J , Escorihuela M J , Font J , Reul N , Gruhier C , Juglea S E , Drinkwater M R , Hahne A , Martín-Neira M and Mecklenburg S . 2010 . The SMOS mission: new tool for monitoring key elements of the global water cycle . Proceedings of the IEEE , 98 ( 5 ): 666 - 687 . [ DOI: 10.1109/jproc.2010.2043032 http://dx.doi.org/10.1109/jproc.2010.2043032 ] The SMOS mission: new tool for monitoring key elements of the global water cycle Larson K M , Braun J J , Small E E , Zavorotny V U , Gutmann E D and Bilich A L . 2010 . GPS multipath and its relation to near-surface soil moisture content . IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing , 3 ( 1 ): 91 - 99 [ DOI: 10.1109/jstars.2009.2033612 http://dx.doi.org/10.1109/jstars.2009.2033612 ] GPS multipath and its relation to near-surface soil moisture content Larson K M , Small E E , Gutmann E , Bilich A , Axelrad P and Braun J . 2008 . Using GPS multipath to measure soil moisture fluctuations: initial results . GPS Solutions , 12 ( 3 ): 173 - 177 [ DOI: 10.1007/s10291-007-0076-6 http://dx.doi.org/10.1007/s10291-007-0076-6 ] Using GPS multipath to measure soil moisture fluctuations: initial results Ma J G , Cao K J , Zhang L and Li B . 2011 . Passive detection of marine targets based on GPS signals reflected from sea surface . Systems Engineering and Electronics , 33 ( 5 ): 987 - 991 Passive detection of marine targets based on GPS signals reflected from sea surface 马建国 , 曹可劲 , 张磊 , 李豹 . 2011 . 基于GPS信号的海面目标被动探测 . 系统工程与电子技术 , 33 ( 5 ): 987- 991 [ DOI: 10.3969/j.issn.1001-506X.2011.05.05 http://dx.doi.org/10.3969/j.issn.1001-506X.2011.05.05 ] 基于GPS信号的海面目标被动探测 Ray R L and Jacobs J M . 2007 . Relationships among remotely sensed soil moisture, precipitation and landslide events . Natural Hazards , 43 ( 2 ): 211 - 222 [ DOI: 10.1007/s11069-006-9095-9 http://dx.doi.org/10.1007/s11069-006-9095-9 ] Relationships among remotely sensed soil moisture, precipitation and landslide events Robock A , Vinnikov K Y , Srinivasan G , Entin J K , Hollinger S E , Speranskaya N A , Liu S X and Namkhai A . 2000 . The global soil moisture data bank . Bulletin of the American Meteorological Society , 81 ( 6 ): 1281 - 1300 [ DOI: 10.1175/1520-0477(2000)081<1281:TGSMDB>2.3.CO;2 http://dx.doi.org/10.1175/1520-0477(2000)081<1281:TGSMDB>2.3.CO;2 ] The global soil moisture data bank Roussel N , Frappart F , Ramillien G , Darrozes J , Baup F and Ha C . 2015 . Detection of soil moisture content changes by using a single geodetic antenna: the case of an agricultural plot // 2015 IEEE International Geoscience and Remote Sensing Symposium . Milan : IEEE : 2008 - 2011 [ DOI: 10.1109/igarss.2015.7326192 http://dx.doi.org/10.1109/igarss.2015.7326192 ] Detection of soil moisture content changes by using a single geodetic antenna: the case of an agricultural plot Roussel N , Frappart F , Ramillien G , Darrozes J , Baup F , Lestarquit L and Ha M C . 2016 . Detection of soil moisture variations using GPS and GLONASS SNR data for elevation angles ranging from 2° to 70° . IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing , 9 ( 10 ): 4781 - 4794 [ DOI: 10.1109/jstars.2016.2537847 http://dx.doi.org/10.1109/jstars.2016.2537847 ] Detection of soil moisture variations using GPS and GLONASS SNR data for elevation angles ranging from 2° to 70° Schaufler G , Kitzler B , Schindlbacher A , Skiba U , Sutton M A and Zechmeister-Boltenstern S . 2010 . Greenhouse gas emissions from European soils under different land use: effects of soil moisture and temperature . European Journal of Soil Science , 61 ( 5 ): 683 - 696 [ DOI: 10.1111/j.1365-2389.2010.01277.x http://dx.doi.org/10.1111/j.1365-2389.2010.01277.x ] Greenhouse gas emissions from European soils under different land use: effects of soil moisture and temperature Seneviratne S I , Corti T , Davin E L , Hirschi M , Jaeger E B , Lehner L , Orlowsky B and Teuling A J . 2010 . Investigating soil moisture-climate interactions in a changing climate: a review . Earth-Science Reviews , 99 ( 3/4 ): 125 - 161 [ DOI: 10.1016/j.earscirev.2010.02.004 http://dx.doi.org/10.1016/j.earscirev.2010.02.004 ] Investigating soil moisture-climate interactions in a changing climate: a review Seuffert G , Wilker H , Viterbo P , Mahfouf J F , Drusch M and Calvet J C . 2003 . Soil moisture analysis combining screen‐level parameters and microwave brightness temperature: a test with field data . Geophysical Research Letters , 30 ( 10 ): 1498 [ DOI: 10.1029/2003GL017128 http://dx.doi.org/10.1029/2003GL017128 ] Soil moisture analysis combining screen‐level parameters and microwave brightness temperature: a test with field data Strode P R R and Groves P D . 2016 . GNSS multipath detection using three-frequency signal-to-noise measurements . GPS Solutions , 20 ( 3 ): 399 - 412 [ DOI: 10.1007/s10291-015-0449-1 http://dx.doi.org/10.1007/s10291-015-0449-1 ] GNSS multipath detection using three-frequency signal-to-noise measurements Tabibi S , Geremia-Nievinski F and Van Dam T . 2017 . Statistical comparison and combination of GPS, GLONASS, and Multi-GNSS multipath reflectometry applied to snow depth retrieval . IEEE Transactions on Geoscience and Remote Sensing , 55 ( 7 ): 3773 - 3785 [ DOI: 10.1109/tgrs.2017.2679899 http://dx.doi.org/10.1109/tgrs.2017.2679899 ] Statistical comparison and combination of GPS, GLONASS, and Multi-GNSS multipath reflectometry applied to snow depth retrieval Vey S , Güntner A , Wickert J , Blume T and Ramatschi M . 2016 . Long-term soil moisture dynamics derived from GNSS interferometric reflectometry: a case study for Sutherland, South Africa . GPS Solutions , 20 ( 4 ): 641 - 654 [ DOI: 10.1007/s10291-015-0474-0 http://dx.doi.org/10.1007/s10291-015-0474-0 ] Long-term soil moisture dynamics derived from GNSS interferometric reflectometry: a case study for Sutherland, South Africa Wan W , Chen X W , Peng X F , Bai W H , Xia J M , Liang H , Zhang X M , Xiong P , Yang T , Cao Y C , Yi C , Zhao L M and Hong Y . 2016 . Overview and outlook of GNSS remote sensing technology and applications . Journal of Remote Sensing , 20 ( 5 ): 858 - 874 Overview and outlook of GNSS remote sensing technology and applications 万玮 , 陈秀万 , 彭学峰 , 白伟华 , 夏俊明 , 梁宏 , 张学民 , 熊攀 , 杨婷 , 曹云昌 , 尹聪 , 赵利民 , 洪阳 . 2016 . GNSS遥感研究与应用进展和展望 . 遥感学报 , 20 ( 5 ): 858- 874 [ DOI: 10.11834/jrs.20166228 http://dx.doi.org/10.11834/jrs.20166228 ] GNSS遥感研究与应用进展和展望 Wan W , Larson K M , Small E E , Chew C C and Braun J J . 2015 . Using geodetic GPS receivers to measure vegetation water content . GPS Solutions , 19 ( 2 ): 237 - 248 [ DOI: 10.1007/s10291-014-0383-7 http://dx.doi.org/10.1007/s10291-014-0383-7 ] Using geodetic GPS receivers to measure vegetation water content Wanders N , Karssenberg D , De Roo A , De Jong S M and Bierkens M F P . 2014 . The suitability of remotely sensed soil moisture for improving operational flood forecasting . Hydrology and Earth System Sciences Discussions , 10 ( 11 ): 13783 - 13816 [ DOI: 10.5194/hessd-10-13783-2013 http://dx.doi.org/10.5194/hessd-10-13783-2013 ] The suitability of remotely sensed soil moisture for improving operational flood forecasting Wang Z M , Liu Z K , An J C and Lin G B . 2018 . Snow depth detection and error analysis derived from SNR of GPS and BDS . Acta Geodaetica et Cartographica Sinica , 47 ( 1 ): 8 - 16 Snow depth detection and error analysis derived from SNR of GPS and BDS 王泽民 , 刘智康 , 安家春 , 林国标 . 2018 . 基于GPS和北斗信噪比观测值的雪深反演及其误差分析 . 测绘学报 , 47 ( 1 ): 8- 16 [ DOI: 10.11947/j.AGCS.2018.20160644 http://dx.doi.org/10.11947/j.AGCS.2018.20160644 ] 基于GPS和北斗信噪比观测值的雪深反演及其误差分析 Wu D , Wang L , Zhang Q , Fu W J , Shu B and Fan L H . 2015 . Implementation and verification analysis of the software of GNSS data quality check . Journal of Geomatics Science and Technology , 32 ( 4 ): 344 - 348 Implementation and verification analysis of the software of GNSS data quality check 吴丹 , 王利 , 张勤 , 付文举 , 舒宝 , 范丽红 . 2015 . GNSS数据质量评估软件的实现及验证分析 . 测绘科学技术学报 , 32 ( 4 ): 344- 348 [ DOI: 10.3969/j.issn.1673-6338.2015.04.004 http://dx.doi.org/10.3969/j.issn.1673-6338.2015.04.004 ] GNSS数据质量评估软件的实现及验证分析 Xu B , Yang T , Tan B H and Chen Y . 2011 . The simulate study of signal detection based on Lomb-Scargle algorithm . Nuclear Electronics and Detection Technology , 31 ( 6 ): 702 - 705 The simulate study of signal detection based on Lomb-Scargle algorithm 徐斌 , 杨涛 , 谭保华 , 陈益 . 2011 . 基于Lomb-Scargle算法的周期信号探测的模拟研究 . 核电子学与探测技术 , 31 ( 6 ): 702- 705 [ DOI: 10.3969/j.issn.0258-0934.2011.06.026 http://dx.doi.org/10.3969/j.issn.0258-0934.2011.06.026 ] 基于Lomb-Scargle算法的周期信号探测的模拟研究 Yang T , Wan W , Chen X W , Chu T X and Hong Y . 2017 . Using BDS SNR observations to measure near-surface soil moisture fluctuations: results from low vegetated surface . IEEE Geoscience and Remote Sensing Letters , 14 ( 8 ): 1308 - 1312 [ DOI: 10.1109/lgrs.2017.2710083 http://dx.doi.org/10.1109/lgrs.2017.2710083 ] Using BDS SNR observations to measure near-surface soil moisture fluctuations: results from low vegetated surface Yang Y X . 1993 . Robust Estimation and Its Applications . Beijing : Bayi Press Robust Estimation and Its Applications 杨元喜 . 1993 . 抗差估计理论及其应用 . 北京 : 八一出版社 抗差估计理论及其应用 Ye X F . 2016 . Research on Processing Method of GNSS Station Environment Error and Its Application Using SNR Data . Wuhan : China University of Geosciences (Wuhan) : 75 - 88 Research on Processing Method of GNSS Station Environment Error and Its Application Using SNR Data 叶险峰 . 2016 . 基于GNSS信噪比数据的测站环境误差处理方法及其应用研究 . 武汉 : 中国地质大学(武汉) : 75 - 88 基于GNSS信噪比数据的测站环境误差处理方法及其应用研究 Yin Z M , Meng F Y , Chen M J , Yan J Q and Wang W . 2016 . Quality analysis of satellite observation . GNSS World of China , 41 ( 1 ): 54 - 59 Quality analysis of satellite observation 尹子明 , 孟凡玉 , 陈明剑 , 闫建巧 , 汪威 . 2016 . 卫星导航数据质量分析 . 全球定位系统 , 41 ( 1 ): 54- 59 [ DOI: 10.13442/j.gnss.1008-9268.2016.01.010 http://dx.doi.org/10.13442/j.gnss.1008-9268.2016.01.010 ] 卫星导航数据质量分析 Yu K G , Ban W , Zhang X H and Yu X W . 2015 . Snow depth estimation based on multipath phase combination of GPS triple-frequency signals . IEEE Transactions on Geoscience and Remote Sensing , 53 ( 9 ): 5100 - 5109 [ DOI: 10.1109/tgrs.2015.2417214 http://dx.doi.org/10.1109/tgrs.2015.2417214 ] Snow depth estimation based on multipath phase combination of GPS triple-frequency signals Zavorotny V U , Larson L M , Braun J J , Small E E , Gutmann E D and Bilich A L . 2010 . A physical model for GPS multipath caused by land reflections: toward bare soil moisture retrievals . IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing , 3 ( 1 ): 100 - 110 [ DOI: 10.1109/jstars.2009.2033608 http://dx.doi.org/10.1109/jstars.2009.2033608 ] A physical model for GPS multipath caused by land reflections: toward bare soil moisture retrievals Zhang S C , Nan Y , Li Z Y , Zhang Q , Dai K Y and Zhao Y H . 2016 . Analysis of tide variation monitored by GNSS-MR . Acta Geodaetica et Cartographica Sinica , 45 ( 9 ): 1042 - 1049 Analysis of tide variation monitored by GNSS-MR 张双成 , 南阳 , 李振宇 , 张勤 , 戴凯阳 , 赵迎辉 . 2016 . GNSS-MR技术用于潮位变化监测分析 . 测绘学报 , 45 ( 9 ): 1042- 1049 [ DOI: 10.11947/j.AGCS.2016.20150498 http://dx.doi.org/10.11947/j.AGCS.2016.20150498 ] GNSS-MR技术用于潮位变化监测分析 10.11834/jrs.20219098 1007-4619（2021）06-1324-14 李婷，张显云，邓小东，李宏达，聂士海.2021.顾及多径环境差异和粗差的GNSS-MR土壤湿度反演.遥感学报，25（6）： 1324-1337 李婷，张显云，邓小东，李宏达，聂士海.2021.顾及多径环境差异和粗差的GNSS-MR土壤湿度反演.遥感学报，25（6）： 1324-1337 Li T，Zhang X Y，Deng X D，Li H D and Nie S H. 2021. GNSS-MR soil moisture retrieval considering the multipath environments differences and gross error. National Remote Sensing Bulletin， 25（6）：1324-1337 Li T，Zhang X Y，Deng X D，Li H D and Nie S H. 2021. GNSS-MR soil moisture retrieval considering the multipath environments differences and gross error. National Remote Sensing Bulletin， 25（6）：1324-1337 国家自然科学基金（编号：41701464，41901225）；贵州大学引进人才科研项目（编号：贵大人基合字（2016）51号） Supported by National Natural Science Foundation of China (No.41701464, 41901225); Guizhou University Introduced Talent Research Project (No. （2016）51) 2019-04-03 2021-06-07 2021-09-18 2020-02-19 遥感学报 6 25