CN105277575A - 一种利用gps信噪比等级数据进行土壤湿度估计的方法 - Google Patents

Abstract

本发明公开了一种利用GPS信噪比等级数据sn_rnx进行土壤湿度估计的方法，该方法根据sn_rnx自身呈矩形波变换的特点对数据进行筛选后根据频谱分析法和最小二乘法提取相位信息，最后建立相位与湿度的函数关系并依据此函数关系反演土壤湿度等步骤；本发明的特点是首次采用了CORS站的GPS信噪比等级数据sn_rnx进行土壤湿度的反演，避免了对原始信噪比数据sn_raw的依赖，这样能够充分发挥CORS站历史数据的作用进行土壤湿度估计，利用率高。

Description

一种利用GPS信噪比等级数据进行土壤湿度估计的方法

技术领域

本发明属于电子信息技术领域，涉及一种基于CORS站的数据反演土壤湿度方法，尤其涉及一种用GPS信噪比等级数据sn_rnx进行土壤湿度估计的方法。

背景技术

土壤湿度即土壤含水量，是土壤的重要物理物质。早已成为水文学、气象学、生态学、农业以及地球科学中一个重要的研究对象，其准确和长期监测是环境科学研究的基础。

目前土壤湿度的测量方法主要有下面几种，它们有着各种优点和缺点。

传统的实测法包括烘干法，土壤湿度计法，中子仪法，电阻法等(张晓虎,李新平.几种常用土壤含水量测定方法的研究进展.陕西农业科学.2008,6:114-117)。这些方法都是接触式法，可以忽略土壤类型、植被以及环境的影响，准确测量土体剖面的含水量，但实测范围有限，无法进行大范围的连续实时性观测，费时费力，数据得不到及时更新。

为了解决传统测量法无法获得大范围的土壤湿度，有学者开展采用微波遥感方法反演土壤湿度(例如赵少华、杨永辉，基于双时相ASAR影像的土壤湿度反演研究，农业工程学报，2008年6月，184-188)利用土壤表面发射反演出土壤水分。但是该方法受限于SAR卫星的重返周期，数据的时间连续性不高，同一地方需几日才能获得一场数据。

目前一些学者开展了利用GPS反射信号测量土壤湿度的研究。美国进行了利用天顶天线和下视天线的实验(毛克彪、王建明、张孟阳，GNSS-R信号反演土壤水分研究分析，遥感信息2009，3)，该方法需要采用特殊设计的延迟接收机。武汉大学曾利用天线和接收机进行过反射信号和直达信号的比值研究来反演土壤湿度(严颂华，张训械，基于GNSS-R信号的土壤湿度反演研究，电波科学学报，2010年第1期)，但该方法主要利用的是直达信号及反射信号功率的比值，容易受到地面粗糙度的影响。另一种方法是采用GPS干涉信号的差分方法，如，采用干涉信号中的信噪比数据来实现土壤湿度的反演，虽然反演得到的土壤湿度精度高(申请号为2014102755487的发明专利：基于北斗基准站接收的低仰角信号的土壤湿度反演方法)，但是信噪比数据在定位时仅仅作为校正使用，因此绝大多数的CORS站并不提供信噪比数据。可见，如何充分利用现成的CORS站点的历史数据反演土壤湿度仍然需要探索。

发明内容

为了克服现有技术的局限性，本发明首次提出利用GPSCORS站的信噪比等级数据sn_rnx，建立其相位与湿度的函数关系，来实现土壤湿度的反演，本方法应用于农业或气象研究。

本发明所采用的技术方案是：一种利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于，包括以下步骤：

1.一种利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于，包括以下步骤：

步骤1：GPS信噪比等级数据、仰角以及方位角的获取；

利用地基CORS站的GPS接收机接收GPS干涉信号，输出RINEX数据格式的观测值O文件和导航电文N文件；然后处理O文件和N文件，可以计算得到各颗卫星的仰角、方位角、UTC时刻等信息，并将上述数据存成txt文档，并将上述数据存成txt文档；最后根据上述生成的txt文档文件和O文件共同的UTC时间信息，对每颗卫星筛选出txt文档中包含该颗卫星对应的UTC时间、仰角、方位角以及O文件包含的信噪比等级数据信息；

步骤2：数据筛选；

先选取低仰角范围的、方位角基本保持不变的，并且GPS信噪比等级数据具有明显变化的卫星及这段范围内的数据；然后将仰角正弦化，并按仰角正弦的从小到大顺序重新排列其对应的GPS信噪比等级数据；

步骤3：频谱分析并获取相位信息；

先对仰角的正弦和GPS信号强度数据进行等间隔均匀插值，以获取足够多的数据样点；然后用均匀采样信号频谱分析法分析频谱并获取主频率f；再对GPS信噪比等级数据进行曲线拟合，拟合得到的相位值即为所求；

步骤4：建立相位与湿度的函数关系并反演土壤湿度。

先按上述步骤处理得到多天的相位信息；然后绘制相位与实测湿度之间的变化关系，并根据两者的关系以及经验模型建立一种函数关系；最后处理某天的数据获得相位，将相位代入上述线性关系可反演出该天的土壤湿度。

作为优选，其特征在于：步骤1中步骤1中所述的处理O文件和N文件的过程，是根据频率L₁上C/A码所测定的伪距算出卫星与GPS接收机的空间位置，包括GPS接收机相对于卫星的仰角和方位角等信息。

作为优选，其特征在于：步骤2中所述的GPS信噪比数据的筛选，其筛选依据是选择具有两个完整周期的方波数据。

作为优选，其特征在于：步骤3中所述的均匀插值，采用的是最近邻插值法。

作为优选，其特征在于：步骤3中所述的频谱分析法，是快速傅里叶变换法FFT。

作为优选，其特征在于：步骤4的具体实现过程中，根据经验模型建立的是一种线性函数关系。

本发明首次提出了利用GPS干涉信号中的信噪比等级数据进行土壤湿度的反演。由于现有的土壤湿度反演方法中多数使用干涉信号中的原始信噪比数据进行，而由于许多CORS站并不提供原始信噪比数据，而只是在需要时才会接收以提供数据，因而本发明能够实现在历史数据的支持下，避免了对原始信噪比数据的依赖，通过使用信噪比等级数据就能够实现土壤湿度的反演。本方法能应用于农业或气象研究，本发明的实施将会扩大CORS站点数据的有效应用范围。

附图说明

图1是本发明实施例的方法流程图。

图2是本发明实施例的GPS卫星接收机的直达信号与反射信号几何关系图。

图3是本发明实施例的GPS信号的数据筛选。

图4是本发明实施例的GPS卫星信噪比等级数据的幅度谱图。

图5是本发明实施例的GPS卫星信噪比等级数据的相位与土壤湿度的线性关系。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明用到的GPS信噪比等级数据sn_rnx实际是接收机在接收的时候对信噪比数据的一个量化，满足以下关系：

sn_rnx＝MIN(MAX(INT(sn_raw/6,)1),9)

其中sn_raw为GPS信噪比数据(引用于RINEX格式的第三版本中第5章中的第7小节，标题为Signalstrengths)。在RINEX格式中，用1～9数字表示信噪比等级。其中1表示可能的最小信噪比等级，5表示良好的sn_raw的阈值，9表示可能的最大信噪比等级，0或空表示未知或未给出。

请见图1，本发明提供一种利用GPS信噪比等级数据进行土壤湿度估计的方法，包括以下步骤：

步骤1：GPS信噪比等级数据、仰角以及方位角的获取；

在GPS基准站附近设置土壤湿度观测站，天线架设高度约为2米，该天线为大地测量中常用的全向性GPS扼流环天线(choke-ring)，且具有近似方位角对称的增益和相位中心的模式。接收天线附近应视野开阔，地形平坦，地表应是有较少的矮草覆盖，四周没有高大建筑物和树木遮挡。

如图2所示，是GPS接收机的接收的直达信号与反射信号的几何关系。图中1表示接收机接收的GPS直达信号、2表示GPS信号射入土壤并反射到接收机的信号、3表示GPS卫星相对于地平面的仰角θ、4表示接收机相对于地平面的高度h、5表示反射信号的有效遥感深度Δh、6表示直达信号与反射信号的几何路径差ΔL，且ΔL＝2Hsinθ，H＝h+Δh。

接收机接收GPS卫星的直达信号与经由地面反射的反射信号合成的干涉信号，信号经接收机处理后，输出观测值O文件和导航电文N文件，其中干涉信号功率可以表示为：

$\begin{array}{c}P={\left|u\right|}^2={\left|u_{dir}+u_{ref}\right|}^2\\ ={\left|u_{dir}\right|}^2+{\left|u_{ref}\right|}^2+2\left|A\left(\mathrm{\θ}\right)\right|\cos \[\frac{4\mathrm{\π}H}{\mathrm{\λ}}\sin \mathrm{\θ}-\mathrm{\ψ}\]\end{array}$

u_dir表示直达信号，u_ref表示反射信号，A与直达信号和反射信号的幅度有关，H表示GPS接收机相对于反射平面的有效高度。ψ与天线特性以及反射面的特性有关，具体表现为与湿度有关。

由于O文件和N文件中包含伪距等信息，本发明根据频率L₁上C/A码所测定的伪距算出卫星与GPS接收机的空间位置，包括GPS接收机相对于卫星的仰角和方位角等信息，并将上述数据信息存成txt文档；然后利用上述生成的txt文件和O文件共同的UTC时间信息，通过计算机对每颗卫星筛选出txt文档中包含该颗卫星对应的UTC时间、仰角、方位角以及O文件包含的GPS信噪比等级数据信息；

步骤2：数据筛选；

由于GPS卫星相对于地面某一站点的周期约为24小时，因此为了准确性，选取数据依据如下：方位角要基本保持不变、仰角要处于在30度以下、并且GPS信噪比等级数据为具有两个完整周期的方波，根据该选择标准选择满足条件的卫星以及对应的仰角和GPS信噪比等级数据，并且记录选取的仰角θ_i，其卫星PRN以及信噪比等级序列E_i数据，图3为GPS01号星仰角、方位角与GPS信噪比等级数据的对应关系图，其中各个图的横坐标都是采样点数，7、8和9分别表示满足上述数据选择标准而所要选择的仰角、方位角以及GPS信噪比等级数据；然后对满足条件的每颗卫星PRN分别将其仰角正弦化sinθ_i，并按仰角正弦的从小到大顺序重新排列对应的GPS信噪比等级数据E_i，得到GPS信噪比等级数据随仰角的变化趋势序列；(x1_i,y1_i),i＝1,2,...,I₁,其中x1_i＝sin(θ_i),y1_i＝E_i。

步骤5：频谱分析以及提取相位信息；

从步骤2中的公式可知，GPS卫星干涉信号功率是仰角的正弦函数sinθ的时间函数，因θ是按时间均匀采样的，但转化成正弦后变成非均匀采样，因此先对仰角的正弦数据和GPS信噪比等级数据进行等间隔均匀采样，以获取足够多的数据样本点为(x_i,y_i),i＝1,2,...,I。然后本发明采用均匀采样频谱分析方法FFT(FastFouriertransform)分析频谱。

FFT是一种DFT的高效算法，称为快速傅立叶变换。由于样本点(x1_i,y1_i)都是离散点，用标准的离散数字信号序列表示为y(n)＝y(x_i),n＝i-1,n＝0,1,...,N-1,N＝I，一共有个N计时观测量。首先DFT的定义为：

$Y\left(k\right)=\underset{n=0}{\overset{N-1}{\Σ}}y\left(n\right)W_N^{kn},k=0,1,...,N-1,W_N=e^{-j\frac{2\mathrm{\π}}{N}};$

W_N为旋转因子，Y(k)为输入序列y(n)对应的N点对应频率点的相对幅度。利用W_N的对称性和周期性，即可按蝶形算法对上式进行一次次分解，算出频域结果。图4中上半部分的图是进行频谱分析后的结果，图4下半部分是对上半部分的部分放大，可以看出20Hz处为最大频谱密度点，我们将这个频率点选取出来作为干涉信号功率模型的主频率f。

求得主频率f后，为此我们将时域信号序列y(x_i)采用三参数的正弦振荡信号模型来表示为：

p＝acos(2πfx-ψ)+C,其中

其中f是正弦振荡的频率，其余3个参数为待求量。因参数ψ在p中以非线性方式出现，我们采用非线性最小二乘拟合法拟合该模型。具体步骤如下：

i.先计算出实际观测值y(x_i)在各点的残差平方和：

$Z(a,\mathrm{\ψ},C)=\underset{i=1}{\overset{N}{\Σ}}\left|{{\mathrm{\α}}_i}^2\right|=\underset{i=1}{\overset{N}{\Σ}}{\left(p\right(x_i,a,\mathrm{\ψ},C)-y(x_i\left)\right)}^2$

ii.然后Z(a,ψ,C)对参数a,ψ,C分别就偏导数并置为零，联立方程解出3个未知数a,ψ,C，ψ即为GPS信噪比等级数据y(x)的相位。

步骤6：建立相位与湿度的函数关系并反演土壤湿度；

为了建立相位与湿度的关系，需要通过GPS接收机接收数天的数据，并且测量对应天数的土壤湿度值。此外，待处理的GPS信噪比等级数据其对应的湿度值应变化较明显，这样可以排除一些因外界的干扰而造成的扰动。因此，本发明选取了湿度值变化明显的约三十天左右的数据，然后对待处理的数据进行上述相同的操作，并提取出ψ_j序列，j＝1,2,...,L，L为处理的总天数。同时将相位ψ_j与对应湿度值mv_j的变化趋势绘制在一张图中，然后将它们的线性关系拟合出来，建立经验模型。如图5所示，因为mv_j与ψ_j有明显的线性关系，建立模型ψ＝a₁mv+a₂，然后采用线性最小二乘法可以求出系数a₁，a₂：

i.先计算出残差平方和 $H(a_1,a_2)=\underset{j=1}{\overset{L}{\Σ}}\left|{{\mathrm{\β}}_j}^2\right|=\underset{j=1}{\overset{L}{\Σ}}{\left(\mathrm{\ψ}\right({\mathrm{mv}}_j,a_1,a_2)-{\mathrm{\ψ}}_j)}^2,$ 然后H(a₁,a₂)对参数a₁，a₂分别就偏导数并置为零其中m＝1,2。联立方程就出2个未知数a₁，a₂；

ii.通过求出的参数a_m确定相位与湿度的线性关系：

图5所示，为某地采样的土壤湿度值和相位的线性线关系，其中11为实际在一个月中采集的土壤湿度值mv及对应的相位值ψ，通过最小二乘法求得拟合的直线为12。

之后即可利用该线性关系，在算得某天的GPS信噪比等级数据的相位值之后可反演该天的土壤湿度。

本发明的特点是首次采用了CORS站的GPS信噪比等级数据进行土壤湿度的反演，避免了对原始信噪比数据的依赖。这样能够充分发挥CORS站历史数据的作用进行土壤湿度估计，利用率高。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (6)

1.一种利用GPS信噪比等级数据sn_rnx进行土壤湿度估计的方法，其特征在于，包括以下步骤：

步骤1：GPS信噪比等级数据、仰角以及方位角的获取；

利用地基CORS站的GPS接收机接收GPS干涉信号，输出RINEX(ReceiverIndependentExchangeFormat)数据格式的观测值O文件和导航电文N文件；然后处理O文件和N文件，可以计算得到各颗卫星的仰角、方位角、UTC时刻等信息，并将上述数据存成txt文档；最后根据上述生成的txt文档文件和O文件共同的UTC时间信息，对每颗卫星筛选出txt文档中包含该颗卫星对应的UTC时间、仰角、方位角以及Rinex格式的O文件包含的信噪比等级数据sn_rnx信息；

步骤2：数据筛选；

先选取低仰角范围的、方位角基本保持不变的，并且GPS信噪比等级数据具有明显变化的卫星及这段范围内的数据；然后将仰角正弦化，并按仰角正弦的从小到大顺序重新排列其对应的GPS信噪比等级数据；

步骤5：频谱分析并获取相位信息；

先对仰角的正弦和GPS信号强度数据进行等间隔均匀插值，以获取足够多的数据样点；然后用均匀采样信号频谱分析法分析频谱并获取主频率f；再对GPS信噪比等级数据进行曲线拟合，拟合得到的相位值即为所求；

步骤6：建立相位与湿度的函数关系并反演土壤湿度；

先按上述步骤处理得到多天的相位信息；然后绘制相位与实测湿度之间的变化关系，并根据两者的关系以及经验模型建立一种函数关系；最后处理某天的数据获得相位，将相位代入上述线性关系可反演出该天的土壤湿度。

2.根据权利要求1所述的利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于：采用的是sn_rnx数据，步骤1中所述的处理O文件和N文件得到仰角等信息，是根据频率L₁上C/A码所测定的伪距算出卫星与GPS接收机的空间位置，包括GPS接收机相对于卫星的仰角和方位角等信息。

3.根据权利要求1所述的利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于：步骤2中所述的GPS信噪比数据的筛选，其筛选依据是选择具有两个完整周期的方波数据。

4.根据权利要求1所述的利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于：步骤3中所述的均匀插值，采用的是最近邻插值法。

5.根据权利要求1所述的利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于：步骤3中所述的频谱分析法，是快速傅里叶变换法FFT。

6.根据权利要求1所述的利用GPS信噪比等级数据进行土壤湿度估计的方法，其特征在于：步骤6的具体实现过程中，根据经验模型建立的是一种线性函数关系。