CN114330163A - 高频地波超视距雷达台风-电离层扰动动力学模型建模方法 - Google Patents

Abstract

本发明提供一种高频地波超视距雷达台风‑电离层扰动动力学模型建模方法，包括：同步实时获取台风‑电离层回波；根据所述台风‑电离层回波获取建立模型所需的参数；所述参数包括：台风中心位置、运动参数及台风轮廓，台风浪场和风场参数，电离层电子浓度、F2层临界频率、电离层多普勒频移、电离层高度参数；根据所述参数建立高频地波超视距雷达台风‑电离层扰动n阶多变量受控的自回归模型CAR(n)；根据所述自回归模型CAR(n)建立高频地波超视距雷达台风‑电离层扰动动力学模型。本发明将促进台风‑电离层扰动动力学特征研究，无论是对于海洋自然灾害和近地空间环境实时监测，还是对于人们认识整个日地空间环境都具有积极的重要科学意义，而且具有广泛的应用前景。

Description

高频地波超视距雷达台风-电离层扰动动力学模型建模方法

技术领域

本发明涉及突发海态与电离层扰动动力学技术领域，尤其涉及一种高频地波超视距雷达台风-电离层扰动动力学模型建模方法。

背景技术

海洋和电离层是和人类生存紧密相关的两大重要组成部分，长期以来科学家们一直在探索研究海洋与电离层之间动力学关系，从而揭开海洋与电离层之间的奥秘。尤其是一直在努力探索海洋灾害性活动对于电离层扰动响应问题，希望通过电离层扰动的观测和分析，间接提取台风、雷暴、火山爆发、地震和海啸等灾害性活动的前兆信息，以实现海洋灾害性活动的早期预警监测。这项重大交叉学科科学课题的研究首先要解决下面两个问题：一是需要构造同步实时获取海洋(台风)-电离层回波的雷达装置，以提供海洋-电离层二维信息；二是需要建立海洋(台风)-电离层之间动力学模型，以提取其动力学特征。前者作为硬件设备，现已由哈尔滨工业大学研制了同步实时获取海洋(台风)-电离层回波信号的高频地波超视距雷达(HFSWR:High Frequency Surface Wave Radar)，可以提供动力学模型建立和动力学特征研究所需要的海洋-电离层同步实时实验数据。而后者，建立海洋(台风)-电离层之间动力学模型，目前的研究基本上处于台风期间电离层扰动的观测分析阶段，尚未开展台风-电离层扰动动力学研究。

发明内容

本发明的目的在于解决上述现有技术存在的缺陷，提供一种高频地波超视距雷达台风-电离层扰动动力学模型建立方法。

高频地波超视距雷达台风-电离层扰动动力学模型建模方法，包括以下步骤：

步骤一：同步实时获取台风-电离层回波；

步骤二：根据所述台风-电离层回波获取建立模型所需的参数；所述参数包括：台风中心位置、运动参数及台风轮廓，台风浪场和风场参数，电离层电子浓度、F2层临界频率、电离层多普勒频移、电离层高度参数；

步骤三：根据所述参数建立高频地波超视距雷达台风-电离层扰动n阶多变量受控的自回归模型CAR(n)；

步骤四：根据所述自回归模型CAR(n)建立高频地波超视距雷达台风-电离层扰动动力学模型。

进一步地，如上所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，步骤二中根据所述台风-电离层回波获取台风中心位置、运动参数及台风轮廓包括：

步骤21，在台风期间，用高频地波超视距雷达录取实验数据，并通过信号处理从实验数据中同步提取海洋和电离层回波信号，从而构造海洋与电离层回波观测时间序列；

步骤22，估计台风中心位置与运动参数：在二维笛卡尔坐标系中建立台风运动模型，状态向量设为

依据高频地波超视距雷达海洋回波风场信息，用多假设交互式多模型IMM算法估计台风中心位置与台风运动参数；

步骤23，估计台风轮廓：台风轮廓看作是不规则扩展目标，分别定义目标状态变量和观测变量为

其中，N_X,k为目标个数，N_Z,k为观测个数，则系统状态方程和观测方程分别表示为

其中，i＝1,2,…,N_X,k-1，系统噪声

观测噪声

F_k为系统矩阵，H_k为观测矩阵；

观测方程在星凸模型下表示为

其中，

表示台风眼位置，Z_k,l为k时刻观测，θ_k,l为目标观测所对应角度，f(θ_k,l)为角度θ_k,l对应的半径，P(θ_k,l)＝[cos(θ_k,l) sin(θ_k,l)]^T是方向矢量，S_k,l为缩放因子，e_k,l～N(0,R)；

依据上述系统状态方程式和观测方程式，用高斯过程回归PHD估计台风轮廓。

进一步地，如上所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，步骤二中根据所述台风-电离层回波获取台风浪场包括：

对于200公里以远台风，用高频地波超视距雷达海洋回波一阶谱反演台风有效波高，即

其中，H_s为有效波高，k₀为雷达波数，α和β为拟合参数，B⁺和B^-分别是正负一阶峰的强度，s为Longuet-Higgins风向分布函数中的参数；

对于200公里以内台风，用高频地波超视距雷达海洋回波一阶谱和二阶谱联合反演台风有效波高，即

其中，ξ为拟合参数，R_w为二阶谱与一阶谱能量比值；

所述风场参数包括：台风等级与风向，其采用传统的方法进行反演。

进一步地，如上所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，步骤二中根据所述台风-电离层回波获取电离层电子浓度包括：设雷达雷达工作频率波长为λ，电离层雷达照射面积为S，由单个电子的散射截面σ_e≈10^-28m²，设照射范围内电子密度为N_e，总散射截面为σ_t为

σ_t＝SN_eσ_e

由此可反演电离层电子浓度N_e，即

其中，R为雷达与电离层的距离，即电离层虚高，P_t为雷达发射功率，P_r为雷达接收功率，G_t为雷达发射增益，G_r为接收天线增益，P_n为噪声功率，L_s为雷达系统损耗，L_p为信号在电离层中的传播损耗，λ为雷达工作频率波长，θ为发射电波仰角，ΔR为距离分辨率；

根据所述台风-电离层回波获取F₂层临界频率、电离层多普勒频移、电离层高度参数用传统方法进行反演获取得到。

进一步地，如上所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，所述步骤三包括：

步骤31、台风-电离层扰动动力学CAR(n)模型为

其中，y(k)＝(y₁(1),...,y_n(k))^T为p×1输出向量,即p个电离层扰动参数，u(k)＝(u₁(k-d),…,u_m(k-d-m))^T为q×1输入变量，即q个台风参数，e(k)＝(e₁(k),…,e_p(k))^T为p×1零均值高斯白噪声，A_i,B_i分别为p×p,p×q模型系数矩阵，d为系统时滞参数；

步骤32，设置CAR(n)模型输入与输出变量：输入变量为台风有效波高、台风轮廓、台风运动参数以及台风等级与风向，输出变量包括电离层电子浓度、F₂层临界频率、电离层多普勒频移、电离层高度；

步骤33，构造反演获得的台风与电离层扰动的观测时间序列；

步骤34，估计CAR(n)建模参数，包括模型参数、模型阶和子阶以及模型的时滞部分：依据台风和电离层扰动观测时间序列数据，由低阶开始对CAR(n)建模，然后逐次增加模型阶数，并用F检验对这些模型进行自动筛选，其估计过程由递推最小二乘法估计器、模型阶的F检验判决器、子阶和时滞的F检验判决器等三个部分组成；

所述子阶包括：自回归AR部分阶、滑动平均MA部分阶和受控部分的阶。

有益效果：

本发明将促进台风-电离层扰动动力学特征研究，无论是对于海洋自然灾害和近地空间环境实时监测，还是对于人们认识整个日地空间环境都具有积极的重要科学意义，而且具有广泛的应用前景。

附图说明

图1示出了FasterR-CNN流程图；

图2示出了HFSWR海洋和电离层回波的RD谱图；

图3示出了有效波高与Q_B关系图；

图4示出了有效波高与风速关系图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明基本依据有下面两个方面：一是台风通过重力波激发的电离层扰动机理。突发海态(如台风、海啸、雷暴等)使原有的海洋-电离层系统远离常态的平衡状态出现能量耗散，其中由台风激发的重力波作为能量传送的主要载体，可以在水平和垂直方向上进行传播，引发动量和能量的耦合，严重影响海洋和电离层之间的能量和动量通量。海平面低层大气的重力波以一定倾斜角向上传播直到电离层高度，使得大气中的中性粒子与电离层等离子体之间发生耦合，改变电离层电子密度空间分布状态，产生大范围长时间的电离层扰动效应，从而达到新的有序非平衡自组织系统。本发明所涉及的电离层扰动响应并不是由于地磁或太阳活动、太阳风等引起的扰动，而是由于受到突发海态的激励而导致电离层结构偏离其常规形态的非平衡态动力学过程，即电离层扰动状态的演变过程。二是HFSWR作为新体制雷达，能够同步实时获取海洋与电离层回波，从而为海洋(台风)-电离层扰动动力学模型建立提供海态与电离层扰动信息。哈尔滨工业大学(威海)研制成功的海洋-电离层同步获取HFSWR(威海)系统为本发明的研究提供了良好的硬件平台。

本发明模型建立具体步骤如下：

第一步：HFSWR系统开机，录取台风期间数据，并经信号处理分别提取海洋和电离层回波信号，从而构造海洋与电离层回波观测时间序列。用深度学习中的目标检测算法Faster R-CNN，从HFSWR回波RD谱中自动提取海洋和电离层回波，其算法过程如图1所示。利用Faster R-CNN算法，从HFSWR距离-多普勒RD谱中同步提取海洋和电离层回波，如图2所示。

第二步：估计台风中心位置与运动参数。为提取随台风运动电离层扰动演变过程规律，需要估计台风中心位置与台风运动参数。其中，台风中心位置，即台风质心，用于确定台风上空电离层位置。在二维笛卡尔坐标系中建立目标运动模型，状态向量设为

依据高频地波超视距雷达海洋回波风场信息，用多假设交互式多模型IMM(Interacting Multiple Model)算法估计台风中心位置与台风运动参数。其中多假设交互式多模型IMM计算流程如下：

流程1：估计值之间相互作用函数计算

多模型中子滤波器模型M_j在k-1时刻的混合状态估值由模型M_i的输出X_i(k-1,k-1)、相应的模型概率w_i(n)以及过渡概率ρ_ij来计算，即

在k-1时刻，子滤波器模型M_j的状态误差协方差阵为

其中，P_i(k-1,k-1)为时刻k-1第i个子滤波器的协方差估值。

流程2：滤波计算根据转换量测卡尔曼滤波算法方程，更新每个子滤波器的状态X_j(k,k)和P_j(k,k)。

流程3：模型概率计算

由残差v_j和协方差更新量B_j计算子模型M_j的似然值，即

其中，m为测量向量的维数。

由贝叶斯定理得到模型概率w_j(k)，即

其中，ρ(z(k)|z^k-1)是归一化常数，它确保所有子模型概率之和等于1。

流程4：联合估计计算

由所有子模型的状态估计加权合获得最终的状态估计值，即

估计协方差阵为

第三步，估计台风轮廓(形状)：台风轮廓(形状)看作不规则扩展目标(Extendedtarget,ET)分别定义目标状态变量和观测变量为

其中，N_X,k为目标个数，N_Z,k为观测个数，则系统状态方程和观测方程表示为

其中，i＝1,2,…,N_X,k-1，系统噪声

观测噪声

F_k为系统矩阵，H_k为观测矩阵。

观测方程在星凸模型下表示为

其中，

表示台风眼位置，Z_k,l为k时刻观测，θ_k,l为目标观测所对应角度，f(θ_k,l)为角度θ_k,l对应的半径，P(θ_k,l)＝[cos(θ_k,l) sin(θ_k,l)]^T是方向矢量，S_k,l为缩放因子，e_k,l～N(0,R)。依据上述系统状态方程式(8)和观测方程式(10)，用高斯过程回归PHD(GaussianProcess Regression Probability Hypothesis Density,GPRPHD)估计台风轮廓(形状)。

第四步：用HFSWR海洋回波一阶谱和二阶谱反演台风有效波高。对于200公里以远台风，用一阶谱反演有效波高，即

其中，H_s为有效波高，k₀为雷达波数，α和β为拟合参数，B⁺和B^-分别是正负一阶峰的强度，s为Longuet-Higgins风向分布函数中的参数。图3为依据P-M谱得到的有效波高与Q_B的关系图，由图3可知，基于一阶谱的有效波高反演方法可精确反演风眼内有效波高，而在风眼周围区域则会进入饱和区，即在风眼周围区域反演的有效波高很高，一般10米以上，依据该特点可进行风眼位置和大小的确定。风眼周围区域的有效波高和风速可由径向流速近似获取。图4表示浪高与风速关系。

对于200公里以内台风，用一阶谱和二阶谱联合反演台风有效波高，即

其中，ξ为拟合参数，R_w为二阶谱与一阶谱能量比值。

用传统方法反演台风风场参数，包括台风等级与风向。

第五步，反演电离层电子浓度、F₂层临界频率、电离层多普勒频移、电离层高度参数，包括以下步骤：

反演电离层电子浓度：设雷达雷达工作频率波长为λ，电离层雷达照射面积为S，由单个电子的散射截面σ_e≈10^-28m²，设照射范围内电子密度为N_e，总散射截面为σ_t为

σ_t＝SN_eσ_e (14)

由此可反演电离层电子浓度N_e，即

其中，R为雷达与电离层的距离，即电离层虚高，P_t为雷达发射功率，P_r为雷达接收功率，G_t为雷达发射增益，G_r为接收天线增益，P_n为噪声功率，L_s为雷达系统损耗，L_p为信号在电离层中的传播损耗，λ为雷达工作频率波长，θ为发射电波仰角，ΔR为距离分辨率。

用传统方法反演F2层电离层临界频率、电离层多普勒频移、电离层高度参数。

第六步，选择HFSWR台风-电离层扰动动力学模型CAR(n)。HFSWR海洋-电离层之间动力学要素关系研究离不开系统动力学模型。在数学领域，微分方程和差分方程是描述动态系统的重要方法。由于任何一个白噪声激发的n阶随机微分方程(即动力学模型)对应的采样模型是n阶自回归滑动平均模型ARMA(n,n-1)(Auto-Regressive and Moving AverageModel)，为此本发明把反演获得的海洋和电离层参数作为观测序列，建立动力学要素时间序列模型来研究动力学关系。海洋和电离层之间动力学要素关系属于多变量关系问题，因此需要建立动态系统受控的自回归滑动平均模型CARMA(Controlled Auto-RegressiveMoving Average Model)。任何CARMA模型均可以充分高阶CAR(n)模型逼近到任意精度。为此，本发明采用比CARMA模型更为简单的CAR(n)模型。采用CAR(n)模型的最大优点是用递推最小二乘法可以获得模型参数的一致估计，而且计算简单。

第七步：建立台风-电离层扰动n阶多变量受控的自回归模型CAR(n)。台风-电离层扰动动力学CAR(n)模型为

其中，y(k)＝(y₁(1),...,y_n(k))^T为p×1输出向量,即p个电离层扰动参数，u(k)＝(u₁(k-d),...,u_m(k-d-m))^T为q×1输入变量，即q个台风参数，e(k)＝(e₁(k),...,e_p(k))^T为p×1零均值高斯白噪声，A_i,B_i分别为p×p,p×q模型系数矩阵，d为系统时滞参数。

第八步，确定模型CAR(n)输入与输出变量。输入变量为台风浪场与风场参数，包括台风有效波高、台风轮廓、台风运动参数以及台风等级与风向，输出变量为电离层扰动参数，包括电离层电子浓度、F2层临界频率、电离层多普勒频移、电离层高度；

第九步，构造反演获得的台风与电离层扰动的观测时间序列；

第十步，估计CAR(n)建模参数，包括模型参数、模型阶和子阶(自回归AR部分阶、滑动平均MA部分阶和受控部分的阶)以及模型的时滞部分：依据台风-电离层扰动观测时间序列数据，由低阶开始对CAR(n)建模，然后逐次增加模型阶数，并用F检验对这些模型进行自动筛选，其估计过程由递推最小二乘法估计器、模型阶的F检验判决器、子阶和时滞的F检验判决器等三个部分组成。建立HFSWR台风-电离层扰动模型的关键之一是模型时滞部分的确定，而且时滞部分阶的确定过程十分复杂，其复杂性主要表现在以下三个方面：一是重力波传播机理的复杂性。台风激励的电离层扰动是通过重力波复杂的传播和耗散动力行为响应的，因此系统时滞不但和重力波非线性特性有关，而且还和重力波外部环境有关，其传播机理十分复杂；二是系统时滞的时变性。不同的突发海态将产生不同的重力波(指能量大小和传播路径)，因此不同的重力波所具有的不同特性将带来电离层扰动响应的随机时滞，并不像传统的时间序列模型具有固定的时滞，这给系统建模的时滞部分确定带来了新的问题；三是时间对齐复杂性。HFSWR台风和电离层数据录取离散随机序列的时间周期不完全相同，即海洋回波提取信号处理积累周期和电离层回波提取信号处理积累周期不相同，因此为建立时间序列模型首先要进行时间对齐处理，使海洋和电离层反演参数时间序列转化为具有同一个周期的、时间状态为离散的随机变量，由此构成台风和电离层扰动的同步离散随机过程。本发明采用时滞不均匀分段的方法(通过实验确定时滞变化范围)，在不同的分割区间分别构造时间序列模型。

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

Claims (6)

1.高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，包括以下步骤：

步骤一：同步实时获取台风-电离层回波；

步骤二：根据所述台风-电离层回波获取建立模型所需的参数；所述参数包括：台风中心位置、运动参数及台风轮廓，台风浪场和风场参数，电离层电子浓度、F2层临界频率、电离层多普勒频移、电离层高度参数；

步骤三：根据所述参数建立高频地波超视距雷达台风-电离层扰动n阶多变量受控的自回归模型CAR(n)；

步骤四：根据所述自回归模型CAR(n)建立高频地波超视距雷达台风-电离层扰动动力学模型。

2.根据权利要求1所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，步骤二中根据所述台风-电离层回波获取台风中心位置、运动参数及台风轮廓包括：

在台风期间，用高频地波超视距雷达录取实验数据，并通过信号处理从实验数据中同步提取海洋和电离层回波信号，从而构造海洋与电离层回波观测时间序列；

估计台风中心位置与运动参数：在二维笛卡尔坐标系中建立台风运动模型，状态向量设为

依据高频地波超视距雷达海洋回波风场信息，用多假设交互式多模型IMM算法估计台风中心位置与台风运动参数；

估计台风轮廓：台风轮廓看作是不规则扩展目标，分别定义目标状态变量和观测变量为

其中，N_X,k为目标个数，N_Z,k为观测个数，则系统状态方程和观测方程分别表示为

其中，i＝1,2,…,N_X,k-1，系统噪声

观测噪声

F_k为系统矩阵，H_k为观测矩阵；

观测方程在星凸模型下表示为

其中，

表示台风眼位置，Z_k,l为k时刻观测，θ_k,l为目标观测所对应角度，f(θ_k,l)为角度θ_k,l对应的半径，P(θ_k,l)＝[cos(θ_k,l) sin(θ_k,l)]^T是方向矢量，S_k,l为缩放因子，e_k,l～N(0,R)；

依据上述系统状态方程式和观测方程式，用高斯过程回归PHD估计台风轮廓。

3.根据权利要求1所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，步骤二中根据所述台风-电离层回波获取台风浪场包括：

对于200公里以远台风，用高频地波超视距雷达海洋回波一阶谱反演台风有效波高，即

其中，H_s为有效波高，k₀为雷达波数，α和β为拟合参数，B⁺和B^-分别是正负一阶峰的强度，s为Longuet-Higgins风向分布函数中的参数；

对于200公里以内台风，用高频地波超视距雷达海洋回波一阶谱和二阶谱联合反演台风有效波高，即

其中，ξ为拟合参数，R_w为二阶谱与一阶谱能量比值；

所述风场参数包括：台风等级与风向，其采用传统的方法进行反演。

4.根据权利要求1所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，步骤二中根据所述台风-电离层回波获取电离层电子浓度包括：设雷达雷达工作频率波长为λ，电离层雷达照射面积为S，由单个电子的散射截面σ_e≈10^-28m²，设照射范围内电子密度为N_e，总散射截面为σ_t为

σ_t＝SN_eσ_e

由此可反演电离层电子浓度N_e，即

其中，R为雷达与电离层的距离，即电离层虚高，P_t为雷达发射功率，P_r为雷达接收功率，G_t为雷达发射增益，G_r为接收天线增益，P_n为噪声功率，L_s为雷达系统损耗，L_p为信号在电离层中的传播损耗，λ为雷达工作频率波长，θ为发射电波仰角，ΔR为距离分辨率；

根据所述台风-电离层回波获取F₂层临界频率、电离层多普勒频移、电离层高度参数用传统方法进行反演获取得到。

5.根据权利要求3所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，所述步骤三包括：

步骤31；台风-电离层扰动动力学CAR(n)模型为

其中，y(k)＝(y₁(1),...,y_n(k))^T为p×1输出向量,即p个电离层扰动参数，u(k)＝(u₁(k-d),...,u_m(k-d-m))^T为q×1输入变量，即q个台风参数，e(k)＝(e₁(k),...,e_p(k))^T为p×1零均值高斯白噪声，A_i,B_i分别为p×p,p×q模型系数矩阵，d为系统时滞参数；

步骤32；设置CAR(n)模型输入与输出变量：输入变量为台风有效波高、台风轮廓、台风运动参数以及台风等级与风向，输出变量包括电离层电子浓度、F₂层临界频率、电离层多普勒频移、电离层高度；

步骤33；构造反演获得的台风与电离层扰动的观测时间序列；

步骤34，估计CAR(n)建模参数，包括模型参数、模型阶和子阶以及模型的时滞部分：依据台风和电离层扰动观测时间序列数据，由低阶开始对CAR(n)建模，然后逐次增加模型阶数，并用F检验对这些模型进行自动筛选，其估计过程由递推最小二乘法估计器、模型阶的F检验判决器、子阶和时滞的F检验判决器等三个部分组成；

所述子阶包括：自回归AR部分阶、滑动平均MA部分阶和受控部分的阶。

6.根据权利要求2所述的高频地波超视距雷达台风-电离层扰动动力学模型建模方法，其特征在于，所述多假设交互式多模型IMM计算方法包括：

步骤1：估计值之间相互作用函数计算

多模型中子滤波器模型M_j在k-1时刻的混合状态估值由模型M_i的输出X_i(k-1,k-1)、相应的模型概率w_i(n)以及过渡概率ρ_ij来计算，即

在k-1时刻，子滤波器模型M_j的状态误差协方差阵为

其中，P_i(k-1,k-1)为时刻k-1第i个子滤波器的协方差估值；

步骤2：滤波计算

根据转换量测卡尔曼滤波算法方程，更新每个子滤波器的状态X_j(k,k)和P_j(k,k)；

步骤3：模型概率计算

由残差v_j和协方差更新量B_j计算子模型M_j的似然值，即

其中，m为测量向量的维数；

由贝叶斯定理得到模型概率w_j(k)，即

其中，ρ(z(k)|z^k-1)是归一化常数，它确保所有子模型概率之和等于1；

步骤4：联合估计计算

由所有子模型的状态估计加权合获得最终的状态估计值，即

估计协方差阵为

Images