CN117194856B - 基于era5和cira的掩星弯曲角优化方法 - Google Patents

Abstract

发明提供了一种基于ERA5和CIRA的掩星弯曲角优化方法，包括以下步骤：S1、读取GNSS中性大气掩星的双频附加相位文件，获取特征信息；S2、根据获取的特征信息，分别利用几何光学方法和全谱法计算弯曲角；S3、对利用几何光学方法和全谱法计算分别得到的弯曲角进行拼接处理，得到双频弯曲角组合；S4、对双频弯曲角组合消电离层影响得到观测弯曲角；S5、对得到的观测弯曲角按权重函数进行计算得到优化弯曲角；S6、进行ABEL反演得到大气折射率廓线。本发明有益效果：基于ERA5和CIRA的掩星弯曲角优化方法，提高大气掩星弯曲角廓线和折射率廓线的精度，提升其在数值预报和再分析系统中的影响，以及数值预报精度。

Description

基于ERA5和CIRA的掩星弯曲角优化方法

技术领域

本发明属于掩星观测领域，尤其是涉及一种基于ERA5和CIRA的掩星弯曲角优化方法。

背景技术

GNSS掩星廓线是数值天气预报和再分析资料重要数据源之一，大气弯曲角廓线和大气折射率廓线可直接应用于数值预报同化系统中，大气折射率由弯曲角经过反演得到。大气弯曲角的精度直接关系到大气折射率产品精度及数值预报同化后的应用效果，故获取高精度的掩星大气弯曲角廓线十分重要。GNSS大气掩星弯曲角在上部（20km以上）误差较大，通常使用CIRA气候态对弯曲角进行优化，但该气候态精度不是很高，导致计算的大气弯曲角和折射率上部（20km以上）误差较大，影响其在数值预报和再分析系统中的应用。

发明内容

有鉴于此，本发明旨在提出一种基于ERA5和CIRA的掩星弯曲角优化方法，高精度的欧洲中期天气预报中心（ECMWF）ERA5月平均模型为37层，高度范围0~45km，CIRA气候态模型高度范围0~120km，联合30年ERA5月平均模型和CIRA，对上部（20km以上）大气弯曲角廓线精度进行优化，进而反演得到高精度的折射率廓线，提高大气弯曲角廓线和折射率廓线的精度。

为达到上述目的，本发明的技术方案是这样实现的：

基于ERA5和CIRA的掩星弯曲角优化方法，包括以下步骤：

S1、读取GNSS中性大气掩星的双频附加相位文件，获取特征信息；

S2、根据获取的特征信息，分别利用几何光学方法和全谱法计算弯曲角廓线；

S3、对利用几何光学方法和全谱法计算分别得到的弯曲角廓线进行拼接处理，得到双频弯曲角廓线组合；

S4、对双频弯曲角廓线组合消电离层影响得到观测弯曲角廓线；

S5、对得到的观测弯曲角廓线按权重函数进行计算得到优化弯曲角廓线；

S6、进行ABEL反演得到大气折射率廓线。

进一步的，在步骤S2中，利用几何光学方法计算弯曲角廓线方法如下：

将地心、LEO、掩星之间相互连线，利用步骤S1中获取的特征信息及预设的弯曲角廓线模型进行逼近计算，通过CIRA模型中温湿压通过气体方程计算得到大气折射率，再积分得到预设的弯曲角廓线模型，由地心、LEO、掩星之间几何关系得到以地心为中心的LEO与掩星之间夹角θ、LEO天顶角Φ₁、掩星天顶角Φ₂及掩星射线方向天顶角增量dΦ；

由几何关系及附加相位延迟和运动引起的多普勒计算得出模型多普勒dfrm信息；

计算双频附加相位对应的附加多普勒dfrm1,dfrm2；

通过模型多普勒dfrm与附加多普勒dfrm1,dfrm2之间差值对无效信号进行截断；

掩星反演过程中采用信号沿直线传播和局部球对称假设，计算曲率中心和曲率半径并将LEO和掩星坐标改正到曲率中心上；

将LEO和掩星的惯性通过角度旋转转换到地固坐标系下，根据几何关系可知，GNSS信号由掩星发射经过大气发生弯曲的角度e计算如公式（1）所示：

e=Φ₁+Φ₂+θ-π（1）。

进一步的，在步骤S2中，利用全谱法计算弯曲角廓线方法如下：

利用附加相位文件输入的信噪比和附加相位进行滤波处理，得到弯曲角廓线是随碰撞参数变化的单调函数。

进一步的，在步骤S3中，获取双频弯曲角廓线组合的方法如下：

以设定高度为界对几何光学和全谱法计算的弯曲角廓线进行拼接，设定高度以下采用全谱法计算的大气弯曲角廓线，设定高度以上采用几何光学计算的大气弯曲角廓线，拼接完成的双频弯曲角廓线分别表示为e1,e2。

进一步的，在步骤S4中，消电离层影响的方法如下：

根据GNSS第一频点和第二频点观测方程中电离层延迟特性进行消电离层组合，组合后的大气观测弯曲角廓线B_o计算公式如公式（2）所示，通过掩星信噪比计算观测弯曲角廓线B_o对应的误差协方差B_{o_var}：

（2）；

其中，,/>分别为双频频率。

进一步的，在步骤S5中，对弯曲角廓线优化方法如下：

根据掩星事件所处月份及经纬度信息读取CIRA气候态的温度、湿度、压强数据，并插值计算其对应的大气折射率N_CIRA，通过正演方法计算出CIRA模型弯曲角廓线B_CIRA；

根据掩星事件所处月份及经纬度信息读取平均化处理后的格网数据，插值并依据公式（3）计算出其对应的大气折射率N_era5，通过正演方法计算出ERA5模型弯曲角廓线B_ERA5，其中常数k₁=77.6K/mb，k₂=3.73×10⁵K²/mb，P和Pw分别为大气总气压和水汽分压，T为大气绝对温度：

N_era5 （3）

通过弯曲角廓线误差协方差分段组建权重函数w如公式（4）所示，其中h为海拔高度，单位km；

使用ERA5和CIRA模型优化弯曲角廓线如公式（5）所示，计算得到优化后的弯曲角廓线B_opt(a)，弯曲角廓线为碰撞参数的函数：

（4）

（5）；

其中，B_o为双频组合后的弯曲角廓线，B_{o_var}为弯曲角廓线B_o对应的误差协方差，B_CIRA和B_ERA5分别为利用CIRA模型和ERA5模型计算得到的弯曲角廓线。

进一步的，在步骤S6中，方法如下：

通过ABEL逆变换将弯曲角廓线反演得到大气折射率N（x），计算如公式（6）所示：

（6）；

其中a为碰撞参数，x为信号弯曲圆弧半径与折射率乘积。

进一步的，本方案公开了一种电子设备，包括处理器以及与处理器通信连接，且用于存储所述处理器可执行指令的存储器，所述处理器用于执行基于ERA5和CIRA的掩星弯曲角优化方法。

进一步的，本方案公开了一种服务器，包括至少一个处理器，以及与所述处理器通信连接的存储器，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述处理器执行，以使所述至少一个处理器执行基于ERA5和CIRA的掩星弯曲角优化方法。

进一步的，本方案公开了一种计算机可读取存储介质，存储有计算机程序，所述计算机程序被处理器执行时实现基于ERA5和CIRA的掩星弯曲角优化方法。

相对于现有技术，本发明所述的基于ERA5和CIRA的掩星弯曲角优化方法具有以下有益效果：

本发明所述的基于ERA5和CIRA的掩星弯曲角优化方法，提高大气掩星弯曲角廓线和折射率廓线的精度，提升其在数值预报和再分析系统中的影响，以及数值预报精度。

附图说明

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为本发明实施例所述的技术路线图示意图；

图2为本发明实施例所述的气候态模型数据格式示意图；

图3为本发明实施例所述的气候态模型数据时间跨度示意图；

图4为本发明实施例所述的气候参数种类示意图；

图5为本发明实施例所述的CIRA气象参数纬度和高度序列示意图。

具体实施方式

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面将参考附图并结合实施例来详细说明本发明。

如图1所示为方法技术路线图：

读取GNSS（中国北斗、美国GPS、俄罗斯GLONASS和欧洲GALILEO）中性大气掩星的双频附加相位文件，利用几何光学方法和全谱法分别计算出大气掩星弯曲角廓线，对弯曲角进行拼接处理，采用双频弯曲角组合消除电离层对中性大气掩星弯曲角的影响得到观测弯曲角。读取CIRA气候态模型数据和30年ERA5月平均模型数据，分别计算CIRA模型弯曲角和ERA5模型弯曲角，采用分段权重函数对观测弯曲角、CIRA模型弯曲角和ERA5模型弯曲角分高度按权重函数进行计算得到优化弯曲角，进而进行ABEL反演得到大气折射率廓线。首先需要说明的是，由于弯曲角是随高度变化的廓线，所以也叫作弯曲角廓线，具体过程如下：

a)读取大气附加相位文件

读取掩星大气附加相位文件，获取掩星发生时间、双频信噪比、惯性坐标下LEO位置、惯性坐标下LEO速度、惯性坐标下掩星位置、惯性坐标下掩星速度、双频附加相位延迟和开环相位模型等信息。

b)几何光学计算弯曲角

地心、LEO、掩星三者之间构成掩平面的几何关系，将地心、LEO、掩星之间相互连线，利用LEO和掩星位置信息及弯曲角模型进行逼近计算，得到以地心为中心的LEO与掩星之间夹角θ、LEO天顶角Φ₁、掩星天顶角Φ₂及掩星射线方向天顶角增量dΦ，由几何关系及附加相位延迟和运动引起的多普勒可计算出模型多普勒dfrm信息。计算双频附加相位对应的附加多普勒dfrm1,dfrm2。通过模型多普勒dfrm与附加多普勒dfrm1,dfrm2之间差值对无效信号进行截断。掩星反演过程中采用信号沿直线传播和局部球对称假设，计算曲率中心和曲率半径并将LEO和掩星坐标改正到曲率中心上。将LEO和掩星的惯性通过角度旋转转换到地固坐标系下。根据几何关系可知，GNSS信号由掩星发射经过大气发生弯曲的角度e计算如公式（1）所示。

e=Φ₁+Φ₂+θ-π（1）

c)全谱法计算弯曲角

10km以下部分水汽较严重，容易造成多路径现象，该部分弯曲角采用全谱法进行计算。利用附加相位文件输入的信噪比和附加相位进行滤波等处理，消除GNSS掩星信号中的多路径问题，避免ABEL反演过程中的多值不单调问题，得到的弯曲角是随碰撞参数变化的单调函数。

d)弯曲角拼接

以10km为界对几何光学和全谱法计算的弯曲角进行拼接，10km以下采用全谱法计算的大气弯曲角，10km以上采用几何光学计算的大气弯曲角，拼接完成的双频弯曲角分别表示为e1,e2。

e)消电离层影响

根据GNSS第一频点和第二频点观测方程中电离层延迟特性进行消电离层组合，组合后的大气观测弯曲角B_o计算公式如公式（2）所示。通过掩星信噪比计算观测弯曲角B_o对应的误差协方差B_{o_var}。

（2）

f)弯曲角优化

CIRA气候态模型格式如图2所示，按照月份（图3所示）、气象参数（图4所示）、纬度序列和高度序列（图5所示）四个维度进行分类存储，月份为1月-12月，气象参数包括纬向平均温度、纬向平均湿度和纬向平均压强，纬度序列为南纬90°~北纬90°且间隔为10°，高度序列为0~120km的37层序列。根据掩星事件所处月份及经纬度信息读取CIRA气候态的温度、湿度、压强数据，并插值计算其对应的大气折射率N_CIRA，通过正演方法计算出CIRA模型弯曲角B_CIRA。

从ECMWF官网下载近30年的再分析资料ERA5月平均模式数据，该数据中包含温度、湿度、压强的四维格网数据，计算出12个月近30年的平均温度、平均湿度、平均压强格网数据，格网分辨率0.25°×0.25°，37层（0~45km）。根据掩星事件所处月份及经纬度信息读取平均化处理后的格网数据，插值并依据公式（3）计算出其对应的大气折射率N_era5，通过正演方法计算出ERA5模型弯曲角B_ERA5。其中常数k₁=77.6K/mb，k₂=3.73×10⁵K²/mb。

（3）

通过弯曲角误差协方差分段组建权重函数w如公式（4）所示，其中h为海拔高度，单位km，可由碰撞参数a与曲率半径等计算得到。使用ERA5和CIRA模型优化弯曲角如公式（5）所示，0~20km弯曲角完全使用观测弯曲角B_o，20km~45km使用ERA530年月平均模型进行优化，45km~120km使用CIRA模型进行优化，计算得到优化后的弯曲角B_opt(a)，弯曲角为碰撞参数的函数。

（4）

（5）

g)ABEL反演折射率

大气弯曲角与大气折射率之间存在关系，通过ABEL逆变换可将弯曲角反演得到大气折射率N（x），计算如公式（6）所示。

（6）

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及方法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

在本申请所提供的几个实施例中，应该理解到，所揭露的方法和系统，可以通过其它的方式实现。例如，以上所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。上述单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本发明实施例方案的目的。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，包括以下步骤：

S1、读取GNSS中性大气掩星的双频附加相位文件，获取特征信息；

S2、根据获取的特征信息，分别利用几何光学方法和全谱法计算弯曲角廓线；

S3、对利用几何光学方法和全谱法计算分别得到的弯曲角廓线进行拼接处理，得到双频弯曲角廓线组合；

S4、对双频弯曲角廓线组合消电离层影响得到观测弯曲角廓线；

S5、对得到的观测弯曲角廓线按权重函数进行计算得到优化弯曲角廓线；

S6、进行ABEL反演得到大气折射率廓线；

在步骤S5中，对弯曲角廓线优化方法如下：

根据掩星事件所处月份及经纬度信息读取CIRA气候态的温度、湿度、压强数据，并插值计算其对应的大气折射率N_CIRA，通过正演方法计算出CIRA模型弯曲角廓线B_CIRA；

根据掩星事件所处月份及经纬度信息读取平均化处理后的格网数据，插值并依据公式(3)计算出其对应的大气折射率N_era5，通过正演方法计算出ERA5模型弯曲角廓线B_ERA5，其中常数k₁＝77.6K/mb，k₂＝3.73×10⁵K²/mb，P和Pw分别为大气总气压和水汽分压，T为大气绝对温度：

通过弯曲角廓线误差协方差分段组建权重函数w如公式(4)所示，其中h为海拔高度，单位km；

使用ERA5和CIRA模型优化弯曲角廓线如公式(5)所示，计算得到优化后的弯曲角廓线B_opt(a)，弯曲角廓线为碰撞参数的函数：

其中，B_o为双频组合后的弯曲角廓线，B_{o_var}为弯曲角廓线B_o对应的误差协方差，B_CIRA和B_ERA5分别为利用CIRA模型和ERA5模型计算得到的弯曲角廓线。

2.根据权利要求1所述的基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，在步骤S2中，利用几何光学方法计算弯曲角廓线方法如下：

将地心、LEO、掩星之间相互连线，利用步骤S1中获取的特征信息及预设的弯曲角廓线模型进行逼近计算，通过气体方程对CIRA中温湿压进行计算得到大气折射率，再积分得到预设的弯曲角廓线模型，由地心、LEO、掩星之间几何关系得到以地心为中心的LEO与掩星之间夹角θ、LEO天顶角Φ₁、掩星天顶角Φ₂及掩星射线方向天顶角增量dΦ；

由几何关系及附加相位延迟和运动引起的多普勒计算得出模型多普勒dfrm信息；

计算双频附加相位对应的附加多普勒dfrm1,dfrm2；

通过模型多普勒dfrm与附加多普勒dfrm1,dfrm2之间差值对无效信号进行截断；

掩星反演过程中采用信号沿直线传播和局部球对称假设，计算曲率中心和曲率半径并将LEO和掩星坐标改正到曲率中心上；

将LEO和掩星的惯性通过角度旋转转换到地固坐标系下，根据几何关系可知，GNSS信号由掩星发射经过大气发生弯曲的角度e计算如公式(1)所示：

e＝Φ₁+Φ₂+θ-π (1)。

3.根据权利要求1所述的基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，在步骤S2中，利用全谱法计算弯曲角廓线方法如下：

利用附加相位文件输入的信噪比和附加相位进行滤波处理，得到弯曲角廓线是随碰撞参数变化的单调函数。

4.根据权利要求1所述的基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，在步骤S3中，获取双频弯曲角廓线组合的方法如下：

以设定高度为界对几何光学和全谱法计算的弯曲角廓线进行拼接，设定高度以下采用全谱法计算的大气弯曲角廓线，设定高度以上采用几何光学计算的大气弯曲角廓线，拼接完成的双频弯曲角廓线分别表示为e1,e2。

5.根据权利要求4所述的基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，在步骤S4中，消电离层影响的方法如下：

根据GNSS第一频点和第二频点观测方程中电离层延迟特性进行消电离层组合，组合后的大气观测弯曲角廓线B_o计算公式如公式(2)所示，通过掩星信噪比计算观测弯曲角廓线B_o对应的误差协方差B_{o_var}：

其中，f₁,f₂分别为双频频率。

6.根据权利要求1所述的基于ERA5和CIRA的掩星弯曲角优化方法，其特征在于，在步骤S6中，方法如下：

通过ABEL逆变换将弯曲角廓线反演得到大气折射率N(x)，计算如公式(6)所示：

其中a为碰撞参数，x为信号弯曲圆弧半径与折射率乘积。

7.一种电子设备，包括处理器以及与处理器通信连接，且用于存储所述处理器可执行指令的存储器，其特征在于：所述处理器用于执行上述权利要求1-6任一所述的基于ERA5和CIRA的掩星弯曲角优化方法。

8.一种服务器，其特征在于：包括至少一个处理器，以及与所述处理器通信连接的存储器，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述处理器执行，以使所述至少一个处理器执行如权利要求1-6任一所述的基于ERA5和CIRA的掩星弯曲角优化方法。

9.一种计算机可读取存储介质，存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现权利要求1-6任一项所述的基于ERA5和CIRA的掩星弯曲角优化方法。