CN105372720A - 一种基于导航卫星掩星资料同化的spsa气象参数解算方法 - Google Patents

Abstract

本发明公开了一种基于导航卫星掩星资料同化的SPSA气象参数解算方法，其特征是按如下步骤进行：首先建立掩星资料同化的代价函数，然后计算代价函数的一阶近似梯度，接着利用近似梯度更新大气状态向量，最后求解当前大气状态向量下的代价函数值，并判断是否满足停止迭代条件，若满足条件，即可获得大气状态向量的解。本发明避免了设计观测算子的切线性算子和伴随算子，极大地简化了代价函数的求解算法，提高了求解大气状态向量的效率和可行性，有广阔的应用前景。

Description

一种基于导航卫星掩星资料同化的SPSA气象参数解算方法

技术领域

本发明属于GPS气象学领域，具体地说是一种基于导航卫星掩星资料同化的SPSA气象参数解算方法。

背景技术

无线电掩星观测技术是GPS气象学的一个重要学科，是利用信号的延迟和弯曲来分析大气参数。GPS导航卫星连续发射L1和L2两个波段的载波信号。载波信号穿过大气层时会发生延迟和弯曲，被低轨道卫星上的卫星接收机接。在信号接收端对两个波段的载波信号进行线性组合，可以消除电离层的影响，并计算出中性大气层引起的信号传输延迟量和信号弯曲角廓线；进一步通过Abel积分变换可以得到大气折射率廓线。弯曲角和大气折射率就是最常用的GPS掩星资料，它们具有精度高、时间分辨率高、垂直分辨率高、全球覆盖、全天候观测、易于在业务应用中普及等优点。这些掩星资料同化进数值天气预报模式中，可以明显改善数值天气预报的初始场，提高数值天气预报的预报能力。掩星资料同化技术中代价函数求解是关键问题。

现有技术中，一般采用LM迭代法或准牛顿迭代法求解代价函数最优解。这两种方法都必须先求解代价函数的一阶导数和二阶导数，然后对矩阵方程进行迭代计算，直至达到停止条件。当观测算子比较复杂时，无法直接计算出一阶导数和二阶导数。需要先设计出观测算子的切线性算子及其伴随算子，再利用LM迭代法或拟牛顿方法求解。该过程大大增加了算法的复杂度，而且切线性算子和伴随算子的设计是个难点。

发明内容

本发明是为克服上述现有技术的不足之处，提出了一种基于导航卫星掩星资料同化的SPSA气象参数求解算法，以期避免设计观测算子的切线性算子和伴随算子，极大地简化代价函数的求解算法，从而提高求解大气状态向量的效率和可行性。

本发明为达到上述发明目的，采用如下技术方案：

本发明一种基于导航卫星掩星资料同化的SPSA气象参数解算方法的特点是按如下步骤进行：

步骤1、建立如式(1)所示的掩星资料同化的代价函数J(x)：

$J\left(x\right)=\frac{1}{2}{(x-x^b)}^T{B}^{-1}(x-x^b)+\frac{1}{2}{(y^o-H(x\left)\right)}^T{R}^{-1}(y^o-H(x\left)\right)---\left(1\right)$

式(1)中，x为大气状态向量，即为待求解的参数；并有：x＝(T,p,q)^T，T表示大气温度，p表示大气压强，q表示大气比湿；x^b表示背景场的状态向量，并有x^b＝(T^b,p^b,q^b)^T；T^b表示背景场的大气温度，p^b表示背景场的大气压强，q^b表示背景场的大气比湿；B为背景场的误差协方差矩阵；y^o是观测场掩星资料；H(x)是观测算子；R是观测场的误差协方差矩阵；

步骤2、定义迭代次数为k，并初始化k＝1；定义阈值为J_min；令x_k＝x^b；

步骤3、利用式(1)获得在x＝x_k处的第k次迭代的代价函数J_k(x)，并判断k≥2且|J_k(x)-J_k-1(x)|≤J_min是否成立，若成立，则停止迭代，获得第k次迭代大气状态向量x_k；否则，执行步骤4；

步骤4、利用式(2)求解式(1)在x＝x_k处的近似导数J_k'：

${J_k}^{\′}=\frac{{J_k}^{+}-{J_k}^{-}}{2\·c_k}\×{\mathrm{\Δ}}_k---\left(2\right)$

并有：

J_k ⁺＝J_k(x+c_k.·Δ_k)(3)

J_k ^-＝J_k(x-c_k.·Δ_k)(4)

式(2)-式(4)中，c_k表示第k次迭代的扰动参数向量，并有：

$c_k=\frac{C}{{(k+1)}^{\mathrm{\γ}}}---\left(5\right)$

式(5)中，c表示初始扰动参数向量；γ表示非负系数；

式(2)-式(4)中，Δ_k表示第k次迭代的m维随机扰动向量，并由MonteCarlo方法产生，并有：

Δ_k＝|Δ_k1,…,Δ_km|^T(6)

步骤5、利用式(7)获得第k+1次迭代的大气状态向量x_k+1：

x_k+1＝x_k+a_k.·J'(7)

式(7)中，a_k表示第k次的迭代步长向量；并有：

$a_k=\frac{a}{{(k+1+A)}^{\mathrm{\α}}}---\left(8\right)$

式(8)中，a表示初始步长向量；A和α表示非负系数；

步骤6、将k+1赋值给k；并返回步骤3顺序执行。

与现有技术相比，本发明的有益效果体现在：

1、本发明通过建立掩星资料同化的代价函数，然后计算代价函数的一阶近似梯度，接着利用近似梯度更新大气状态向量，最后求解当前大气状态向量下的代价函数值，并判断是否满足停止迭代条件，若满足条件，即可获得大气状态向量的解，从而提高了求解大气状态参数的效率和可行性，具有广阔的应用前景。

2、本发明通过SPSA算法求解代价函数，简化了代价函数的求解算法，并通过代价函数的一阶近似梯度更新大气状态向量，不直接求解一阶和二阶导数，避免了设计观测算子的切线性算子和伴随算子。

3、本发明中SPSA算法的扰动参数和迭代步长均为与大气状态向量同维度的向量，由根据大气参数向量各个值对代价函数的影响度和方差确定，能更快的收敛到大气状态参数向量的最优估计。

附图说明

图1为现有技术无线电掩星观测原理示意图；

图2为本发明SPSA求解代价函数的流程图。

具体实施方式

本实施例中，掩星资料同化的原理如图1所示，GPS导航卫星连续发射L1和L2两个波段的载波信号穿过大气层发生延迟和弯曲后被低轨道卫星上的卫星接收机接。信号接收端可以计算出中性大气层引起的载波信号的传输延迟量和弯曲角廓线，并进一步得到大气折射率廓线。将弯曲角和大气折射率等掩星资料同化进数值天气预报模式中，可以解算大气状态参数向量。如图2所示，一种基于导航卫星掩星资料同化的SPSA气象参数求解算法是按如下步骤进行：

步骤1、当前数值天气预报采用的业务系统为同化系统，同化系统几乎可以同化所有类型的观测资料，使得大量的非常规的观测资料能够得以更好的应用，弥补常规观测资料的不足。掩星资料同化就是将大气折射率和弯曲角等掩星资料同化进数值天气模式中。同化方法的本质是构建模式的代价函数，利用代价函数的极小化来求解大气状态的最优解，属于复杂的多维函数寻优问题。假设观测和背景场具有无偏的高斯误差条件下，可以利用贝叶斯定理推导出代价函数；则建立如式(1)所示的掩星资料同化的代价函数J(x)：

$J\left(x\right)=\frac{1}{2}{(x-x^b)}^T{B}^{-1}(x-x^b)+\frac{1}{2}{(y^o-H(x\left)\right)}^T{R}^{-1}(y^o-H(x\left)\right)---\left(1\right)$

式(1)中，x为大气状态向量，即为待求解的参数；并有：x＝(T,p,q)^T，T表示大气温度，p表示大气压强，q表示大气比湿；x^b表示背景场的状态向量，并有x^b＝(T^b,p^b,q^b)^T；T^b表示背景场的大气温度，p^b表示背景场的大气压强，q^b表示背景场的大气比湿；B为背景场的误差协方差矩阵；y^o是观测场掩星资料；R是观测场的误差协方差矩阵；H(x)是观测算子；

观测算子是将背景场的大气状态向量映射到观测场，同观测值进行比较。根据同化数据的不同，主要分为折射率算子和弯曲角算子，每一种观测算子又可分为局地算子和非局地算子。例如，大气折射率正掩算子的形式如下，H(x)＝IH_N(x)，H_N是利用背景场上的温度，大气压和比湿计算大气折射率。假设忽略非理想气体的影响，背景场的大气折射率N的计算公式H_N可表达为：e为水汽分压，k₁，k'₂，k₃是经验常数，由实验测定。I是将背景场上的大气折射率插值到观测场，使背景场和观测场数据在同一坐标系和同一势高下进行比较。

步骤2、定义迭代次数为k，并初始化k＝1；定义阈值为J_min；令x_k＝x^b；

步骤3、利用式(1)获得在x＝x_k处的第k次迭代的代价函数J_k(x)，并判断k≥2且|J_k(x)-J_k-1(x)|≤J_min是否成立，若成立，则停止迭代，获得第k次迭代大气状态向量x_k；否则，执行步骤4；

步骤4、SPSA算法利用代价函数的一阶近似导数更新未知大气向量的值，所以首先要利用式(2)求解式(1)在x＝x_k处的近似导数J_k'：

${J_k}^{\′}=\frac{{J_k}^{+}-{J_k}^{-}}{2\·c_k}\×{\mathrm{\Δ}}_k---\left(2\right)$

若只求一次一阶近似导数，在某些情况下，一阶近似导数与真实导数相比，会存在较大的误差。为提高算法在噪声条件下的性能，可采用梯度平均法进行导数估计，，假设求n次导数，则有J'＝n^-1∑J'k，并有：

J_k ⁺＝J_k(x+c_k.·Δ_k)(3)

J_k ^-＝J_k(x-c_k.·Δ_k)(4)

式(2)-式(4)中，c_k表示第k次迭代的扰动参数向量，并有：

$c_k=\frac{c}{{(k+1)}^{\mathrm{\γ}}}---\left(5\right)$

式(5)中，c表示初始扰动参数向量，在常规的SPSA算法中，初始扰动参数一般为一个常数值，但是由于大气状态参数向量中各个参数相同变化值对代价函数的影响明显不同，所以，根据不同参数对代价函数的影响，确定初始扰动参数为一个和大气状态参数向量同维度的向量；γ表示非负系数；

式(2)-式(4)中，Δ_k表示第k次迭代的m维随机扰动向量，并由MonteCarlo方法产生，并有：

Δ_k＝|Δ_k1,…,Δ_km|^T(6)

步骤5、利用式(7)获得第k+1次迭代的大气状态向量x_k+1：

x_k+1＝x_k+a_k.·J'(7)

式(7)中，a_k表示第k次的迭代步长向量；并有：

$a_k=\frac{a}{{(k+1+A)}^{\mathrm{\α}}}---\left(8\right)$

式(8)中，a表示初始步长向量；变分同化代价函数中大气状态参数向量的解与背景场和观测场的误差协方差矩阵相关。若在每次更新估计值时，每一维参数的迭代步长值大小一样，则大气状态参数向量逼近最优解比较慢。对步长进行改进，改进策略为将初始步长参数变为一个向量，每一维的值根据背景场和观测场的误差协方差矩阵确定，已达到算法能准确迅速地收敛到最优解的目的。A和α表示非负系数；

步骤6、将k+1赋值给k；并返回步骤3顺序执行。

Claims (1)

1.一种基于导航卫星掩星资料同化的SPSA气象参数解算方法，其特征是按如下步骤进行：

步骤1、建立如式(1)所示的掩星资料同化的代价函数J(x)：

$J\left(x\right)=\frac{1}{2}{(x-x^b)}^T{B}^{-1}(x-x^b)+\frac{1}{2}{(y^o-H(x\left)\right)}^T{R}^{-1}(y^o-H(x\left)\right)---\left(1\right)$

式(1)中，x为大气状态向量，即为待求解的参数；并有：x＝(T,p,q)^T，T表示大气温度，p表示大气压强，q表示大气比湿；x^b表示背景场的状态向量，并有x^b＝(T^b,p^b,q^b)^T；T^b表示背景场的大气温度，p^b表示背景场的大气压强，q^b表示背景场的大气比湿；B为背景场的误差协方差矩阵；y^o是观测场掩星资料；H(x)是观测算子；R是观测场的误差协方差矩阵；

步骤2、定义迭代次数为k，并初始化k＝1；定义阈值为J_min；令x_k＝x^b；

步骤3、利用式(1)获得在x＝x_k处的第k次迭代的代价函数J_k(x)，并判断k≥2且|J_k(x)-J_k-1(x)|≤J_min是否成立，若成立，则停止迭代，获得第k次迭代大气状态向量x_k；否则，执行步骤4；

步骤4、利用式(2)求解式(1)在x＝x_k处的近似导数J_k'：

${J_k}^{\′}=\frac{{J_k}^{+}-{J_k}^{-}}{2\·c_k}\×{\mathrm{\Δ}}_k---\left(2\right)$

并有：

J_k ⁺＝J_k(x+c_k.·Δ_k)(3)

J_k ^-＝J_k(x-c_k.·Δ_k)(4)

式(2)-式(4)中，c_k表示第k次迭代的扰动参数向量，并有：

$c_k=\frac{c}{{(k+1)}^{\mathrm{\γ}}}---\left(5\right)$

式(5)中，c表示初始扰动参数向量；γ表示非负系数；

式(2)-式(4)中，Δ_k表示第k次迭代的m维随机扰动向量，并由MonteCarlo方法产生，并有：

Δ_k＝|Δ_k1,…,Δ_km|^T(6)

步骤5、利用式(7)获得第k+1次迭代的大气状态向量x_k+1：

x_k+1＝x_k+a_k.·J'(7)

式(7)中，a_k表示第k次的迭代步长向量；并有：

$a_k=\frac{a}{{(k+1+A)}^{\mathrm{\α}}}---\left(8\right)$

式(8)中，a表示初始步长向量；A和α表示非负系数；

步骤6、将k+1赋值给k；并返回步骤3顺序执行。