CN109145251A - 一种改进型同步扰动随机逼近算法的大气参数求解方法 - Google Patents

Abstract

本发明公开了一种改进型同步扰动随机逼近算法的大气参数求解方法，其特征是按如下步骤进行：首先建立掩星资料同化的大气参数求解代价函数；然后采用同步扰动随机逼近算法求解代价函数的最优解，以代价函数的近似梯度引导算法的搜索过程，避免了复杂的代价函数梯度求解。在求解过程中，对搜索空间网格化划分可改进算法的全局搜索能力；同时对初始步长的自适应改进可以提高算法的收敛速度。本发明简化了掩星资料同化技术中代价函数的求解，提高了求解大气参数向量的效率和可行性，有广阔的应用前景。

Description

一种改进型同步扰动随机逼近算法的大气参数求解方法

技术领域

本发明属于GPS气象学领域，具体地说是一种改进型同步扰动随机逼近算法的大气参数求解方法。

背景技术

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

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

发明内容

本发明是为克服上述现有技术的不足之处，提出了一种改进型同步扰动随机逼近算法的大气参数求解方法，以期简化掩星资料同化技术中代价函数的求解，提高求解大气参数向量的效率和可行性。

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

本发明一种改进型同步扰动随机逼近算法的大气参数求解方法的特点是按如下步骤进行：

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

式(1)中，x为大气参数向量，即为待求解的参数，并有：x＝(t,p,q)^T，其中t表示大气温度，且t∈[t_i,t_f]，t_i表示大气温度t的取值范围的下限，t_f表示大气温度t的取值范围的上限；p表示大气压强，且p∈[p_i,p_f]，p_i表示大气压强p的取值范围的下限，p_f表示大气压强p的取值范围的上限；q表示大气比湿，且q∈[q_i,q_f]，q_i表示大气比湿q的取值范围的上限，q_f表示大气比湿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、假设大气参数向量x的三个元素t,p,q在空间中相互垂直，以点O(0,0,0)为原点、以大气温度t为x轴、大气压强p为y轴、大气比湿q为z轴建立关于大气参数向量x的空间直角坐标系O-xyz，并根据大气参数向量x的三个元素t,p,q的取值范围，构建改进型同步扰动随机逼近算法的搜索空间，所述搜索空间即为空间直角坐标系O-xyz上的一个立方体；

步骤3、将所述算法搜索空间的每一维均分为N份，从而得到N³个网格；定义m,n,s分别为任意一个网格在大气温度t、大气压强p、大气比湿q三个维度上的序号，m＝1,···,N；n＝1,···,N；s＝1,···,N；并利用式(2)得到在所述三维空间坐标系O-xyz中空间网格的中心坐标x_m,n,s：

步骤4、利用改进型同步扰动随机逼近算法在搜索空间的每个网格中求解最优解：

步骤4.1、定义网格序号为g，g＝1,2,…,N³；并初始化g＝1；定义代价函数阈值为J_min；

步骤4.2、判断s是否等于N，否，s＝s+1,m＝1,n＝1，执行步骤4.3；是，停止迭代，得到代价函数最优解为：

x＝argmin(J(x_g)),g＝1,2,...,N³ (3)

步骤4.3、判断n是否等于N，否，n＝n+1,m＝1，执行步骤4.4；是，执行步骤4.2；

步骤4.4、判断m是否等于N，否，m＝m+1，执行步骤4.5；是，执行步骤4.3；

步骤4.5、定义迭代次数为k，并初始化k＝1；令当前网格的初始迭代点坐标x_k＝x_m,n,s，即网格中心点；令当前第g个网格的参数向量为x_g；

步骤4.6、利用式(1)获得在第k个迭代点处的大气参数向量x_k的代价函数J(x_k)，并判断k≥2且|J(x_k)-J(x_k-1)|≤J_min是否成立，若成立，则获得第k个迭代点处的大气参数向量x_k，即为当前第g个网格的最优解，并赋值给第g个网格的参数向量x_g，将g+1赋值给g，执行步骤4.4；否则，执行步骤4.7；

步骤4.7、利用式(4)求解式(1)在第k个迭代点处的大气参数向量x_k的近似梯度J′_k：

并有：

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

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

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

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

式(4)-式(6)中，Δ_k表示第k次迭代的d维随机扰动向量，并有：

步骤4.8、利用式(9)获得第k+1个迭代点处的大气参数向量x_k+1：

x_k+1＝x_k-a_k.·J′ (9)

式(9)中，a_k表示第k个迭代点处的步长向量，并有：

式(10)中，A和α表示非负系数；a表示初始步长向量，并有a＝(a(t),a(p),a(q))；a(t)表示大气温度t的初始步长向量，并通过式(11)获得；a(p)大气压强p的初始步长向量，并通过式(12)获得；a(q)大气比湿q的初始步长向量，并通过式(13)获得；

式(11)中，为背景场大气温度t^b的方差，为观测场大气温度t的方差；式(12)中，为背景场大气压强p^b的方差，为观测场大气压强p的方差；式(13)中，为背景场大气比湿q^b的方差，为观测场大气比湿q的方差；

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

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

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

2、本发明根据大气参数：大气温度、大气压强、大气比湿及其取值范围，构建大气参数向量的三维搜索空间，并对搜索空间进行网格化划分，将这些网格中心点作为一系列的算法搜索初始迭代点，从而避免了算法陷入局部最优，改进了算法的全局搜索能力。

3、本发明通过改进型同步扰动随机逼近算法求解代价函数，简化了代价函数的求解算法，并通过代价函数的近似梯度更新大气参数向量，不直接求解一阶和二阶导数，避免了设计观测算子的切线性算子和伴随算子。

4、本发明中改进型同步扰动随机逼近算法的扰动参数和迭代步长均与大气参数向量同维度，由大气参数向量各个值对代价函数的影响度和方差确定，且每层大气参数所对应的初始步长更新根据背景场和观测场的误差协方差矩阵自适应整定，能更快地收敛到大气参数向量的最优估计。

附图说明

图1为本发明网格化搜索空间划分示意图(N＝3)；

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

图3为本发明改进型同步扰动随机逼近算法求解代价函数的流程图。

具体实施方式

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

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

式(1)中，x为大气参数向量，即为待求解的参数，并有：x＝(t,p,q)^T，其中t表示大气温度，且t∈[t_i,t_f]，t_i表示大气温度t的取值范围的下限，t_f表示大气温度t的取值范围的上限；p表示大气压强，且p∈[p_i,p_f]，p_i表示大气压强p的取值范围的下限，p_f表示大气压强p的取值范围的上限；q表示大气比湿，且q∈[q_i,q_f]，q_i表示大气比湿q的取值范围的上限，q_f表示大气比湿q的取值范围的下限；x^b表示背景场的参数向量，并有x^b＝(t^b,p^b,q^b)^T；t^b表示背景场的大气温度，p^b表示背景场的大气压强，q^b表示背景场的大气比湿；B为背景场的误差协方差矩阵；y^o是观测场掩星资料；H(x)是观测算子；R是观测场的误差协方差矩阵；

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

大气参数求解问题就是在t,p,q的取值范围内搜索最优的参数组合x＝(t,p,q)^T，可以使代价函数J(x)达到最小。

步骤2、本发明采用同步扰动随机逼近算法求解上述问题。

假设大气参数向量x的三个元素t,p,q在空间中相互垂直，以点O(0,0,0)为原点、以大气温度t为x轴、大气压强p为y轴、大气比湿q为z轴建立关于大气参数向量x的空间直角坐标系O-xyz，并根据大气参数向量x的元素t,p,q的取值范围，构建改进型同步扰动随机逼近算法的搜索空间，搜索空间即为空间直角坐标系O-xyz上的一个立方体；

步骤3、考虑到同步扰动随机逼近算法的搜索能力受初始迭代点的选取影响较大，单次随机选取易导致算法陷入局部最优，错失对全局最优解的求解，因此提出以下算法改进策略。

算法改进策略一：搜索空间网格化划分

如图1所示，将算法搜索空间的每一维均分为N份，从而得到N³个网格；定义m,n,s分别为任意一个网格在大气温度t、大气压强p、大气比湿q三个维度上的序号，m＝1,···,N；n＝1,···,N；s＝1,···,N；并利用式(2)得到在三维空间坐标系O-xyz中空间网格的中心坐标x_m,n,s：

将这些网格中心点作为一系列的算法搜索初始迭代点。

步骤4、利用改进型同步扰动随机逼近算法在搜索空间的每个网格中求解最优解：

同步扰动随机逼近算法(SPSA)是由美国霍普金斯大学教授Spall在1987年提出的一种全新的随机逼近算法。这种逼近算法具有算式结构简单、不需要被控对象明确的数学模型、对包含噪声的测量数据有较好的处理效果等特点。在使用SPSA算法设计控制器时，不必进行针对被控系统的数学建模，只要可以进行损失函数值的测量(即使该函数估计值带有随机噪声也无妨)。最重要的是不需要重复对每个待优化参数单独扰动然后求得损失函数值进而获得一个偏导数，而是不论待优化参数的数量，仅需同时扰动所有的待优化参数，再计算两次损失函数值就可得到迭代过程中的逼近梯度。本发明中，改进型同步扰动随机逼近算法步骤如下：

步骤4.1、定义网格序号为g，g＝1,2,…,N³；并初始化g＝1；定义代价函数阈值为J_min；

步骤4.2、判断s是否等于N，否，s＝s+1,m＝1,n＝1，执行步骤4.3；是，停止迭代，得到代价函数最优解为：

x＝argmin(J(x_g)),g＝1,2,...,N³ (3)

步骤4.3、判断n是否等于N，否，n＝n+1,m＝1，执行步骤4.4；是，执行步骤4.2；

步骤4.4、判断m是否等于N，否，m＝m+1，执行步骤4.5；是，执行步骤4.3；

步骤4.5、定义迭代次数为k，并初始化k＝1；令当前网格的初始迭代点坐标x_k＝x_m,n,s，即网格中心点；令当前第g个网格的参数向量为x_g；

步骤4.6、利用式(1)获得在第k个迭代点处的大气参数向量x_k的代价函数J(x_k)，并判断k≥2且|J(x_k)-J(x_k-1)|≤J_min是否成立，若成立，则获得第k个迭代点处的大气参数向量x_k，即为当前第g个网格的最优解，并赋值给第g个网格的参数向量x_g，将g+1赋值给g，执行步骤4.4；否则，执行步骤4.7；

步骤4.7、利用式(4)求解式(1)在第k个迭代点处的大气参数向量x_k的近似梯度J′_k：

并有：

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

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

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

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

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

步骤4.8、利用式(9)获得第k+1个迭代点处的大气参数向量x_k+1：

x_k+1＝x_k-a_k.·J′ (9)

式(9)中，a_k表示第k个迭代点处的步长向量，并有：

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

变分同化代价函数中大气参数参数向量的解与背景场和观测场的误差协方差矩阵相关。若在每次更新估计值时，每一维参数的迭代步长值大小一样，则大气参数参数向量逼近最优解比较慢。

以大气压强为例，假设其最优解为 分别为第1层到第u层的大气压强最优解。背景场的大气压强向量为 分别为第1层到第u层的背景场大气压强。大气压强分析增量的向量为Δp中各元素在数值上差别很大，而且势高越大，分析增量的值越小。在确定初始步长时，若不考虑这种差异性，同步扰动随机逼近算法也会出现收敛速度慢的情况。为了加速收敛过程，我们对每次的迭代步长进行优化。

改进策略为将初始步长参数变为一个向量，各元素根据相应参数的实际情况确定，且每层大气参数所对应的初始步长更新根据背景场和观测场的误差协方差矩阵自适应整定，以达到算法能准确迅速地收敛到最优解的目的。改进策略详述如下：

算法改进策略二：初始步长参数自适应整定

式(10)中,有a＝(a(t),a(p),a(q))；a(t)表示大气温度t的初始步长向量，并通过式(11)获得；a(p)大气压强p的初始步长向量，并通过式(12)获得；a(q)大气比湿q的初始步长向量，并通过式(13)获得；初始步长向量a(t),a(p),a(q)均根据相应的方差进行调整；

式(11)中，为背景场大气温度t^b的方差，为观测场大气温度t的方差；式(12)中，为背景场大气压强p^b的方差，为观测场大气压强p的方差；式(13)中，为背景场大气比湿q^b的方差，为观测场大气比湿q的方差；

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

Claims (1)

1.一种改进型同步扰动随机逼近算法的大气参数求解方法，其特征是按如下步骤进行：

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

式(1)中，x为大气参数向量，即为待求解的参数，并有：x＝(t,p,q)^T，其中t表示大气温度，且t∈[t_i,t_f]，t_i表示大气温度t的取值范围的下限，t_f表示大气温度t的取值范围的上限；p表示大气压强，且p∈[p_i,p_f]，p_i表示大气压强p的取值范围的下限，p_f表示大气压强p的取值范围的上限；q表示大气比湿，且q∈[q_i,q_f]，q_i表示大气比湿q的取值范围的上限，q_f表示大气比湿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、假设大气参数向量x的三个元素t,p,q在空间中相互垂直，以点O(0,0,0)为原点、以大气温度t为x轴、大气压强p为y轴、大气比湿q为z轴建立关于大气参数向量x的空间直角坐标系O-xyz，并根据大气参数向量x的三个元素t,p,q的取值范围，构建改进型同步扰动随机逼近算法的搜索空间，所述搜索空间即为空间直角坐标系O-xyz上的一个立方体；

步骤3、将所述算法搜索空间的每一维均分为N份，从而得到N³个网格；定义m,n,s分别为任意一个网格在大气温度t、大气压强p、大气比湿q三个维度上的序号，m＝1,…,N；n＝1,…,N；s＝1,…,N；并利用式(2)得到在所述三维空间坐标系O-xyz中空间网格的中心坐标x_m,n,s：

步骤4、利用改进型同步扰动随机逼近算法在搜索空间的每个网格中求解最优解：

步骤4.1、定义网格序号为g，g＝1,2,…,N³；并初始化g＝1；定义代价函数阈值为J_min；

步骤4.2、判断s是否等于N，否，s＝s+1,m＝1,n＝1，执行步骤4.3；是，停止迭代，得到代价函数最优解为：

x＝argmin(J(x_g)),g＝1,2,...,N³ (3)

步骤4.3、判断n是否等于N，否，n＝n+1,m＝1，执行步骤4.4；是，执行步骤4.2；

步骤4.4、判断m是否等于N，否，m＝m+1，执行步骤4.5；是，执行步骤4.3；

步骤4.5、定义迭代次数为k，并初始化k＝1；令当前网格的初始迭代点坐标x_k＝x_m,n,s，即网格中心点；令当前第g个网格的参数向量为x_g；

步骤4.6、利用式(1)获得在第k个迭代点处的大气参数向量x_k的代价函数J(x_k)，并判断k≥2且|J(x_k)-J(x_k-1)|≤J_min是否成立，若成立，则获得第k个迭代点处的大气参数向量x_k，即为当前第g个网格的最优解，并赋值给第g个网格的参数向量x_g，将g+1赋值给g，执行步骤4.4；否则，执行步骤4.7；

步骤4.7、利用式(4)求解式(1)在第k个迭代点处的大气参数向量x_k的近似梯度J′_k：

并有：

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

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

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

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

式(4)-式(6)中，Δ_k表示第k次迭代的d维随机扰动向量，并有：

步骤4.8、利用式(9)获得第k+1个迭代点处的大气参数向量x_k+1：

x_k+1＝x_k-a_k.·J′ (9)

式(9)中，a_k表示第k个迭代点处的步长向量，并有：

式(10)中，A和α表示非负系数；a表示初始步长向量，并有a＝(a(t),a(p),a(q))；a(t)表示大气温度t的初始步长向量，并通过式(11)获得；a(p)大气压强p的初始步长向量，并通过式(12)获得；a(q)大气比湿q的初始步长向量，并通过式(13)获得；

式(11)中，为背景场大气温度t^b的方差，为观测场大气温度t的方差；式(12)中，为背景场大气压强p^b的方差，为观测场大气压强p的方差；式(13)中，为背景场大气比湿q^b的方差，为观测场大气比湿q的方差；

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