CN113959329A - 一种基于多星数据融合的积雪深度反演方法 - Google Patents

Abstract

本发明涉及积雪深度监测技术领域，尤其涉及一种基于多星数据融合的积雪深度反演方法，包括获得监测点的GNSS观测数据；对GNSS数据进行预处理，获取观测时段各卫星的信噪比、高度角数据；分离出卫星高度角在5°~20°的信噪比和高度角数据，通过二次多项式拟合去除信噪比趋势项以获得信噪比残差序列；对信噪比残差序列进行Lomb‑Scargle频谱分析获得序列的主频率；获取各卫星反演的积雪深度；由获得的各卫星反演雪深结果建立GNSS‑IR雪深反演多星融合监测MARS模型；输出反演结果。本申请通过MARS理论建立多星数据融合积雪反演模型，能有效综合各卫星的反演结果，自动剔除对结果影响较大的数据，获得最优的反演卫星组合，大大提高了反演结果的可靠性和精度。

Description

一种基于多星数据融合的积雪深度反演方法

技术领域

本发明涉及积雪深度监测技术领域，尤其涉及一种基于多星数据融合的积雪深度反演方法。

背景技术

积雪是水文系统中的重要组成部分，它对淡水资源的平衡以及气候的调节过程都起着重要作用。因此持续、准确地监测地面积雪深度变化具有重要意义。传统的积雪监测方式主要是在地面上依靠人工手动测量或机器自动测量，但这些传统的方法存在外业工作量大、监测设备成本高、时空分辨率低等不足。

随着全球导航定位系统的研究深入，GNSS-IR(Global Navigation SatelliteSystem Interferometric Reflectometry)技术作为一种新型遥感监测技术逐渐发展起来。与传统方法相比，GNSS-IR技术具有反应迅速、精度高、成本低、全天时等优点，因此具有广阔的应用前景。但测站周围环境以及卫星自身的差异会影响卫星反演的结果，使得部分卫星反演的可靠性和精度降低，因此采用单颗卫星进行反演不确定性大，难以准确反映测站附近的雪深信息。

因此需要研发出一种基于多星数据融合的积雪深度反演方法来解决上述问题。

发明内容

本发明的目的就在于为了解决上述问题设计了一种基于多星数据融合的积雪深度反演方法。

本发明通过以下技术方案来实现上述目的：

一种基于多星数据融合的积雪深度反演方法，包括以下步骤：

S1、获得监测点的GNSS观测数据；

S2、对GNSS数据进行预处理，获取观测时段各卫星的信噪比、高度角数据；

S3、分离出卫星高度角在5°~20°的信噪比和高度角数据，通过二次多项式拟合去除信噪比趋势项以获得信噪比残差序列；

S4、对信噪比残差序列进行Lomb-Scargle频谱分析获得序列的主频率；

S5、获取各卫星反演的积雪深度；

S6、由获得的各卫星反演雪深结果建立GNSS-IR雪深反演多星融合监测MARS模型；

S7、输出反演结果。

具体地，步骤S2中获取的信噪比观测值表示为：

其中，SNR为信噪比观测值，A_c为合成信号，A_d为直射信号振幅，A_m为反射信号振幅，Ψ为直射信号和反射信号夹角。

具体地，步骤S2中还包括将信噪比线性化，计算公式如下：

其中，SNR_dB-Hz为接收机接收到的数值，SNR_vv为线性化后的信噪比数据。

具体地，步骤S3中信噪比残差序列的计算公式如下：

其中，h为垂直反射距离，λ为载波波长，θ为卫星高度角。

具体地，步骤S4中Lomb-Scargle公式可表示为：

其中，

表示离散观测序列的均值：

；σ²表示离散观测序列的方差：

；ω为角频率；t_i为离散观测值的时间；τ为相位；S(ω)是ω的功率；

其中，τ的计算公式为：

。

具体地，步骤S5中用以下公式获取积雪深度：

其中，λ为载波波长，H为天线相位中心到裸土的高度，h为积雪表面到天线相位中心的垂直距离，f为利用步骤5中Lomb-Scargle公式获取的反射信号的主频，

为积雪深度。

具体地，步骤S6中，MARS模型的建立过程首先是前向逐步，使用基函数将输入的雪深数据划分为不同区间，再在每个区间分别进行拟合，得到一个过拟合的模型，基函数定义为：

其中，S_km的取值为-1或1，表示在节点的左侧或右侧的样条函数，v(k,m)为独立变量的标识，t表示不同节点，x_v(k,m)表示不同节点左右的函数，x为雪深数据， t_km表示具体的某个节点位置，x-t_km和t_km-x分别为节点t右侧和左侧区域的截断样条函数，‘+’表示对负值取0，则由上两式得到MARS模型：

其中

为MARS模型输出的雪深反演值，a₀为常量参数，a_m为第m个基函数的系数，M为基函数的数量。

本发明的有益效果在于：

本申请与现有单星反演模型相比，通过MARS理论建立多星数据融合积雪反演模型，能有效综合各卫星的反演结果，剔除对结果影响较大的数据，获得最优的反演卫星组合，大大提高了反演结果的可靠性和精度。

附图说明

图1为本申请的流程框图；

图2本发明实施例中接收机接收到的一段信噪比数据图；

图3本发明实施例中信噪比序列去除趋势项后的残差序列；

图4本发明实施例中对残差序列进行Lomb-Scargle频谱分析得到的频谱图；

图5本发明实施例中部分卫星信号反射点轨迹图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明的描述中，需要理解的是，术语“上”、“下”、“内”、“外”、“左”、“右”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，或者是本领域技术人员惯常理解的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的设备或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

此外，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，“设置”、“连接”等术语应做广义理解，例如，“连接”可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接连接，也可以通过中间媒介间接连接，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

下面结合附图，对本发明的具体实施方式进行详细说明。

本实施例中，进行雪深反演所使用的卫星数据来自美国板块边界观测项目(PlateBoundary Observation,PBO)，测站编号为P351（43.8744°N，114.7192°W），该站站高为2m，2014年6月27日至2019年8月7日之间，卫星采样率为15s。；实测雪深数据来自美国农业部SNOTEL网络，监测站编号为490(43.87497°N，114.71363°W)。数据采集时间为2015年1月1日（年积日1）~2015年4月10日（年积日110），反演时利用年积日1~77的雪深数据建立多元自适应回归样条模型（MARS），利用年积日78~110的雪深数据进行反演。

如图1所示，一种基于多星数据融合的积雪深度反演方法，具体包括：

S1、利用GNSS接收机，接收监测点的GNSS观测数据；

S2、利用TEQC软件，对收集的GPS观测数据进行预处理，获取观测时段各卫星的信噪比、高度角数据，如图2所示，此时接收机接收到的信噪比观测值可表示为：

其中，SNR为信噪比观测值，A_c为合成信号，A_d为直射信号振幅，A_m为反射信号振幅，Ψ为直射信号和反射信号夹角。

S3、获取高度角在5°-20°范围完整弧段的信噪比和高度角数据，如图5所示，这个图可以直观地看出，并非所有记录的卫星数据都适合用来进行雪深反演，比如在右上角的卫星数据，除了6号卫星，其他卫星在5°-20°都有缺失，在反演时应该放弃这些数据，（同理正下方的3 6 12 也应该放弃）；利用以下公式将信噪比线性化：

其中，SNR_dB-Hz为接收机接收到的数值，SNR_vv为线性化后的信噪比数据。

进行二次多项式拟合，去除信噪比趋势项，获得信噪比残差序列（如图3所示）：

其中，h为垂直反射距离，λ为载波波长，θ为卫星高度角。

S4、对信噪比残差序列进行Lomb-Scargle频谱分析（得到的频谱图如图4所示），获取反射信号振幅峰值处频率，Lomb-Scargle公式可表示为：

其中，

表示离散观测序列的均值：

；σ²表示离散观测序列的方差：

；ω为角频率；S(ω)是ω的功率，即某一频率信号的功率谱值。

S5、获取峰值处频率f后，依次用以下公式获取积雪深度：

其中，λ为载波波长，H为天线相位中心到裸土的高度，h为积雪表面到天线相位中心的垂直距离。

S6、利用获取的各卫星的雪深数据，建立MARS模型。MARS模型的建立过程首先是前向逐步，即使用基函数将输入的雪深数据划分为不同区间，再在每个区间分别进行拟合，从而得到一个过拟合的模型，基函数定义为：

其中，x为雪深数据，t为节点位置，x-t_km和t_km-x分别为节点t右侧和左侧区域的样条函数，‘+’表示对负值取0，则由上两式可以得到MARS模型：

其中

为MARS模型输出的雪深反演值，a₀为常量参数，a_m为第m个基函数的系数，M为基函数的数量。

S7、输出反演结果。

前向逐步建立的模型是过拟合的，因此需要通过后向剪枝剔除对模型贡献度较小的基函数以获取最优基函数组合的模型，确定最优模型的原则是广义交叉验证（GCV）：

其中，λ为模型中项的个数，M(λ)为模型中有效参数的个数，N为基函数的个数，

为每一步最佳估计的模型数值。

最后通过交互验证确定最优模型，确定最优模型的指标为误差均方根(RMSECV)和q²：

本实施例中，以P351测站为试验站点，将单星反演模型、基于MRAS的多星融合模型进行对比：

模型	均方根误差/m	平均绝对误差/M	相关系数
				PRN1	0.051032	0.040255	0.745742
PRN10	0.053031	0.041496	0.826361
				PRN22	0.117893	0.1107	0.810472
PRN30	0.095782	0.085909	0.804396
				MARS	0.022294	0.0169	0.960303

由上表可以看出，单星反演模型不确定性较大、精度不够高，难以准确反映测站周围积雪深度。显然，利用本文提出的基于MARS理论的多星数据融合方法能显著提升反演的精度和可靠性，平均绝对误差(MAE)为0.0169，均方根误差(RMSE)为0.022294，并且与SNOTEL提供的当地实测雪深具有极高的相关性，相关系数为0.960303。

以上仅为本发明的较佳实施例，并不用于限制本发明，因此任何对本发明进行修改或等同替换而一切不脱离发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种基于多星数据融合的积雪深度反演方法，其特征在于，包括以下步骤：

S1、获得监测点的GNSS观测数据；

S2、对GNSS数据进行预处理，获取观测时段各卫星的信噪比、高度角数据；

S3、分离出卫星高度角在5°~20°的信噪比和高度角数据，通过二次多项式拟合去除信噪比趋势项以获得信噪比残差序列；

S4、对信噪比残差序列进行Lomb-Scargle频谱分析获得序列的主频率；

S5、获取各卫星反演的积雪深度；

S6、由获得的各卫星反演雪深结果建立GNSS-IR雪深反演多星融合监测MARS模型；

S7、输出反演结果。

2.根据权利要求1所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S2中获取的信噪比观测值表示为：

其中，SNR为信噪比观测值，A_c为合成信号，A_d为直射信号振幅，A_m为反射信号振幅，Ψ为直射信号和反射信号夹角。

3.根据权利要求2所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S2中还包括将信噪比线性化，计算公式如下：

其中，SNR_dB-Hz为接收机接收到的数值，SNR_vv为线性化后的信噪比数据。

4.根据权利要求1所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S3中信噪比残差序列的计算公式如下：

其中，h为垂直反射距离，λ为载波波长，θ为卫星高度角。

5.根据权利要求1所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S4中Lomb-Scargle公式可表示为：

其中，

表示离散观测序列的均值：

；σ²表示离散观测序列的方差：

；ω为角频率；t_i为离散观测值的时间；τ为相位；

S(ω)是ω的功率；

其中，τ的计算公式为：

。

6.根据权利要求5所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S5中用以下公式获取积雪深度：

其中，λ为载波波长，H为天线相位中心到裸土的高度，h为积雪表面到天线相位中心的垂直距离，f为利用步骤5中Lomb-Scargle公式获取的反射信号的主频，

为积雪深度。

7.根据权利要求1所述的一种基于多星数据融合的积雪深度反演方法，其特征在于，步骤S6中，MARS模型的建立过程首先是前向逐步，使用基函数将输入的雪深数据划分为不同区间，再在每个区间分别进行拟合，得到一个过拟合的模型，基函数定义为：

其中，S_km的取值为-1或1，表示在节点的左侧或右侧的样条函数，v(k,m)为独立变量的标识，t表示不同节点，x_v(k,m)表示不同节点左右的函数，x为雪深数据， t_km表示具体的某个节点位置，x-t_km和t_km-x分别为节点t右侧和左侧区域的截断样条函数，‘+’表示对负值取0，则由上两式得到MARS模型：

其中

为MARS模型输出的雪深反演值，a₀为常量参数，a_m为第m个基函数的系数，M为基函数的数量。

Images