基于广义延拓逼近的GNSS-IR土壤湿度反演方法

基于广义延拓逼近的gnss-ir土壤湿度反演方法
技术领域
1.本发明涉及全球卫星导航系统干涉反射遥感技术、生态环境监测及气象灾害防防治领域，具体来说，即一种基于广义延拓逼近模型的gnss-ir(global navigation satellite system-interferometric reflection，全球导航卫星系统干涉反射测量)土壤湿度反演方法。


背景技术：

2.全土壤湿度及地表土壤含水量是人类生产与生活中的重要环境信息，稳定、准确地预测土壤湿度变化对气象地质灾害和水循环研究有重大的意义。全球导航卫星系统干涉反射测量(global navigation satellite system interferometric reflectometry,gnss-ir)是一种新型的非合作式被动遥感技术，利用导航卫星的直射信号与反射信号的干涉信息来反演地表反射面的土壤湿度，很好地发挥导航卫星全天时天候、无需发射源、覆盖范围广、时空分辨率高等优势。目前国内外学者对gnss-ir土壤湿度监测进行了大量研究。
3.现有的研究在gnss系统、频段及卫星的选取，反射信号的分离与反演模型改进等方面均涌现了许多成果。然而仍存在许多需要改进的问题，首先在计算卫星高度角时，由于国际gnss服务中心(international gnss service,igs)发布的精密星历仅限00:00:00-23:45:00时间段，且采样间隔长达15分钟，想要得到任意时刻的卫星位置，需要对精密星历进行内插和外推处理。gnss-ir相关研究中通常使用teqc软件(translation,editing and quality checking，teqc)处理直接得到卫星位置坐标，其采用的拉格朗日插值法在高阶时容易产生“龙格”现象，且外推效果不理想。其次由于观测过程中卫星方位不断改变，其对应的镜面反射点位置也不断改变，因此在使用最小二乘正弦拟合时，直接用天线高度计算信噪比振荡频率是存在很大误差的。最后在现有的反演模型中，当卫星信噪比数据缺失或产生的相对延迟相位异常时，gnss-ir反演土壤湿度难以保持准确及稳定。综上所述，为了实现全面、准确的土壤湿度监测，gnss-ir反演技术仍有许多难题亟待突破。


技术实现要素：

4.针对上述现有技术的一些不足，本文采用工程科学中的广义延拓插值与外推模型对gnss-ir反演土壤湿度流程中的三个环节进行改进：利用广义延拓法对精密星历进行内插及外推，使得每个观测历元都能得到可靠的卫星位置坐标，可以计算得到相应的高精度卫星高度角；采用广义延拓插值逼近法和归一化lomb-scargle谱分析法(generalised lomb-scargle,gls)相结合的频率估计方法来优化相对延迟相位的计算；利用广义延拓法补充缺失的相对延迟相位数据，增加反演模型的稳定性。
5.为了实现上述目的，本发明提出一种基于广义延拓逼近模型的gnss-ir土壤湿度反演方法，本发明采取的具体技术方案具体包括以下步骤：
6.步骤一：采用连续运行的卫星定位服务参考站的测量型gnss接收机提供的观测数据，提取各卫星的信噪比数据及时间信息，以及测站周围对应时段的土壤湿度的实测数据。
7.步骤二：计算卫星高度角；
8.根据步骤一的观测数据，利用广义延拓插值逼近法对每天00:00:00～23:45:00精密星历数据进行内插处理，内插采用的逼近模型如下：
[0009][0010]
其中，a1，a2，a3为插值逼近函数的待定系数，xi为第i个历元的时间标志，yi为第i个历元的wgs-84地心地固坐标系下的卫星在轨位置坐标，(x
k-1
，y
k-1
)与(xk，yk)为插值点；
[0011]
利用广义延拓外推法对23:45:15～24:00:00时间段的精密星历进行外推处理，采用的外推模型如下：
[0012][0013]
以拉格朗日乘子法求得系数a1，a2，a3后，以最后m个数据的平均值作为插值点，n为数据总数，得到卫星位置坐标的广义延拓逼近函数，将卫星观测文件的采样间隔转化为对应时间戳代入，得到每历元的卫星在轨位置坐标，通过每历元的卫星在轨位置坐标计算卫星高度角与方位角；
[0014]
基于步骤一提取的接收机位置坐标与卫星在轨位置的相对关系计算卫星高度角与方位角，最后根据卫星编号和时间戳将每历元的snr、卫星高度角、方位角一一对应成组。
[0015]
步骤三：根据设置的卫星高度角范围，筛选有效的卫星信噪比数据，分离卫星信噪比反射分量；
[0016]
根据设置的卫星高度角范围5
°
～20
°
，筛选出有效的卫星信噪比数据；采用局部加权回归分离卫星反射信号时选取的窗口长度根据数据长度在0.5～0.8范围内变化，分离卫星信噪比反射分量。
[0017]
步骤四：对卫星信噪比反射分量进行重采样，用广义延拓逼近的频率估计方法计算信噪比序列的振荡频率，并利用非线性最小二乘拟合法估计特征参数，即相对延迟相位。
[0018]
基于步骤三提取的卫星反射分量信噪比的时间序列，根据对应历元的卫星高度角正弦值大小对信噪比数据进行重排序；采用基于广义延拓插值模型改进的归一化lomb-scargle谱分析法对重排序后的卫星反射分量信噪比序列进行频率估计，得到反射信号信噪比的振荡频率；
[0019]
先利用归一化lomb-scargle谱分析法谱分析生成信号关于测试频率的离散功率谱，在离散功率谱峰值点(x
p
,y
p
)的左侧或者右侧建立两个单元域，分别向左右拓展生成延拓域，相应的广义延拓插值逼近模型如下：
[0020][0021]
在解得待定系数(a
11
,a
12
,a
13
)和(a
21
,a
22
,a
23
)后根据广义延拓逼近函数在两单元域内找到对应的极大值点(x
e1
,y
max1
),(x
e2
,y
max2
)；再以两极大值点为边界点作第三个单元域，构建新的一元广义延拓模型如下：
[0022][0023]
在第三单元域δe3内得到极大值(x
e3
,y
max3
)；将三个极大值与离散功率谱峰值中的最大值对应的测试频率作为最终的频率估计结果f
max
,如下：
[0024][0025]
再将f
max
作为已知量通过非线性最小二乘正弦拟合法求出特征参数，采用的公式如下：
[0026][0027]
其中，snrr为卫星信噪比反射分量，ar为反射信号振幅，θ为卫星高度角，作为输入变量，为反射信号的相对延迟相位，即特征参数。
[0028]
步骤五：广义延拓插值法对测试样本中缺失的时间段的相对延迟相位进行预测。
[0029]
当某时间段内缺失卫星观测数据或无法计算出相对延迟相位时，利用广义延拓插值逼近法通过对相邻实测数据进行拟合和插值，得到缺损时间段的特征参数的预测值。
[0030]
采用广义延拓插值法对测试样本中缺失的时间段的相对延迟相位进行预测：
[0031]
先根据测试集缺失数据相对应的年积日确定需要插值或者外推的位置，在满足条件的两相对延迟相位间(doyn，doy
n 1
)建立单元域，并分别往外扩展一个数据点作为延拓域，使用的广义延拓逼近模型为：
[0032][0033]
xj表示(doyn，doy
n 1
)，yj表示的xj对应的信噪比；采用的拟合多项式u(x)＝a1 a2x a3x2，代表年积日为x时的相对延迟相位预测值；将目标时间段的年积日设为(doy1,
…
,doy
end
)。
[0034]
步骤六：建立基于随机森林的多星联合反演模型，并对待测时间段的土壤湿度进行反演。
[0035]
根据各卫星相对延迟相位与土壤湿度间的相关性大小，选取合适的卫星。将这些卫星的实际观测数据与对应的土壤湿度实测值作为训练集，使用随机森林回归建立多卫星组合的土壤湿度反演模型。将广义延拓法产生的相对延迟相位数据作为测试集，预测数据丢失时间段内的土壤湿度。
[0036]
与现有技术相比，本发明提出的基于广义延拓逼近模型的gnss-ir土壤湿度反演方法，利用广义延拓插值及外推法将历元间隔为15分钟的igs精密星历数据插值到15秒，可以准确高效地提取卫星高度角；利用广义延拓插值逼近法结合归一化gls谱分析法估计反射信号信噪比振荡频率，减小了直接使用天线高度计算频率带来的误差，有助于非线性最小二乘拟合法求解特征参数；可以应对gnss-ir土壤湿度反演中观测数据缺失的问题，缺失数据量在一定范围内时，仍能得到较准确的反演结果。
附图说明
[0037]
图1是卫星信号镜面反射几何关系图
[0038]
图2是本发明的流程图
[0039]
图3是本发明方法的具体实现步骤
[0040]
图4是广义延拓法内插精密星历示意图
[0041]
图5是广义延拓法外推精密星历示意图
[0042]
图6是广义延拓逼近的gls频率估计法的单元域划分示意图
[0043]
图7是非线性最小二乘拟合法相对延迟相位示意图
[0044]
图8是本发明方法的土壤湿度反演效果图
[0045]
图9是本发明方法的土壤湿度反演精度与数据缺失量的关系图
具体实施方式
[0046]
现结合附图对本发明具体实施案例做进一步描述。
[0047]
如图3所示，本发明利用测量型gnss接收机所采集的卫星信号反演地面近地表土壤湿度。本实例以本发明方案为前提，以板块边界观测计划(the plate boundary observation，pbo)p041测站2011年的卫星观测数据作为测试数据进行实例描述。
[0048]
图1为gnss-ir土壤湿度反演过程中的卫星信号镜面反射几何关系图，可以看出，gnss卫星发射的信号经由地面反射被gnss天线接收。当地面较平坦时，θ近似为卫星高度角；h为天线距离反射面的垂直高度。直射信号与反射信号的传播路径差d为：
[0049]
d＝2h sin θ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0050]
根据信号相位叠加的矢量关系，卫星直射与反射信号的载波相位差随着卫星高度角的变化而变化，卫星信号snr观测值在去除了直射分量后,信噪比反射分量snrr与高度角正弦值之间存在着一种正余弦关系：
[0051][0052]
其中，ar为反射分量的相对幅值；为反射分量的相对延迟相位。可见，ar与即为能表征地表土壤湿度的特征参数，参考前人的研究本发明中选择反射信号的相对延迟相位作为研究对象。
[0053]
图2所示，为本发明的实施流程，具体操作步骤如下：
[0054]
步骤一：导航卫星数据采集及预处理
[0055]
信噪比数据解析与卫星高度角的计算。采用p041测站所获取的观测文件：rinex2.11版本格式的观测o文件和igs提供的精密星历：sp3格式的过后精密星历。利用matlab程序解析并提取各卫星载波信噪比(l2频段)、卫星在轨位置坐标、接收机位置坐标及历元时间标签。利用广义延拓插值逼近法分别对00:00:00～23:45:00和23:45:15～24:00:00时间段的精密星历提供的卫星轨道位置进行内插和外推处理，处理后的精密星历如图4和图5所示。最后利用接收机与卫星的相对位置计算卫星高度角与方位角，得到各卫星每历元的信噪比及高度角信息。
[0056]
步骤二：反射信号信噪比的提取及重排序。
[0057]
基于步骤一得到的数据，以合理的预设高度角范围筛选满足镜面反射条件的有效卫星信号。使用局部加权回归法拟合各个卫星信噪比的变化趋势，并作为直射信号分量剔除，将残差序列作为卫星反射信号分量信噪比时间序列。最后将信噪比时间序列根据每历元的卫星高度角正弦值大小进行重排列，得到反射信号信噪比与卫星高度计正弦值之间的相干序列。
[0058]
步骤三：特征参数的计算及预测。
[0059]
利用广义延拓改进的归一化lomb-scargle谱分析法对重排列后的卫星反射信号snr序列进行频率估计，将其峰值f
max
作为信噪比振荡频率，所采用的一元广义延拓插值模型单元域划分规则如图6所示。采用非线性最小二乘拟合法进行正弦拟合，得到所需的特征参数，即相对延迟相位，拟合效果如图7所示。对于缺失卫星观测数据或无法计算出相对延迟相位的时间段，利用广义延拓插值逼近法通过对相邻实测数据进行拟合和插值，得到缺损的特征参数的预测值。
[0060]
步骤四：基于随机森林的多星联合土壤湿度反演
[0061]
将2011年第16～283天的相对延迟相位数据作为训练输入样本，相应的土壤湿度参考值作为训练输出样本；将2011年第284～334天的相对延迟相位作为预测输入样本。选取了5颗卫星的观测数据用于试验分析，分别是prn6、12、18、22、28号卫星，利用随机森林回归方法建立了各卫星的相对延迟相位与土壤湿度之间的反演模型，当设定的数据缺失天数为16天时，土壤湿度反演结果如图8所示，其中图8(a)为未丢失数据时的反演效果；图8(b)表示缺失数据的长度及分布情况；图8(c)为缺失数据补充后的反演效果；图8(d)表示土壤湿度反演结果的误差。
[0062]
如图9所示，本发明可以应对不同长度及分布的观测数据缺失情况，缺失数据量到达测试集的60％仍能保持0.85以上的相关系数。
[0063]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。