CN106093339A - 一种基于gnss干涉功率峰谷值的土壤湿度测量方法 - Google Patents

Abstract

一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其步骤如下：1：信号接收处理；2：相关功率估计；3：剔除直射天线增益的影响；4：剔除直射信号的影响；5：归一化干涉功率降噪；6：提取归一化干涉功率峰谷值，估计反射系数；7：介电常数反演；根据反射系数与土壤相对介电常数和卫星仰角的关系，反演土壤相对介电常数；8：土壤湿度反演；根据土壤相对介电常数与土壤湿度间的关系，反演土壤湿度；通过以上步骤本发明给出了利用该峰谷值进行土壤相对介电常数与土壤湿度反演的方法，很大程度上减少了测量过程对经验模型的依赖，其算法复杂度低，时间分辨率高，普适性强。

Description

一种基于GNSS干涉功率峰谷值的土壤湿度测量方法

技术领域

本发明提供一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，它具体涉及一种基于全球卫星导航系统(以下简称GNSS)干涉功率峰谷值的土壤湿度测量方法，属于土壤湿度测量技术领域。

技术背景

在利用GNSS作为信号源进行土壤湿度测量的方法中，一般使用经土壤反射后的信号与直射信号的功率比值作为土壤反射率的估计，从而推算出土壤相对介电常数，最终得到土壤湿度，这种方法需要两根天线，一根为右旋圆极化(以下简称RHCP)天线，用于接收GNSS直射信号，一根为左旋圆极化(以下简称LHCP)天线，用于接收土壤反射信号，由于经过土壤反射的信号通常会比直射信号弱，因此接收反射信号的天线必须使用高增益天线，而高增益天线的造价比较高，从而增加了制造成本。近些年来出现了一种利用GNSS干涉信号进行土壤湿度测量的技术，称为GNSS-IR技术，该技术理论上只需要一根天线即可完成土壤湿度测量，如专利(大面积土壤湿度测量方法，专利申请号：2010102073992)使用单根线极化天线进行土壤湿度测量，但由于GNSS卫星发射的信号是RHCP极化的信号，若用单根线极化天线则接收能量将损失一半左右，从而使得信噪比降低，造成土壤湿度测量误差，另外该方法使用干涉信号拐点进行土壤湿度测量，其时间分辨率较低，一颗卫星的上升段信号观测值只能给出一个土壤湿度值。再如专利“基于北斗基准站接收的低仰角信号的土壤湿度反演方法”(专利申请号：2014102755487)，使用普通接收机天线接收直射与反射信号，这样能够克服线极化天线的缺点，但其土壤湿度反演算法利用的是干涉信号的初相，其算法复杂度高，时间分辨率较低，一颗卫星的上升段也只能给出一个土壤湿度值，且该方法利用的是干涉信号相位与土壤湿度之间的纯 经验模型，普适性较差。

发明内容

针对现有技术的不足，本发明目的在于，提供一种算法复杂度低，时间分辨率高的体积土壤湿度干涉测量方法，即利用GNSS-IR技术从GNSS接收机相关器输出的相关功率中提取归一化干涉功率，利用归一化干涉功率的振荡峰谷值进行土壤湿度反演，理论表明该方法能够给出反演的明确表达式，极大的减少了对经验模型的依赖，普适性提高；由于一颗卫星上升段归一化干涉功率曲线存在多个峰值与谷值，因此能够输出多个土壤湿度值，提高了时间分辨率。

本发明一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其步骤如下:

步骤1：信号接收处理

单根低增益RHCP天线接收直射信号与土壤反射信号，信号通过射频模块下变频到数字中频信号，经接收机捕获跟踪后，从相关器输出相关值；

步骤2：相关功率估计

通常相关器输出含有较大噪声，为了更加准确的估计相关功率，运用一定的信号估计算法估计相关功率，形成相关功率观测值序列；

步骤3：剔除直射天线增益的影响

取卫星上升段0～90度范围内的相关功率观测值，从中剔除天线增益对干涉信号的影响；

所述的“天线增益”可通过测量及仿真得到；

步骤4：剔除直射信号的影响

估计直射信号功率，然后从步骤3得到的观测值序列中将其剔除，得到含有噪声的归一化干涉功率；

步骤5：归一化干涉功率降噪

取步骤4中卫星低仰角范围内的归一化干涉功率观测值序列进行滤波或降噪，得到最终用于反演的归一化干涉功率；

步骤6：提取归一化干涉功率峰谷值，估计反射系数

运用峰谷值搜索算法搜索归一化干涉功率的峰谷值，利用这些功率值和对应的卫星仰角与天线增益信息计算反射系数；

步骤7：介电常数反演

根据反射系数与土壤相对介电常数和卫星仰角的关系，反演土壤相对介电常数；

步骤8：土壤湿度反演

根据土壤相对介电常数与土壤湿度间的关系，反演土壤湿度。

其中，在步骤1中所述的“相关器输出相关值”可用复信号表示为：

上式假设仅有一条土壤镜面反射路径，其中i、q表示相关器IQ两路输出，A为直射信号到达天线前端时的平均功率，G(·)为天线增益，θ为卫星仰角，R(·)为归一化码自相关函数，τ_i、分别为码相位和载波相位，i取0、1、x时分别表示直射信号、反射信号、接收机本地信号，|Γ|为右旋-右旋(RR)反射系数模值，为反射系数相角，为复数噪声项，其实部与虚部可用独立同分布的高斯白噪声描述；

其中，在步骤2中所述的“相关功率”，其计算方法如下：

在天线架设高度较低的情况下，直射与反射信号间的多径延时可忽略，因此τ₀≈τ₁≈τ_x，R(·)≈1，在忽略噪声的情况下，相关功率可用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。

其中，在步骤2中所述的“信号估计算法”，是以极大似然(ML)估计算法为优选，该算法假设不考虑导航电文位的影响，在这种情况下相关功率观测值的估计式如下：

${\hat{P}}_i=\frac{1}{M^2}{\left|\underset{k=1}{\overset{M}{\Σ}}{\tilde{p}}_{i,k}\right|}^2---\left(3\right)$

其中为相关功率估计值，M为参与估计的相关器输出值的个数。

其中，在步骤4中所述的“估计直射信号功率”的方法为：在剔除直射天线增益的影响后，取90度仰角附近无明显震荡的多个观测值的平均值得到，于是所述“归一化干涉功率”估计值可用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。

其中，在步骤5中所述的“低仰角范围”，是以35度以下为优选，所述“滤波或降噪”方法，是以“Lomb-Scargle算法辅助RLS滤波算法”为优选，Lomb-Scargle算法用来估计振荡频率。

其中，在步骤6中所述的“峰谷值搜索算法”，是以“可变尺度峰值检测算法”(Automatic Multiscale-based Peak Detection以下简称AMPD算法)为优选，其步骤如下：

步骤6.1：设步骤5得到的观测值序列记为X＝[x₁,x₂,...,x_N-1,x_N]，建立一个L×N的矩阵M，矩阵M的每个元素的定义如下：

其中L为窗的个数， 为不小于z的最小整数，窗的长度w_k为{w_k＝2k|k＝1,2...,L}，r是在[0,1]之间均匀分布的随机数，α为一个常数(可取1)。

步骤6.2：对矩阵M每一行进行求和运算：

${\mathrm{\γ}}_k=\underset{i=1}{\overset{N}{\Σ}}{m}_{k,i},k\∈\{1,2,\mathrm{...},L\}---\left(6\right)$

得到一个L维列向量，求出该列向量的最小值所对应的下标，记为λ＝argmin(γ_k)并取矩阵M的前λ行，构成一个新矩阵M_r。

步骤6.3：对矩阵M_r求每一列的标准差：

${\mathrm{\σ}}_i=\frac{1}{\mathrm{\λ}-1}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{\[{(m_{k,i}-\frac{1}{\mathrm{\λ}}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{m}_{k,i})}^2\]}^{\frac{1}{2}},i\∈\{1,2,...,N\}---\left(7\right)$

则所有σ_i＝0的元素所对应的下标就是原始序列峰值所在位置，于是可从序列中提取峰值。

步骤6.4：将步骤6.1中公式(5)改为：

重复步骤6.2至6.3即可得到原始序列谷值所在位置，于是可从序列中提取谷值。

其中，在步骤6中所述“估计反射系数”，其计算公式如下：

式(4)当相位差时，归一化干涉功率达到局部峰值P_peak， 时，归一化干涉功率达到局部谷值P_valley，于是有：

$\left|{\mathrm{\Γ}}_1\right|=(\sqrt{P_{peak}}-1)\sqrt{\frac{G\left({\mathrm{\θ}}_1\right)}{G(-{\mathrm{\θ}}_1)}}---\left(9\right)$

$\left|{\mathrm{\Γ}}_2\right|=(1-\sqrt{P_{valley}})\sqrt{\frac{G\left({\mathrm{\θ}}_2\right)}{G(-{\mathrm{\θ}}_2)}}---\left(10\right)$

其中θ₁为局部峰值对应的卫星高度角，θ₂为局部谷值对应的卫星高度角，Γ₁、Γ₂为相应的反射系数，G(·)为天线增益。

其中，在步骤7中所述“反射系数与土壤相对介电常数和卫星仰角的关系”，可用下式描述：

$\mathrm{\Γ}=\frac{(1-{\mathrm{\ϵ}}_r){\cos }^2\mathrm{\θ}}{(sin\mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})({\mathrm{\ϵ}}_r\sin \mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})}---\left(11\right)$

上式假设土壤为良介质，其中ε_r为土壤相对介电常数，θ为卫星仰角。

其中，在步骤8中所述“土壤相对介电常数与土壤湿度间的关系”，可利用现有的经验模型，如Wang模型，其表达式如下：

ε_r＝3.1+17.36·SM+63.12·SM² (12)

其中ε_r为土壤相对介电常数，SM为体积土壤湿度单位为cm³/cm³，该经 验模型也可以自行通过其他仪器分别同时测量土壤相对介电常数和土壤湿度，通过数据拟合的方法建立两者之间的关系模型。

通过以上步骤实现了从GNSS直射信号与反射信号的干涉效应中提取归一化干涉功率峰谷值，并给出了利用该峰谷值进行土壤相对介电常数与土壤湿度反演的方法，以上绝大多数步骤都有明确的理论支撑，很大程度上减少了测量过程对经验模型的依赖。

和现有技术相比，本发明具有如下有以下有益效果：

提供了一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其算法复杂度低，时间分辨率高，普适性强。

附图说明

图1本发明所述测量方法流程图。

图2 GNSS干涉场景。

图3相关功率估计值。

图4未经滤波的归一化干涉功率。

图5 Lomb-Scargle算法辅助RLS滤波算法框图。

图6经降噪后的归一化干涉功率曲线及AMPD算法峰谷值提取结果。

具体实施方式

下面结合附图对本发明实施例做进一步描述：

本发明一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其流程如图1所示，包括如下步骤:

步骤1：信号接收处理

单根低增益RHCP天线接收直射信号与土壤反射信号，其场景如图2所示，信号通过射频模块下变频到数字中频信号，经接收机捕获跟踪后，从相关器输出 相关值；

步骤2：相关功率估计

通常相关器输出含有较大噪声，为了更加准确的估计相关功率，运用一定的信号估计算法估计相关功率，形成相关功率观测值序列；

步骤3：剔除直射天线增益的影响

取卫星上升段0～90度范围内的相关功率观测值，从中剔除天线增益对干涉信号的影响。

所述的“天线增益”可通过测量及仿真得到；

步骤4：剔除直射信号的影响

估计直射信号功率，然后从步骤3得到的观测值序列中将其剔除，得到含有噪声的归一化干涉功率；

步骤5：归一化干涉功率降噪

取步骤4中卫星低仰角范围内的归一化干涉功率观测值序列进行滤波或降噪，得到最终用于反演的归一化干涉功率；

步骤6：提取归一化干涉功率峰谷值，估计反射系数

运用峰谷值搜索算法搜索归一化干涉功率的峰谷值，利用这些功率值和对应的卫星仰角与天线增益信息计算反射系数。

步骤7：介电常数反演

根据反射系数与土壤相对介电常数和卫星仰角的关系，反演土壤相对介电常数；

步骤8：土壤湿度反演

根据土壤相对介电常数与土壤湿度间的关系，反演土壤湿度。

其中，在步骤1中所述的“相关器输出相关值”可用复信号表示为：

上式假设仅有一条土壤镜面反射路径，其中i、q表示相关器IQ两路输出，A为直射信号到达天线前端时的平均功率，G(·)为天线增益，θ为卫星仰角，R(·)为归一化码自相关函数，τ_i、分别为码相位和载波相位，i取0、1、x时分别表示直射信号、反射信号、接收机本地信号，|Γ|为右旋-右旋(RR)反射系 数模值，为反射系数相角，为复数噪声项，其实部与虚部可用独立同分布的高斯白噪声描述；

其中，在步骤2中所述的“相关功率”计算方法如下：

在天线架设高度较低的情况下，直射与反射信号间的多径延时可忽略，因此τ₀≈τ₁≈τ_x，R(·)≈1，在忽略噪声的情况下，相关功率可用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。

其中，在步骤2中所述的“信号估计算法”，是以极大似然(ML)估计算法为优选，该算法假设不考虑导航电文位的影响，在这种情况下相关功率观测值的估计式如下：

${\hat{P}}_i=\frac{1}{M^2}{\left|\underset{k=1}{\overset{M}{\Σ}}{\tilde{p}}_{i,k}\right|}^2---\left(15\right)$

其中为相关功率估计值，M为参与估计的相关器输出值的个数，假设在此期间干涉信号的各参量近似不变，由于(15)式的平均作用，干涉信号的噪声将降低，如图3所示。

其中，在步骤4中所述的“估计直射信号功率”A的方法为：在剔除直射天线增益的影响后，取90度仰角附近无明显震荡的多个观测值的平均值得到，于是所述“归一化干涉功率”估计值可用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。“归一化干涉功率”形如图4所示。

其中，在步骤5中所述的“低仰角范围”，是以35度以下为优选，所述“滤波或降噪”方法，是以“Lomb-Scargle算法辅助RLS滤波算法”为优选， Lomb-Scargle算法用来估计振荡频率，算法原理如下：

式(16)可写成如下形式：其中A、B、C为待估计的参数，其通过RLS滤波算法进行估计，f可预先使用Lomb-Scargle算法估计，整个滤波过程可用图5表示。

其中，在步骤6中所述“峰谷值搜索算法”，是以“可变尺度峰值检测算法”(Automatic Multiscale-based Peak Detection以下简称AMPD算法)为优选，其步骤如下：

步骤6.1：设步骤5得到的观测值序列记为X＝[x₁,x₂,...,x_N-1,x_N]，建立一个L×N的矩阵M，矩阵M的每个元素的定义如下：

其中L为窗的个数， 为不小于z的最小整数，窗的长度w_k为{w_k＝2k|k＝1,2...,L}，r是在[0,1]之间均匀分布的随机数，α为一个常数(可取1)。

步骤6.2：对矩阵M每一行进行求和运算：

${\mathrm{\γ}}_k=\underset{i=1}{\overset{N}{\Σ}}{m}_{k,i},k\∈\{1,2,...,L\}---\left(18\right)$

得到一个L维列向量，求出该列向量的最小值所对应的下标，记为λ＝argmin(γ_k)并取矩阵M的前λ行，构成一个新矩阵M_r。

步骤6.3：对矩阵M_r求每一列的标准差：

${\mathrm{\σ}}_i=\frac{1}{\mathrm{\λ}-1}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{\[{(m_{k,i}-\frac{1}{\mathrm{\λ}}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{m}_{k,i})}^2\]}^{\frac{1}{2}},i\∈\{1,2,...,N\}---\left(19\right)$

则所有σ_i＝0的元素所对应的下标就是原始序列峰值所在位置，于是可从序列中提取峰值。

步骤6.4：将步骤6.1中公式(17)改为：

重复步骤6.2至6.3即可得到原始序列谷值所在位置，于是可从序列中提取谷值。

其中步骤6所述“反射系数”计算公式如下：

式(16)当相位差时，归一化干涉功率达到局部峰值P_peak， 时，归一化干涉功率达到局部谷值P_valley，于是有：

$\left|{\mathrm{\Γ}}_1\right|=(\sqrt{P_{peak}}-1)\sqrt{\frac{G\left({\mathrm{\θ}}_1\right)}{G(-{\mathrm{\θ}}_1)}}---\left(21\right)$

$\left|{\mathrm{\Γ}}_2\right|=(1-\sqrt{P_{valley}})\sqrt{\frac{G\left({\mathrm{\θ}}_2\right)}{G(-{\mathrm{\θ}}_2)}}---\left(22\right)$

其中θ₁为局部峰值对应的卫星高度角，θ₂为局部谷值对应的卫星高度角，Γ₁、Γ₂为相应的反射系数，G(·)为天线增益，峰谷值搜索结果如图6所示,从图中可以看出卫星上升段归一化干涉功率上升段存在多个峰谷值，均可用于反演，因此时间分辨率大大提高。

其中，在步骤7中所述“反射系数与土壤相对介电常数和卫星仰角的关系”，可用下式描述：

$\mathrm{\Γ}=\frac{(1-{\mathrm{\ϵ}}_r){\cos }^2\mathrm{\θ}}{(sin\mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})({\mathrm{\ϵ}}_r\sin \mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})}---\left(23\right)$

上式假设土壤为良介质，其中ε_r为土壤相对介电常数，θ为卫星仰角。

其中，在步骤8中所述“土壤相对介电常数与土壤湿度间的关系”，可利用现有的经验模型，如wang模型，其表达式如下：

ε_r＝3.1+17.36·SM+63.12·SM² (24)

其中ε_r为土壤相对介电常数，SM为体积土壤湿度单位为cm³/cm³，该经验模型也可以自行通过其他仪器分别同时测量土壤相对介电常数和土壤湿度，通过数据拟合的方法建立两者之间的关系模型。

Claims (10)

1.一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：其步骤如下:

步骤1：信号接收处理

单根低增益RHCP天线接收直射信号与土壤反射信号，信号通过射频模块下变频到数字中频信号，经接收机捕获跟踪后，从相关器输出相关值；

步骤2：相关功率估计

通常相关器输出含有较大噪声，为了更加准确的估计相关功率，运用一定的信号估计算法估计相关功率，形成相关功率观测值序列；

步骤3：剔除直射天线增益的影响

取卫星上升段0～90度范围内的相关功率观测值，从中剔除天线增益对干涉信号的影响；所述的天线增益是通过测量及仿真得到；

步骤4：剔除直射信号的影响

估计直射信号功率，然后从步骤3得到的观测值序列中将其剔除，得到含有噪声的归一化干涉功率；

步骤5：归一化干涉功率降噪

取步骤4中卫星低仰角范围内的归一化干涉功率观测值序列进行滤波或降噪，得到最终用于反演的归一化干涉功率；

步骤6：提取归一化干涉功率峰谷值，估计反射系数

运用峰谷值搜索算法搜索归一化干涉功率的峰谷值，利用这些功率值和对应的卫星仰角与天线增益信息计算反射系数；

步骤7：相对介电常数反演

根据反射系数与土壤相对介电常数和卫星仰角的关系，反演土壤相对介电常数；

步骤8：土壤湿度反演

根据土壤相对介电常数与土壤湿度间的关系，反演土壤湿度；

通过以上步骤实现了从GNSS直射信号与反射信号的干涉效应中提取归一化干涉功率峰谷值，并给出了利用该峰谷值进行土壤相对介电常数与土壤湿度反演的方法，以上绝大多数步骤都有明确的理论支撑，很大程度上减少了测量过程对经验模

型的依赖，其算法复杂度低，时间分辨率高，普适性强。

2.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤1中所述的“相关器输出相关值”可用复信号表示为：

上式假设仅有一条土壤镜面反射路径，其中i、q表示相关器IQ两路输出，A为直射信号到达天线前端时的平均功率，G(·)为天线增益，θ为卫星仰角，R(·)为归一化码自相关函数，τ_i、分别为码相位和载波相位，i取0、1、x时分别表示直射信号、反射信号、接收机本地信号，|Γ|为右旋-右旋(RR)反射系数模值，为反射系数相角，为复数噪声项，其实部与虚部用独立同分布的高斯白噪声描述。

3.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤2中所述的“相关功率”，其计算方法如下：

在天线架设高度较低的情况下，直射与反射信号间的多径延时被忽略，因此τ₀≈τ₁≈τ_x，R(·)≈1，在忽略噪声的情况下，相关功率用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。

4.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤2中所述的“信号估计算法”，是以极大似然即ML估计算法为优选，该算法假设不考虑导航电文位的影响，在这种情况下相关功率观测值的估计式如下：

${\hat{P}}_i=\frac{1}{M^2}{\left|\underset{k=1}{\overset{M}{\Σ}}{\tilde{p}}_{i,k}\right|}^2---\left(3\right)$

其中为相关功率估计值，M为参与估计的相关器输出值的个数。

5.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤4中所述的“估计直射信号功率”的方法为：在剔除直射天线增益的影响后，取90度仰角附近无明显震荡的多个观测值的平均值得到，于是所述“归一化干涉功率”估计值用下式表示：

其中 为反射信号相对于直射信号的路径延迟导致的相位延迟，为反射系数相角，θ为卫星高度角，H为天线相位中心到地面的高度，Γ为反射系数。

6.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤5中所述的“低仰角范围”，是以35度以下为优选，所述“滤波或降噪”方法，是以“Lomb-Scargle算法辅助RLS滤波算法”为优选，Lomb-Scargle算法用来估计振荡频率。

7.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤6中所述“峰谷值搜索算法”，是以“可变尺度峰值检测算法”即“AMPD算法”为优选，其步骤如下：

步骤6.1：设步骤5得到的观测值序列记为X＝[x₁,x₂,...,x_N-1,x_N]，建立一个L×N的矩阵M，矩阵M的每个元素的定义如下：

其中L为窗的个数， 为不小于z的最小整数，窗的长度w_k为{w_k＝2k|k＝1,2...,L}，r是在[0,1]之间均匀分布的随机数，α为一个常数(可取1)；

步骤6.2：对矩阵M每一行进行求和运算：

${\mathrm{\γ}}_k=\underset{i=1}{\overset{N}{\Σ}}{m}_{k,i},k\∈\{1,2,...,L\}---\left(6\right)$

得到一个L维列向量，求出该列向量的最小值所对应的下标，记为λ＝argmin(γ_k)，并取矩阵M的前λ行，构成一个新矩阵M_r；

步骤6.3：对矩阵M_r求每一列的标准差：

${\mathrm{\σ}}_i=\frac{1}{\mathrm{\λ}-1}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{\[{(m_{k,i}-\frac{1}{\mathrm{\λ}}\underset{k=1}{\overset{\mathrm{\λ}}{\Σ}}{m}_{k,i})}^2\]}^{\frac{1}{2}},i\∈\{1,2,...,N\}---\left(7\right)$

则所有σ_i＝0的元素所对应的下标就是原始序列峰值所在位置，于是可从序列中提取峰值；

步骤6.4：将步骤6.1中公式(5)改为：

重复步骤6.2至6.3即可得到原始序列谷值所在位置，于是可从序列中提取谷值。

8.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤6中所述“估计反射系数”，其计算公式如下：

式(4)当相位差时，归一化干涉功率达到局部峰值P_peak，时，归一化干涉功率达到局部谷值P_valley，于是有：

$\left|{\mathrm{\Γ}}_1\right|=(\sqrt{P_{peak}}-1)\sqrt{\frac{G\left({\mathrm{\θ}}_1\right)}{G(-{\mathrm{\θ}}_1)}}---\left(9\right)$

$\left|{\mathrm{\Γ}}_2\right|=(1-\sqrt{P_{valley}})\sqrt{\frac{G\left({\mathrm{\θ}}_2\right)}{G(-{\mathrm{\θ}}_2)}}---\left(10\right)$

其中θ₁为局部峰值对应的卫星高度角，θ₂为局部谷值对应的卫星高度角，Γ₁、Γ₂为相应的反射系数，G(·)为天线增益。

9.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤7中所述“反射系数与土壤相对介电常数和卫星仰角的关系”，用下式描述：

$\mathrm{\Γ}=\frac{(1-{\mathrm{\ϵ}}_r){\cos }^2\mathrm{\θ}}{(sin\mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})({\mathrm{\ϵ}}_{r}sin\mathrm{\θ}+\sqrt{({\mathrm{\ϵ}}_r-{\cos }^2\mathrm{\θ})})}---\left(11\right)$

上式假设土壤为良介质，其中ε_r为土壤相对介电常数，θ为卫星仰角。

10.根据权利要求1所述的一种基于GNSS干涉功率峰谷值的土壤湿度测量方法，其特征在于：在步骤8中所述“土壤相对介电常数与土壤湿度间的关系”，利用现有的经验模型，如Wang模型，其表达式如下：

ε_r＝3.1+17.36·SM+63.12·SM² (12)

其中ε_r为土壤相对介电常数，SM为体积土壤湿度单位为cm³/cm³，该经验模型也可以自行通过其他仪器分别同时测量土壤相对介电常数和土壤湿度，通过数据拟合的方法建立两者之间的关系模型。