两种GNSS辅助WRF的实时PPP改善方法

本发明涉及gnss精密单点定位，尤其涉及两种gnss辅助wrf的实时ppp改善方法。

背景技术：

1、精密单点定位技术(precise point positioning,ppp)是指利用卫星实时或事后的精密轨道和钟差产品，对观测方程中的各项误差进行精确改正，同时采取合适的参数估计策略，在仅需一台gnss接收机的情况下，实现全球范围内精密绝对定位。但gnss信号穿过50km以下的中性大气层时，会被其中的中性大气干扰产生折射和弯曲，这种影响被称为对流层延迟。倾斜对流层延迟(slant tropospheric delay,std)可表示为ztd与映射函数(mapping functions,mf)的乘积加上大气水平梯度改正项与其映射函数乘积。ztd可分为天顶干延迟(zenith hydrostatic delay,zhd)和天顶湿延迟(zenith wet delay,zwd)。在传统ppp中，由于zhd较稳定且易于模型化，故常采用模型法进行改正，如saastamoinen模型、hopfiled模型等。随后，将模型计算的zhd作为ppp参数估计的先验值进行解算。这类模型往往需要gnss站点实测气象参数支持，但大多数gnss测站均未配备气象探测装置，因此气象参数往往难以获取。

2、数值天气预报(numerical weather prediction,nwp)模式通过求解一系列数学物理方程描述未来一段时间内的大气运动状态及天气变化过程，即可对温度、大气压强、相对湿度等气象变量进行精确预报。而数值天气预报结果的优劣在很大程度上取决于模式的初始湿度场，其分析质量会直接影响数值天气预报的准确性。为解决数值预报模式存在的初值问题，模式融合地面站点观测数据、卫星遥感数据等资料的数值同化技术应运而生，并成为辅助数值预报模式获取高精度、高分辨率气象参数的有效手段。nwp模型对于气象参数的获取不受时间和地点的约束，且同化gnss ztd可显著提高模式预报气象参数的精度，可利用其预报的气象参数解决gnss测站无法获取实测气象参数的缺陷，获取高精度的zhd。

3、但是，现有的技术难题是gnss ppp参数估计中计算先验zhd的站点气象参数难以获取以及计算的先验zhd精度较低，使用效果不够理想。

技术实现思路

1、本发明的目的是解决现有的gnss ppp参数估计中计算先验zhd的站点气象参数难以获取以及计算的先验zhd精度较低，使用效果不够理想的缺点，从而提出两种gnss辅助wrf的实时ppp改善方法。

2、为了实现上述目的，本发明采用了如下技术方案：

3、两种gnss辅助wrf的实时ppp改善方法，包括以下步骤：

4、s1、gnss ppp反演天顶对流层延迟：

5、gnss信号从卫星发射传播到接收机的途中会受到各类误差的干扰，这些误差会不同程度的影响定位精度，在充分顾及各项误差的情况下，gnss ppp观测方程具体表达如下：

6、

7、式中，表示载波相位观测值，i为卫星信号频率，r和s分别代表接收机和卫星系统，表示星地几何距离，tr,i为接收机钟差，表示卫星钟差，为电离层延迟，为对流层延迟，br,i和分别表示接收机和卫星端的伪距延迟，br,i和则分别为接收机端和卫星端的载波相位硬件延迟，λi为信号频率的波长，表示整周模糊度，为伪距的残差项，为载波相位的残差项；

8、上式中的对流层延迟即为std，其可以表示为天顶方向对流层延迟与映射函数的乘积加上大气水平梯度改正项与其映射函数乘积，ztd主要包括zhd和zwd两部分。因此，std具体表示如下：

9、std＝mfh×zhd+mfw×zwd+mfg·[gns cosα+gewsinα]+εt (2)

10、式中，mfh表示zhd对应的干映射函数，mfw表示zwd对应的湿映射函数，mfg表示水平梯度映射函数，gns和gew分别表示水平梯度改正的南北方向分量与东西方向分量，α为测站到卫星的方位角，εt表示对流层延迟残差；

11、在传统ppp解算中，利用gnss站点实测气压或经验气压模型计算zhd，并将其作为先验值与其他先验参数一同代入线性化后的观测方程中，利用卡尔曼滤波方法对ppp函数模型进行解算，获取ztd估值。

12、s2：gnss辅助wrf的三维变分同化：

13、数值同化技术将gnss ppp反演的ztd实际观测值及误差信息作用于数值天气预报模式，以期获得有所改善的大气状态估计值。以wrf模式的三维变分同化系统进行同化可理解为观测场与分析场、分析场与背景误差场之间的二次泛函极小化问题，其目标函数的表达式如下：

14、

15、式中，x为待估计的大气状态向量；xb为大气状态的先验信息；yo为实际观测值向量；h为观测算子；h(x)为投影函数；b为xb的误差协方差矩阵；r为实际观测值的误差协方差矩阵；目标函数式亦可表示为：

16、j(x)＝jb(x)+jo(x) (4)

17、式中，jb(x)表示分析量x偏离背景场xb的程度；jo(x)表示由投影函数计算的分析量x偏离实际观测值的程度；目标函数j(x)的极小值是通过逐步迭代极小化方法得到的，最优分析场为代价函数极小化时的解。

18、基于wrf-da三维变分同化系统对gnss ztd进行同化时，将gnss ztd观测资料加入式jo(x)中，即：

19、

20、式中，hztd为gnss ztd的观测算子，r为gnss ztd的观测误差协方差矩阵。

21、s3：基于wrf预报的gnss站点气象参数获取：

22、s3.1：分层水蒸气分压获取：

23、wrf模式对gnss ztd观测资料进行三维变分同化有效提高模式的预报能力，可以输出气压、温度以及水汽混合比等高精度的气象预报参数；不同层的水蒸气分压e(hpa)可由水蒸气混合比q(kg·kg-1)转换得到：

24、

25、式中，p为空气总压(hpa)，即空气基态压力与扰动压力之和，其计算公式如下：

26、p＝pbase+pper (7)

27、式中，pbase、pper分别为空气基态压力和扰动压力(hpa)；

28、s3.2：分层温度获取：

29、将分层的空气总电位温度转换为空气温度可用如下公式表示：

30、

31、式中，tθ为模式输出的潜在温度(k)；p为大气压强(hpa)；rd为干空气气体常数，其值为287.053j·k-1·kg-1；cp为恒压下的比热容；rd和cp的比值为泊松常数，rd/cp＝2/7；

32、s3.3：气压获取：

33、gnss站点的气压可根据wrf模式输出的格网点气压参数插值得到；由于gnss站点坐标的高度为大地高，wrf模式输出的高度为位势高度，为使两者高度统一，需将gnss站点的大地高改正到位势高；因此，采用压高公式对wrf模式获取的气压进行修正，得到gnss站点的表面气压值，其计算公式如下：

34、

35、式中，hs、hg分别为gnss站点高度和wrf模式格网点高度(m)；ps(hs)、ps(hg)分别表示gnss站点高度和wrf模式格网点高度对应的表面气压(hpa)；gm为重力加速度，取值为9.80665m/s2，其计算公式如下：

36、

37、式中，为gnss站点纬度。

38、s4：高精度先验对流层干延迟zhd获取：

39、依据wrf模式提供的不同气象参数，可通过以下两种方式计算zhd。

40、s4.1：基于wrf模式预报气压的zhd计算：

41、利用步骤s3联合saastamoinen模型可获取高精度的zhd值，计算公式如下：

42、

43、式中，ps为基于wrf模式预报的气象参数计算的gnss站点气压(hpa)，h为测站高度(m)，为gnss站点纬度；

44、将公式(11)计算的zhd值取代传统ppp中的先验zhd值作为卡尔曼滤波估计的初值，进行高精度的解算；

45、s4.2：基于wrf模式预报气象参数的zhd积分获取方法：

46、利用数值天气预报模式资料反演对流层延迟可采用分层积分法，依据探空数据计算对流层延迟的方法，基于wrf模式输出参数分层计算zhd，在5hpa以下部分采用直接积分法计算，其具体公式如下：

47、

48、式中，dp为天顶干延迟不同分层的大气压强(hpa)；plevel为分层气压积分限，其值为5hpa；ps为gnss站点高度的大气压(hpa)；nd为大气干折射率，其表达式为：

49、

50、式中，pd为干空气大气压(hpa)；e为水蒸气分压(hpa)；t为大气温度(k)；

51、对于5hpa以上部分的zhd分量，以模式顶层的气象参数作为输入参数，采用saastamoinen模型进行计算，其具体公式如下：

52、

53、式中，ptop为模式顶层高度大气压(hpa)；htop为模式顶层高度(m)；

54、根据以上两部分计算的zhd分量，总的zhd可由下式计算得到：

55、zhdgrid＝zhdlevel+zhdtop (15)

56、而由于gnss站点一般不与wrf模型的网格点重合，因此需采用双线性插值(水平方向)和高度归算(垂直方向)获得gnss站点的zhd。

57、同理，将采用插值和高度归算获得的gnss站点zhd值作为卡尔曼滤波估计的初值，进行ppp解算。

58、s5：基于两种先验zhd的实时ppp解算：

59、利用wrf模式预报气压、水汽混合比等气象参数，得出两种不同的zhd值，考虑到其对ppp结果的改善较为精确，故将其作为载波相位观测方程滤波的zhd先验值进行解算，计算过程如下：

60、s5.1：gnss ppp载波相位观测方程线性化

61、

62、s5.2：由观测方程构建状态方程

63、xk+1＝φxk+wk (17)

64、s5.3：滤波过程的观测方程矩阵

65、lk+1＝ak+1·xk+1+vk+1 (18)

66、上述过程中，为每次滤波的观测方程矩阵，为滤波过程中的系数矩阵，为待估参数，为测量噪声；

67、在上述过程中给定先验zhd与其它先验参数，随后利用卡尔曼滤波进行参数估计，将得到改善后的ppp定位结果，该方法能够明显提高定位精度及加快收敛速度。

68、本发明中，所述的两种gnss辅助wrf的实时ppp改善方法，首先，依据传统ppp技术进行ztd反演，并将其同化到wrf模式中预报相关气象参数；其次，根据预报产品分层格网点上的温度、压强及水汽混合比等气象参数，基于经验模型或积分方法计算zhd。最后，将估算的两种zhd分别代入ppp的观测方程中，达到提高ppp定位精度和加快收敛速度的目的。

69、本发明通过利用数值天气预报模式预报的气象参数提出两种方法估计高精度的zhd初值，以达到提高ppp定位精度及收敛速度目的，值得推广使用。