CN116466376A - 数值预报模式辅助的实时ppp改善方法 - Google Patents
数值预报模式辅助的实时ppp改善方法 Download PDFInfo
- Publication number
- CN116466376A CN116466376A CN202310249580.7A CN202310249580A CN116466376A CN 116466376 A CN116466376 A CN 116466376A CN 202310249580 A CN202310249580 A CN 202310249580A CN 116466376 A CN116466376 A CN 116466376A
- Authority
- CN
- China
- Prior art keywords
- gnss
- wrf
- delay
- ppp
- atmospheric
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 239000005436 troposphere Substances 0.000 claims abstract description 18
- 230000008569 process Effects 0.000 claims abstract description 17
- 238000004364 calculation method Methods 0.000 claims abstract description 14
- 238000012937 correction Methods 0.000 claims abstract description 12
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000005516 engineering process Methods 0.000 claims abstract description 5
- 230000002706 hydrostatic effect Effects 0.000 claims description 21
- 230000014509 gene expression Effects 0.000 claims description 18
- 239000011159 matrix material Substances 0.000 claims description 18
- 230000003068 static effect Effects 0.000 claims description 18
- 238000011160 research Methods 0.000 claims description 16
- 238000004458 analytical method Methods 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 15
- 238000013507 mapping Methods 0.000 claims description 15
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 10
- 238000011835 investigation Methods 0.000 claims description 6
- 230000001902 propagating effect Effects 0.000 claims description 6
- 230000001934 delay Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 239000005433 ionosphere Substances 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000008054 signal transmission Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 abstract description 4
- 230000007547 defect Effects 0.000 description 1
- 230000008570 general process Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01W—METEOROLOGY
- G01W1/00—Meteorology
- G01W1/10—Devices for predicting weather conditions
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Environmental & Geological Engineering (AREA)
- General Physics & Mathematics (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Life Sciences & Earth Sciences (AREA)
- Atmospheric Sciences (AREA)
- Biodiversity & Conservation Biology (AREA)
- Ecology (AREA)
- Environmental Sciences (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开了数值预报模式辅助的实时PPP改善方法，主要涉及GNSS气象学技术领域，针对现有的GNSS数据处理过程中，经验模型改正和未知参数估计往往存在一定偏差，导致PPP定位精度较低，使用效果不够理想的问题，现提出如下方案，包括以下步骤：S1：GNSS PPP反演天顶对流层延迟ZTD；S2：GNSS ZTD的三维变分同化；S3：GNSS辅助WRF模式的气象参数预测；S4：GNSS卫星信号斜路径对流层延迟估算；S5：WRF辅助的实时PPP改善方法。本发明针对PPP定位观测方程中包含的对流层延迟作为待估参数的现状，提出一种基于WRF同化ZTD预报的气象参数计算高精度STD方法，并将其直接改正到实时PPP定位观测方程中，以减少观测方程中的未知参数，提升PPP技术的定位精度及其收敛速度。
Description
技术领域
本发明主要涉及GNSS气象学技术领域，尤其涉及数值预报模式辅助的实时PPP改善方法。
背景技术
对流层延迟由对流层中的干空气成分和水汽造成，且与GNSS站点高度参数有较强的相关性，是影响精密单点定位(Precise Point Positioning,PPP)精度的主要因素之一,也是GNSS气象学中重要的研究对象。精确的斜路径对流层延迟(Slant TroposphericDelay,STD)估计有助于改善PPP的定位精度和收敛性。传统PPP定位中的STD可以表示为流体静力学延迟(Zenith Hydrostatic Delay,ZHD)和非静力学延迟(Zenith Wet Delay,ZWD)与其相应的映射函数(Mapping Functions,MF)的乘积，再加上大气水平梯度改正项乘以其映射函数。由于ZHD较稳定且易于模型化，故在一般处理中常采用经验模型进行改正，如Saastamoinen模型、Hopfiled模型等；而ZWD不确定性较强，难以建立精确模型改正，因此利用参数估计方法确定其大小。
但是，在传统GNSS数据处理中，经验模型改正和未知参数估计往往存在一定偏差，导致PPP定位精度较低，使用效果不够理想。
发明内容
本发明的目的是解决现有的GNSS数据处理过程中，经验模型改正和未知参数估计往往存在一定偏差，导致PPP定位精度较低，使用效果不够理想的缺点，因而提出数值预报模式辅助的实时PPP改善方法。
为了实现上述目的，本发明采用了如下技术方案：
数值预报模式辅助的实时PPP改善方法，包括以下步骤：
S1：GNSS PPP反演天顶对流层延迟ZTD：
GNSS接收机接收到的信号主要包括载波相位观测值与测码伪距观测值/>GNSS数据处理主要利用两种观测值组成非差观测方程；GNSS非差观测方程可以很好地保留接收机钟差、卫星钟差、电离层延迟、对流层延迟等误差信息，具体公式如下：
式中，表示站星几何距离，fi为卫星信号频率，i表示信号频率标志，r和s分别代表接收机和卫星系统，/>为测码伪距观测值，/>表示载波观测值，vr,i、kr,i和Kr,i分别表示接收机端的钟差、伪距延迟以及载波相位硬件延迟，属于接收机端误差；/>以及/>分别表示卫星端的钟差、伪距延迟以及载波相位硬件延迟，属于卫星端误差；/>和/>分别表示电离层延迟和对流层延迟，属于信号传输误差；λi为载波信号频率的波长，/>表示整周模糊度；/>和/>分别表示伪距的残差项载波相位的残差项；
公式(1)中的STD，可以表示为流体静力学延迟与非静力学延迟与其相应的映射函数的乘积，还需要顾及大气各项异性的情况，因此，STD的计算公式具体表示如下：
式中，ε为卫星高度角，为测站到卫星的方位角，Mh(ε)为天顶流体静力学延迟ZHD对应的映射函数，Mw(ε)为天顶非流体静力学延迟ZWD对应的映射函数，Mg(ε)为水平梯度映射函数，Gns和Gew分别表示水平梯度改正的南北方向分量与东西方向分量，γt表示对流层延迟残差；
在GNSS数据处理过程中，首先利用GNSS站点实测气压解算ZHD，其次将其作为先验值与其他先验参数一同代入线性化后的观测方程中，最后利用卡尔曼滤波方法对PPP函数模型进行解算，获取ZTD估值；
S2：GNSS ZTD的三维变分同化：
数值同化技术将GNSS ZTD及其误差信息作用于WRF模式，可获得改善后的大气状态估计值。基于WRF模式的三维变分同化，可理解为观测场与分析场、分析场与背景误差场之间的二次泛函极小化问题，其目标函数的表达式如下：
式中，x为待估计的大气状态向量，xb为大气状态的先验信息(背景场变量)，yo为实际观测值向量，H为观测算子；H(x)为关联背景场与实际观测值的投影函数，B为背景场xb的误差协方差矩阵，R为实际观测值的误差协方差矩阵，目标函数式亦可表示为：
J(x)＝Jb(x)+Jo(x) (4)
式中，Jb(x)表示分析量x偏离背景场xb的程度，Jo(x)表示由投影函数计算的分析量x偏离实际观测值的程度；目标函数J(x)的极小值是通过逐步迭代极小化方法得到的，最优分析场为代价函数极小化时的解，即待估计的真实大气状态x可表示为：
x＝xb+[B-1+HTR-1H]-1·HR-1[yo-H(xb)] (5)
基于WRF-DA三维变分同化系统对GNSS ZTD进行同化时，将GNSS ZTD观测资料加入式Jo(x)中，即：
式中，HZTD为GNSS ZTD的观测算子，R为GNSS ZTD的观测误差协方差矩阵；
S3：GNSS辅助WRF模式的气象参数预测：
WRF模式同化GNSS ZTD可提高预报气象参数的精度，其输出参数包含与计算对流层延迟密切相关的温度、气压、相对湿度、水汽混合比等；
大气压强可由WRF模式输出的空气基准态气压和扰动气压计算得到，其表达式如下：
P＝Pbase+Pper (7)
式中，Pbase、Pper分别表示模式输出的空气基准态气压和扰动气压(hPa)；
基于计算的大气压强P(hPa)，水蒸气分压e(hPa)可由模式输出的水蒸气混合比Q(kg·kg-1)转换得到：
大气温度可通过转换分层的空气总电位温度得到，其表达式为：
式中，Tθ为模式输出的潜在温度(K)，Rd为干空气气体常数，其值为287.053J·K-1·kg-1，cp为恒压下的比热容，Rd和cp的比值为泊松常数，Rd/cp＝2/7；
S4：GNSS卫星信号斜路径对流层延迟估算：
GNSS卫星信号的斜路径对流层延迟可基于卫星信号高度角、方位角及其在WRF三维网格中的截距对大气折射率进行积分获得，即获得GNSS卫星信号的SHD和SWD，进而求得STD，具体估算步骤如下：
S4.1：建立局部空间直角坐标系：
由于WRF模式所划分的网格在Z方向的分辨率一般为500m～1km，故可以忽略由垂线偏差所造成的影响，对椭球上的法线与垂线不加以严格区分；以研究区域左下角(西南角)的WRF格网点为坐标原点，北方向为X方向，东方向为Y方向，位势高方向为Z方向，建立局部空间直角坐标系；在此坐标系的基础上，计算出GNSS站点的空间直角坐标，并根据WRF网格空间分辨率、垂直分层计算每个WRF格网点的空间直角坐标；
在数值预报模式中，由于不同研究区域的下垫面情况不同，即使在一个很小的区域，地面的海拔也不一样，所以不适于按照位势高度来进行垂直分层；因此，WRF模式采用气压追随坐标eta在垂直方向上进行分层表示，取值范围从1到0，其表达式为：
eta＝(Plevel-Ptop)/(Pbot-Ptop) (10)
式中，Plevel为模式每一分层的气压(hPa)；Ptop表示模式顶层气压(hPa)；Pbot表示地表气压(hPa)；从eta坐标的表达式可以看出，地表的eta值为1；顶层的eta值为0；
而eta坐标对应的离地位势高度，可经过以下公式进行转换得到，即：
gmp＝(PH+PHB)/9.81-HGT (11)
式中，PH、PHB分别为模式输出的扰动位势高度和基准态位势高度；HGT为研究区域的地面海拔高度；位势高度和几何高度在数值上比较接近，则根据eta坐标对应的离地位势高度可以计算出WRF不同分层之间的高度；
S4.2：计算GNSS信号射线在WRF三维网格中的截距：
WRF模式研究区域内所有与X轴垂直的平面称为X型面，所有与Y轴垂直的平面称为Y型面，与Z轴垂直的平面称为Z型面；根据网格划分和垂直分层情况，WRF模式研究区域被离散化成若干X、Y和Z型面；GNSS信号射线穿过研究区域，与每一个型面有且只有一个交点；因此，计算GNSS信号射线与WRF网格交点即求解信号射线与所有X、Y、Z型面的交点坐标；若GNSS站点坐标为G，(x0,y0,z0)射线方向向量为S(m,n,p)，则信号射线与X型面的交点(xxi,yxi,zxi)可表示为：
GNSS信号射线与Y型面的交点(xyi,yyi,zyi)可表示为：
GNSS信号射线与Z型面的交点(xzi,yzi,zzi)可表示为：
则任意一条GNSS信号射线在网格(i,j,k)内的截距可表示为
式中，(xq,yq,zq)，(xq-1,yq-1,zq-1)表示GNSS信号射线与WRF网格(i,j,k)的两个交点坐标；
S4.3：STD估计：
(a)SHD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对网格内的大气干折射率进行积分，得到不同网格内的斜对流层静力学延迟，以WRF模式输出的气象参数直接积分计算大气层的静力学延迟，其具体公式如下：
SHD＝10-6·∫Nddl (16)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nd为WRF模式三维网格中的大气干折射率，其计算公式为：
式中，Pd为干大气压(hPa)，e为水蒸气分压(hPa)，T为大气温度(K)；
根据WRF模式垂直分层情况，将研究区域离散化为若干个独立的三维单元网格，并以每个网格中心点的大气干折射率代表三维网格的大气干折射率，则GNSS信号在斜路径上的对流层静力学延迟可表达为以下形式：
SHD＝∑ijk(dijk·Nd(ijk)) (18)
式中，不考虑地球曲率时，dijk表示GNSS信号射线在(i,j,k)网格的截距；考虑地球曲率时，dijk表示(i,j,k)网格的内插系数；Nd(ijk)表示(i,j,k)网格的大气干折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气干折射率插值得到；
(b)SWD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对不同网格内的大气湿折射率进行积分，同理可得其计算公式如下：
SWD＝10-6·∫Nwdl (19)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nw为WRF模式三维网格中的大气湿折射率，其计算公式为：
式中，Rw＝461.495J·K-1·kg-1为水蒸气的气体常数，T为大气温度(k)，k1＝77.6890K/hPa、k2＝71.2952K/hPa、k3＝375463K2/hPa为折射率常数；
同理可得，GNSS信号在斜路径上的非静力学延迟可表达为以下形式：
SWD＝∑ijk(dijk·Nw(ijk)) (21)
式中，Nw(ijk)表示(i,j,k)网格的大气湿折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气湿折射率插值得到；
最终，斜路径对流层延迟可由基于WRF模式输出气象参数计算的斜对流层静力学延迟和非静力学延迟求和得到，其具体表达式如下：
STD＝SHD+SWD (22)
式中，SHD、SWD分别为斜路径对流层静力学延迟和非静力学延迟；
S5：WRF辅助的实时PPP改善方法：
将基于WRF模式预报气象参数估算的STD加入到PPP定位观测方程中，对传统PPP定位方法进行改进，因此，顾及WRF模式预报气象参数改善STD的PPP载波相位观测方程表达如下：
式中，表示顾及STD项的载波相位观测值，/>为经过改正后的载波相位残差，其余各项不变；
通过对上式线性化，并联合多测站多卫星载波相位观测值构建GNSS PPP的函数模型和误差模型，利用卡尔曼滤波方法对PPP函数模型进行实时解算；
利用卡尔曼滤波方法对PPP载波相位观测方程进行解算时，通过将基于WRF模式预报气象参数估算的STD直接应用到载波相位观测方程中，减少了待估参数的个数；因此，滤波过程的各项权阵、滤波增益矩阵、协方差阵以及状态值矩阵都随之减少一阶，使得PPP解算过程更为简单高效，同时可以得到更高精度的PPP定位结果，并加快其收敛速度。
本发明中，所述的数值预报模式辅助的实时PPP改善方法，数值天气预报(Numerical Weather Prediction,NWP)模式可对未来一段时间内的大气温度、压强、相对湿度等气象参数进行预报，其获取气象参数不受时间和地点的约束，且同化ZTD可对气象参数进行更高精度的预报。而采用NWP模型预报气象参数反演的STD精度对最终测站坐标的精度和收敛速度有重要影响，利用数值预报模式反演STD可以减少对流层对GNSS卫星信号传播的影响，有效改善采用传统方法进行参数估计所产生的误差；
本发明针对PPP定位观测方程中包含的对流层延迟作为待估参数的现状，提出一种基于WRF同化ZTD预报的气象参数计算高精度STD方法，并将其直接改正到实时PPP定位观测方程中，以减少观测方程中的未知参数，进一步提升PPP技术的定位精度及其收敛速度，使用效果好。
附图说明
图1为本发明提出的数值预报模式辅助的实时PPP改善方法的流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。
参照图1，本方案提供的一种实施例：数值预报模式辅助的实时PPP改善方法，包括以下步骤：
S1：GNSS PPP反演天顶对流层延迟ZTD：
GNSS接收机接收到的信号主要包括载波相位观测值与测码伪距观测值/>GNSS数据处理主要利用两种观测值组成非差观测方程；GNSS非差观测方程可以很好地保留接收机钟差、卫星钟差、电离层延迟、对流层延迟等误差信息，具体公式如下：
式中，表示站星几何距离，fi为卫星信号频率，i表示信号频率标志，r和s分别代表接收机和卫星系统，/>为测码伪距观测值，/>表示载波观测值，vr,i、kr,i和Kr,i分别表示接收机端的钟差、伪距延迟以及载波相位硬件延迟，属于接收机端误差；/>以及/>分别表示卫星端的钟差、伪距延迟以及载波相位硬件延迟，属于卫星端误差；/>和/>分别表示电离层延迟和对流层延迟，属于信号传输误差；λi为载波信号频率的波长，/>表示整周模糊度；/>和/>分别表示伪距的残差项载波相位的残差项；
公式(1)中的STD，可以表示为流体静力学延迟与非静力学延迟与其相应的映射函数的乘积，还需要顾及大气各项异性的情况，因此，STD的计算公式具体表示如下：
式中，ε为卫星高度角，为测站到卫星的方位角，Mh(ε)为天顶流体静力学延迟ZHD对应的映射函数，Mw(ε)为天顶非流体静力学延迟ZWD对应的映射函数，Mg(ε)为水平梯度映射函数，Gns和Gew分别表示水平梯度改正的南北方向分量与东西方向分量，γt表示对流层延迟残差；
在GNSS数据处理过程中，首先利用GNSS站点实测气压解算ZHD，其次将其作为先验值与其他先验参数一同代入线性化后的观测方程中，最后利用卡尔曼滤波方法对PPP函数模型进行解算，获取ZTD估值；
S2：GNSS ZTD的三维变分同化：
数值同化技术将GNSS ZTD及其误差信息作用于WRF模式，可获得改善后的大气状态估计值。基于WRF模式的三维变分同化，可理解为观测场与分析场、分析场与背景误差场之间的二次泛函极小化问题，其目标函数的表达式如下：
式中，x为待估计的大气状态向量，xb为大气状态的先验信息(背景场变量)，yo为实际观测值向量，H为观测算子；H(x)为关联背景场与实际观测值的投影函数，B为背景场xb的误差协方差矩阵，R为实际观测值的误差协方差矩阵，目标函数式亦可表示为：
J(x)＝Jb(x)+Jo(x) (4)
式中，Jb(x)表示分析量x偏离背景场xb的程度，Jo(x)表示由投影函数计算的分析量x偏离实际观测值的程度；目标函数J(x)的极小值是通过逐步迭代极小化方法得到的，最优分析场为代价函数极小化时的解，即待估计的真实大气状态x可表示为：
x＝xb+[B-1+HTR-1H]-1·HR-1[yo-H(xb)] (5)
基于WRF-DA三维变分同化系统对GNSS ZTD进行同化时，将GNSS ZTD观测资料加入式Jo(x)中，即：
式中，HZTD为GNSS ZTD的观测算子，R为GNSS ZTD的观测误差协方差矩阵；
S3：GNSS辅助WRF模式的气象参数预测：
WRF模式同化GNSS ZTD可提高预报气象参数的精度，其输出参数包含与计算对流层延迟密切相关的温度、气压、相对湿度、水汽混合比等；
大气压强可由WRF模式输出的空气基准态气压和扰动气压计算得到，其表达式如下：
P＝Pbase+Pper (7)
式中，Pbase、Pper分别表示模式输出的空气基准态气压和扰动气压(hPa)；
基于计算的大气压强P(hPa)，水蒸气分压e(hPa)可由模式输出的水蒸气混合比Q(kg·kg-1)转换得到：
大气温度可通过转换分层的空气总电位温度得到，其表达式为：
式中，Tθ为模式输出的潜在温度(K)，Rd为干空气气体常数，其值为287.053J·K-1·kg-1，cp为恒压下的比热容，Rd和cp的比值为泊松常数，Rd/cp＝2/7；
S4：GNSS卫星信号斜路径对流层延迟估算：
GNSS卫星信号的斜路径对流层延迟可基于卫星信号高度角、方位角及其在WRF三维网格中的截距对大气折射率进行积分获得，即获得GNSS卫星信号的SHD和SWD，进而求得STD，具体估算步骤如下：
S4.1：建立局部空间直角坐标系：
由于WRF模式所划分的网格在Z方向的分辨率一般为500m～1km，故可以忽略由垂线偏差所造成的影响，对椭球上的法线与垂线不加以严格区分；以研究区域左下角(西南角)的WRF格网点为坐标原点，北方向为X方向，东方向为Y方向，位势高方向为Z方向，建立局部空间直角坐标系；在此坐标系的基础上，计算出GNSS站点的空间直角坐标，并根据WRF网格空间分辨率、垂直分层计算每个WRF格网点的空间直角坐标；
在数值预报模式中，由于不同研究区域的下垫面情况不同，即使在一个很小的区域，地面的海拔也不一样，所以不适于按照位势高度来进行垂直分层；因此，WRF模式采用气压追随坐标eta在垂直方向上进行分层表示，取值范围从1到0，其表达式为：
eta＝(Plevel-Ptop)/(Pbot-Ptop) (10)
式中，Plevel为模式每一分层的气压(hPa)；Ptop表示模式顶层气压(hPa)；Pbot表示地表气压(hPa)；从eta坐标的表达式可以看出，地表的eta值为1；顶层的eta值为0；
而eta坐标对应的离地位势高度，可经过以下公式进行转换得到，即：
gmp＝(PH+PHB)/9.81-HGT (11)
式中，PH、PHB分别为模式输出的扰动位势高度和基准态位势高度；HGT为研究区域的地面海拔高度；位势高度和几何高度在数值上比较接近，则根据eta坐标对应的离地位势高度可以计算出WRF不同分层之间的高度；
S4.2：计算GNSS信号射线在WRF三维网格中的截距：
WRF模式研究区域内所有与X轴垂直的平面称为X型面，所有与Y轴垂直的平面称为Y型面，与Z轴垂直的平面称为Z型面；根据网格划分和垂直分层情况，WRF模式研究区域被离散化成若干X、Y和Z型面；GNSS信号射线穿过研究区域，与每一个型面有且只有一个交点；因此，计算GNSS信号射线与WRF网格交点即求解信号射线与所有X、Y、Z型面的交点坐标；若GNSS站点坐标为G，(x0,y0,z0)射线方向向量为S(m,n,p)，则信号射线与X型面的交点(xxi,yxi,zxi)可表示为：
GNSS信号射线与Y型面的交点(xyi,yyi,zyi)可表示为：
GNSS信号射线与Z型面的交点(xzi,yzi,zzi)可表示为：
则任意一条GNSS信号射线在网格(i,j,k)内的截距可表示为
式中，(xq,yq,zq)，(xq-1,yq-1,zq-1)表示GNSS信号射线与WRF网格(i,j,k)的两个交点坐标；
S4.3：STD估计：
(a)SHD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对网格内的大气干折射率进行积分，得到不同网格内的斜对流层静力学延迟，以WRF模式输出的气象参数直接积分计算大气层的静力学延迟，其具体公式如下：
SHD＝10-6·∫Nddl (16)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nd为WRF模式三维网格中的大气干折射率，其计算公式为：
式中，Pd为干大气压(hPa)，e为水蒸气分压(hPa)，T为大气温度(K)；
根据WRF模式垂直分层情况，将研究区域离散化为若干个独立的三维单元网格，并以每个网格中心点的大气干折射率代表三维网格的大气干折射率，则GNSS信号在斜路径上的对流层静力学延迟可表达为以下形式：
SHD＝∑ijk(dijk·Nd(ijk)) (18)
式中，不考虑地球曲率时，dijk表示GNSS信号射线在(i,j,k)网格的截距；考虑地球曲率时，dijk表示(i,j,k)网格的内插系数；Nd(ijk)表示(i,j,k)网格的大气干折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气干折射率插值得到；
(b)SWD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对网格内的大气湿折射率进行积分，同理可得其计算公式如下：
SWD＝10-6·∫Nwdl (19)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nw为WRF模式三维网格中的大气湿折射率，其计算公式为：
式中，Rw＝461.495J·K-1·kg-1为水蒸气的气体常数，T为大气温度(k)，k1＝77.6890K/hPa、k2＝71.2952K/hPa、k3＝375463K2/hPa为折射率常数；
同理可得，GNSS信号在斜路径上的非静力学延迟可表达为以下形式：
SWD＝∑ijk(dijk·Nw(ijk)) (21)
式中，Nw(ijk)表示(i,j,k)网格的大气湿折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气湿折射率插值得到；
最终，斜路径对流层延迟可由基于WRF模式输出气象参数计算的斜对流层静力学延迟和非静力学延迟求和得到，其具体表达式如下：
STD＝SHD+SWD (22)
式中，SHD、SWD分别为斜路径对流层静力学延迟和非静力学延迟；
S5：WRF辅助的实时PPP改善方法：
将基于WRF模式预报气象参数估算的STD加入到PPP定位观测方程中，对传统PPP定位方法进行改进，因此，顾及WRF模式预报气象参数改善STD的PPP载波相位观测方程表达如下：
式中，表示顾及STD项的载波相位观测值，/>为经过改正后的载波相位残差，其余各项不变；
通过对上式线性化，并联合多测站多卫星载波相位观测值构建GNSS PPP的函数模型和误差模型，利用卡尔曼滤波方法对PPP函数模型进行实时解算；
利用卡尔曼滤波方法对PPP载波相位观测方程进行解算时，通过将基于WRF模式预报气象参数估算的STD直接应用到载波相位观测方程中，减少了待估参数的个数；因此，滤波过程的各项权阵、滤波增益矩阵、协方差阵以及状态值矩阵都随之减少一阶，使得PPP解算过程更为简单高效，同时可以得到更高精度的PPP定位结果，并加快其收敛速度
以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。
Claims (1)
1.数值预报模式辅助的实时PPP改善方法，其特征在于，包括以下步骤：
S1：GNSS PPP反演天顶对流层延迟ZTD：
GNSS接收机接收到的信号主要包括载波相位观测值与测码伪距观测值/>GNSS数据处理主要利用两种观测值组成非差观测方程；GNSS非差观测方程可以很好地保留接收机钟差、卫星钟差、电离层延迟、对流层延迟等误差信息，具体公式如下：
式中，表示站星几何距离，fi为卫星信号频率，i表示信号频率标志，r和s分别代表接收机和卫星系统，/>为测码伪距观测值，/>表示载波观测值，vr,i、kr,i和Kr,i分别表示接收机端的钟差、伪距延迟以及载波相位硬件延迟，属于接收机端误差；/>以及/>分别表示卫星端的钟差、伪距延迟以及载波相位硬件延迟，属于卫星端误差；/>和/>分别表示电离层延迟和对流层延迟，属于信号传输误差；λi为载波信号频率的波长，/>表示整周模糊度；/>和/>分别表示伪距的残差项载波相位的残差项；
公式(1)中的STD，可以表示为流体静力学延迟与非静力学延迟与其相应的映射函数的乘积，还需要顾及大气各项异性的情况，因此，STD的计算公式具体表示如下：
式中，ε为卫星高度角，为测站到卫星的方位角，Mh(ε)为天顶流体静力学延迟ZHD对应的映射函数，Mw(ε)为天顶非流体静力学延迟ZWD对应的映射函数，Mg(ε)为水平梯度映射函数，Gns和Gew分别表示水平梯度改正的南北方向分量与东西方向分量，γt表示对流层延迟残差；
在GNSS数据处理过程中，首先利用GNSS站点实测气压解算ZHD，其次将其作为先验值与其他先验参数一同代入线性化后的观测方程中，最后利用卡尔曼滤波方法对PPP函数模型进行解算，获取ZTD估值；
S2：GNSS ZTD的三维变分同化：
数值同化技术将GNSS ZTD及其误差信息作用于WRF模式，可获得改善后的大气状态估计值；基于WRF模式的三维变分同化，可理解为观测场与分析场、分析场与背景误差场之间的二次泛函极小化问题，其目标函数的表达式如下：
式中，x为待估计的大气状态向量，xb为大气状态的先验信息(背景场变量)，yo为实际观测值向量，H为观测算子；H(x)为关联背景场与实际观测值的投影函数，B为背景场xb的误差协方差矩阵，R为实际观测值的误差协方差矩阵，目标函数式亦可表示为：
J(x)＝Jb(x)+Jo(x) (4)
式中，Jb(x)表示分析量x偏离背景场xb的程度，Jo(x)表示由投影函数计算的分析量x偏离实际观测值的程度；目标函数J(x)的极小值是通过逐步迭代极小化方法得到的，最优分析场为代价函数极小化时的解，即待估计的真实大气状态x可表示为：
x＝xb+[B-1+HTR-1H]-1·HR-1[yo-H(xb)] (5)
基于WRF-DA三维变分同化系统对GNSS ZTD进行同化时，将GNSS ZTD观测资料加入式Jo(x)中，即：
式中，HZTD为GNSS ZTD的观测算子，R为GNSS ZTD的观测误差协方差矩阵；
S3：GNSS辅助WRF模式的气象参数预测：
WRF模式同化GNSS ZTD可提高预报气象参数的精度，其输出参数包含与计算对流层延迟密切相关的温度、气压、相对湿度、水汽混合比等；
大气压强可由WRF模式输出的空气基准态气压和扰动气压计算得到，其表达式如下：
P＝Pbase+Pper (7)
式中，Pbase、Pper分别表示模式输出的空气基准态气压和扰动气压(hPa)；
基于计算的大气压强P(hPa)，水蒸气分压e(hPa)可由模式输出的水蒸气混合比Q(kg·kg-1)转换得到：
大气温度可通过转换分层的空气总电位温度得到，其表达式为：
式中，Tθ为模式输出的潜在温度(K)，Rd为干空气气体常数，其值为287.053J·K-1·kg-1，cp为恒压下的比热容，Rd和cp的比值为泊松常数，Rd/cp＝2/7；
S4：GNSS卫星信号斜路径对流层延迟估算：
GNSS卫星信号的斜路径对流层延迟可基于卫星信号高度角、方位角及其在WRF三维网格中的截距对大气折射率进行积分获得，即获得GNSS卫星信号的SHD和SWD，进而求得STD，具体估算步骤如下：
S4.1：建立局部空间直角坐标系：
由于WRF模式所划分的网格在Z方向的分辨率一般为500m～1km，故可以忽略由垂线偏差所造成的影响，对椭球上的法线与垂线不加以严格区分；以研究区域左下角(西南角)的WRF格网点为坐标原点，北方向为X方向，东方向为Y方向，位势高方向为Z方向，建立局部空间直角坐标系；在此坐标系的基础上，计算出GNSS站点的空间直角坐标，并根据WRF网格空间分辨率、垂直分层计算每个WRF格网点的空间直角坐标；
在数值预报模式中，由于不同研究区域的下垫面情况不同，即使在一个很小的区域，地面的海拔也不一样，所以不适于按照位势高度来进行垂直分层；因此，WRF模式采用气压追随坐标eta在垂直方向上进行分层表示，取值范围从1到0，其表达式为：
eta＝(Plevel-Ptop)/(Pbot-Ptop) (10)
式中，Plevel为模式每一分层的气压(hPa)；Ptop表示模式顶层气压(hPa)；Pbot表示地表气压(hPa)；从eta坐标的表达式可以看出，地表的eta值为1，顶层的eta值为0；
而eta坐标对应的离地位势高度，可经过以下公式进行转换得到，即：
gmp＝(PH+PHB)/9.81-HGT(11)
式中，PH、PHB分别为模式输出的扰动位势高度和基准态位势高度；HGT为研究区域的地面海拔高度；位势高度和几何高度在数值上比较接近，则根据eta坐标对应的离地位势高度可以计算出WRF不同分层之间的高度；
S4.2：计算GNSS信号射线在WRF三维网格中的截距：
WRF模式研究区域内所有与X轴垂直的平面称为X型面，所有与Y轴垂直的平面称为Y型面，与Z轴垂直的平面称为Z型面；根据网格划分和垂直分层情况，WRF模式研究区域被离散化成若干X、Y和Z型面；GNSS信号射线穿过研究区域，与每一个型面有且只有一个交点；因此，计算GNSS信号射线与WRF网格交点即求解信号射线与所有X、Y、Z型面的交点坐标；若GNSS站点坐标为G，(x0,y0,z0)射线方向向量为S(m,n,p)，则信号射线与X型面的交点(xxi,yxi,zxi)可表示为：
GNSS信号射线与Y型面的交点(xyi,yyi,zyi)可表示为：
GNSS信号射线与Z型面的交点(xzi,yzi,zzi)可表示为：
则任意一条GNSS信号射线在网格(i,j,k)内的截距可表示为
式中，(xq,yq,zq)，(xq-1,yq-1,zq-1)表示GNSS信号射线与WRF网格(i,j,k)的两个交点坐标；
S4.3：STD估计：
(a)SHD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对网格内的大气干折射率进行积分，得到不同网格内的斜对流层静力学延迟，以WRF模式输出的气象参数直接积分计算大气层的静力学延迟，其具体公式如下：
SHD＝10-6·∫Nddl (16)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nd为WRF模式三维网格中的大气干折射率，其计算公式为：
式中，Pd为干大气压(hPa)，e为水蒸气分压(hPa)，T为大气温度(K)；
根据WRF模式垂直分层情况，将研究区域离散化为若干个独立的三维单元网格，并以每个网格中心点的大气干折射率代表三维网格的大气干折射率，则GNSS信号在斜路径上的对流层静力学延迟可表达为以下形式：
SHD＝∑ijk(dijk·Nd(ijk)) (18)
式中，不考虑地球曲率时，dijk表示GNSS信号射线在(i,j,k)网格的截距；考虑地球曲率时，dijk表示(i,j,k)网格的内插系数；Nd(ijk)表示(i,j,k)网格的大气干折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气干折射率插值得到；
(b)SWD估计：
基于GNSS信号射线穿过WRF模式不同网格的截距，对网格内的大气湿折射率进行积分，同理可得其计算公式如下：
SWD＝10-6·∫Nwdl (19)
式中，l表示GNSS信号射线由卫星传播到GNSS接收机的路径；Nw为WRF模式三维网格中的大气湿折射率，其计算公式为：
式中，Rw＝461.495J·K-1·kg-1为水蒸气的气体常数，T为大气温度(k)，k1＝77.6890K/hPa、k2＝71.2952K/hPa、k3＝375463K2/hPa为折射率常数；
同理可得，GNSS信号在斜路径上的非静力学延迟可表达为以下形式：
SWD＝∑ijk(dijk·Nw(ijk)) (21)
式中，Nw(ijk)表示(i,j,k)网格的大气湿折射率，其由GNSS信号射线穿过的WRF三维网格顶点的大气湿折射率插值得到；
最终，斜路径对流层延迟可由基于WRF模式输出气象参数计算的斜对流层静力学延迟和非静力学延迟求和得到，其具体表达式如下：
STD＝SHD+SWD(22)
式中，SHD、SWD分别为斜路径对流层静力学延迟和非静力学延迟；
S5：WRF辅助的实时PPP改善方法：
将基于WRF模式预报气象参数估算的STD加入到PPP定位观测方程中，对传统PPP定位方法进行改进，因此，顾及WRF预报气象参数改善STD的PPP载波相位观测方程表达如下：
式中，表示顾及STD项的载波相位观测值，/>为经过改正后的载波相位残差，其余各项不变；
通过对上式线性化，并联合多测站多卫星载波相位观测值构建GNSS PPP的函数模型和误差模型，利用卡尔曼滤波方法对PPP函数模型进行实时解算；
利用卡尔曼滤波方法对PPP载波相位观测方程进行解算时，通过将基于WRF模式预报气象参数估算的STD直接应用到载波相位观测方程中，减少了待估参数的个数；因此，滤波过程的各项权阵、滤波增益矩阵、协方差阵以及状态值矩阵都随之减少一阶，使得PPP解算过程更为简单高效，同时可以得到更高精度的PPP定位结果，并加快其收敛速度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310249580.7A CN116466376A (zh) | 2023-03-15 | 2023-03-15 | 数值预报模式辅助的实时ppp改善方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310249580.7A CN116466376A (zh) | 2023-03-15 | 2023-03-15 | 数值预报模式辅助的实时ppp改善方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116466376A true CN116466376A (zh) | 2023-07-21 |
Family
ID=87176177
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310249580.7A Pending CN116466376A (zh) | 2023-03-15 | 2023-03-15 | 数值预报模式辅助的实时ppp改善方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116466376A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116626730A (zh) * | 2023-07-24 | 2023-08-22 | 山东科技大学 | 一种顾及nwp的海上区域cors增强ppp方法 |
-
2023
- 2023-03-15 CN CN202310249580.7A patent/CN116466376A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116626730A (zh) * | 2023-07-24 | 2023-08-22 | 山东科技大学 | 一种顾及nwp的海上区域cors增强ppp方法 |
| CN116626730B (zh) * | 2023-07-24 | 2023-10-10 | 山东科技大学 | 一种顾及nwp的海上区域cors增强ppp方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114518586B (zh) | 一种基于球谐展开的gnss精密单点定位方法 | |
| CN110031877B (zh) | 一种基于grnn模型的区域nwp对流层延迟改正方法 | |
| CN109883313B (zh) | 一种基于单频gnss定位技术监测高铁大桥动态挠度的方法 | |
| CN109597040B (zh) | 一种星载sar影像无场几何定标方法 | |
| CN111596315B (zh) | 一种用于实时监测双频多星座星基增强系统性能的方法 | |
| CN111796309B (zh) | 导航卫星单频数据同步确定大气水汽和总电子含量的方法 | |
| CN116466376A (zh) | 数值预报模式辅助的实时ppp改善方法 | |
| CN115061167B (zh) | 一种适用于短距离大高差rtk的对流层延迟改正方法 | |
| CN113671505A (zh) | 一种基于系统几何误差补偿的合成孔径雷达立体定位方法 | |
| CN110146904B (zh) | 一种适用于区域电离层tec的精确建模方法 | |
| CN112711022B (zh) | 一种GNSS层析技术辅助的InSAR大气延迟改正方法 | |
| Abdelazeem et al. | MGR-DCB: a precise model for multi-constellation GNSS receiver differential code bias | |
| CN115980317B (zh) | 基于修正相位的地基gnss-r数据土壤水分估算方法 | |
| CN113009531A (zh) | 一种小尺度高精度的低空对流层大气折射率模型 | |
| CN115755103B (zh) | 一种抗差自适应的gnss水汽层析方法 | |
| CN115099159B (zh) | 基于神经网络且顾及地表差异的modis水汽反演方法 | |
| Liu et al. | The Impact of Different Mapping Function Models and Meteorological Parameter Calculation Methods on the Calculation Results of Single‐Frequency Precise Point Positioning with Increased Tropospheric Gradient | |
| CN114252875B (zh) | 一种成像高度计数据的高精度网格化方法 | |
| CN114910939B (zh) | 短距离大高差rtk中对流层延迟实测气象改正方法 | |
| CN115616615A (zh) | 一种PPP-B2b增强的低成本单频GNSS接收机精密定位方法 | |
| CN115755115A (zh) | 基于gnss对流层层析技术的ppp改善方法 | |
| CN115902968A (zh) | 基于北斗三号geo播发增强信息的ppp终端定位方法 | |
| CN112528213B (zh) | 基于低地球轨道卫星全球电离层总电子含量多层解析方法 | |
| CN115808686A (zh) | 基于波形重跟踪质量加权的垂线偏差格网化方法及其系统 | |
| CN116203610A (zh) | 两种gnss辅助wrf的实时ppp改善方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |