[发明专利]基于时间序列的高频地波雷达海态数据融合方法

权利要求书

1.基于时间序列的高频地波雷达海态数据融合方法，其特征在于：所述方法具体过程为：

步骤一：对高频地波雷达收到的回波进行海态反演处理，得到载频1海态数据和载频2海态数据，将载频1海态数据和载频2海态数据融合，得到融合后数据；

所述海态数据为流速、流向、风速、风向、海浪的高度；

步骤二：对融合后数据做平稳性检验，对不满足平稳性的序列做差分，直到融合后数据的所有序列满足平稳性；

步骤三：对满足平稳性的序列进行模型的拟合，得到时间序列方程；

步骤四：将步骤二融合后满足平稳性的序列减掉时间序列方程拟合出的序列后得到残差序列，观察残差序列是否符合正态分布，当残差序列满足正态分布时，时间序列方程ARIMA(p,i,q)有效；当残差序列不满足正态分布时，时间序列方程ARIMA(p,i,q)无效；

步骤五：按时间顺序将海态数据代入有效时间序列方程ARIMA(p,i,q)中，即得到下一个时刻预测的值；

步骤六：根据有效时间序列方程ARIMA(p,i,q)中的系数得到自适应卡尔曼滤波的状态转移方程与噪声方程；

步骤七：给定初值，对海态数据做自适应卡尔曼滤波，得到海态滤波数据；

步骤八：根据实际情况以一定的长度m滑动时间窗口，1≤m≤n，对待测高频地波雷达收到的回波进行海态反演处理，得到载频1海态数据和载频2海态数据，将载频1海态数据和载频2海态数据融合，得到融合后数据，重复步骤二～步骤七，得到海态滤波数据；

所述步骤一中对高频地波雷达收到的回波进行海态反演处理，得到载频1海态数据和载频2海态数据，将载频1海态数据和载频2海态数据融合，得到融合后数据；

所述n为数据时间窗口长度；

具体过程为：

载频1海态数据对应的序列为Data₁(i)，i＝1,2,3...n；

载频2海态数据对应的序列为Data₂(i)，i＝1,2,3...n；

载频1海态数据的信噪比为Snr₁(i),i＝1,2,3...n；

载频2海态数据的信噪比为Snr₂(i),i＝1,2,3,...n；

所述n为数据时间窗口长度；

将载频1海态数据和载频2海态数据融合，得到融合后数据Data(i),i＝1,2,3,...n，融合后数据Data(i)表达式为：

式中，*为乘号；

所述步骤二中对融合后数据做平稳性检验，对不满足平稳性的序列做差分，直到融合后数据的所有序列满足平稳性为止；具体步骤如下：

采用时序图方法对融合后数据Data(i)做平稳性检验；

若融合后数据的所有序列都满足平稳性，得到序列Data_stable(i),i＝1,2,3,...n，其中Data_stable(i)＝Data(i),i＝1,2,3,...n；

若融合后数据的某序列不满足平稳性，则对不满足平稳性的序列Data(i),i＝1,2,3...n做差分，差分后的序列为Data_diff(i),i＝1,2,3,...n-1，其中

Data_diff(i)＝Data(i+1)-Data(i),i＝1,2,3,..n-1；

对差分后的序列继续做平稳性检验，若不满足平稳性继续续做差分，直至满足平稳性得到序列Data_stable(i),i＝1,2,3,...m′，记录差分的次数i′，i′＝n-m′；

其中Data_stable(i)＝Data_diff(i),i＝1,2,3,...m′；

所述步骤三中对满足平稳性的序列进行模型的拟合，得到时间序列方程；具体步骤如下：

步骤三一、将满足平稳性的序列Data_stable(i),i＝1,2,3,...n拟合成时间序列方程ARIMA(p,i,q)；

式中，p为AR模型参数的个数，q为MA模型参数的个数；

步骤三二、使用AIC定阶准则对ARIMA(p,i,q)模型进行定阶，具体过程为：

AIC定阶准则如下：

选用不同的p,q，并计算相应AIC值，使AIC达到极小值时对应的系数个数为p,q取值最佳个数；

步骤三三、对确定p、q后的ARIMA(p,i,q)模型的系数取值进行估计，得到确定p、q后的ARIMA(p,i,q)模型的系数取值；

步骤三四、基于系数个数和系数取值，带入时间序列方程ARIMA(p,i,q)表达式，得到确定的时间序列方程ARIMA(p,i,q)；

所述步骤三一中时间序列方程ARIMA(p,i,q)表达式为：

x(t)＝φ₁x(t-1)+φ₂x(t-2)+...+φ_px(t-p)+a(t)+θ₁a(t-1)+θ₂a(t-2)+...θ_qa(t-q)

式中，x(t)为t时刻的海态数据，t为时间，φ₁、φ₂、φ_p为海态数据系数，a(t)、a(t-1)、a(t-2)、a(t-q)为噪声，θ₁、θ₂、θ_q为噪声系数，p为AR模型参数的个数，q为MA模型参数的个数；

所述步骤三二中

其中，Q(μ,φ₁,...φ_p,θ₁,...θ_q)为ARIMA(p,i,q)模型的剩余平方和，N为序列的长度，为序列的残差方差，μ为序列的均值；

所述步骤六中根据有效时间序列方程ARIMA(p,i,q)中的系数得到自适应卡尔曼滤波的状态转移方程与噪声方程；具体步骤如下：

通过有效时间序列方程ARIMA(p,i,0)，得到时间序列方程模型如下：

x(t)＝φ₁x(t-1)+φ₂x(t-2)+...+φ_px(t-p)+a(t)

使t＝t+1，得到

令x₁(t)＝x(t),x₂(t)＝x(t-1),...,x_p(t)＝x(t-p+1)，则有

由于x₂(t+1)＝x₁(t),x₃(t+1)＝x₂(t),...,x_p(t+1)＝x_p-1(t)

可以得到

式中，x₁(t)为t时刻x(t)取值，x_p(t)为p时刻x(t)取值；x₁(t+1)为t+1时刻海态数据值，x₂(t+1)为t时刻海态数据值，x_p(t+1)为t-p+2时刻海态数据值；

式中，Φ为状态转移矩阵，Γ为噪声矩阵；

所述步骤七中给定初值，对海态数据做自适应卡尔曼滤波，得到海态滤波数据；具体步骤如下：

(1)根据t-1时刻的状态得到t时刻的预测状态X(t|t-1)＝ΦX(t-1)

式中，X(t|t-1)为t时刻的预测状态，X(t-1)为t-1时刻的状态；

(2)给出协方差矩阵的预测P(t|t-1)＝ΦP(t-1)Φ^T+ΓQ(t)Γ^T

式中，Q(t)为系统噪声协方差矩阵，P(t|t-1)为根据t-1时刻的值预测t时刻的值，P(t-1)为t-1时刻的值，T为转置；

(3)为自适应卡尔曼滤波引入调节更新参量：

式中，b为遗忘因子；

(4)计算预测值与两个观测值的距离，选择距离近的观测值作为卡尔曼滤波的观测值；

观测误差为：

v(t)＝Z(t)-H(t)X(t|t-1)

式中，H(t)为观测矩阵，Z(t)为观测值；

(5)更新观测噪声协方差矩阵：

R(t)＝(1-d_t)R(t-1)+

d_t([I-H(t)K(t-1)]v(t)v^T(t)[I-H(t)K(t-1)]^T+H(t)P(t-1)H^T(t))

式中，I为单位矩阵，R(t-1)为t-1时刻的观测噪声协方差矩阵，K(t-1)为t-1时刻的卡尔曼增益，v(t)为观测误差；

(6)计算卡尔曼增益K(t)＝P(t|t-1)H^T(t)S-¹(t)

式中，S(t)＝H(t)P(t|t-1)H^T(t)+R(t)，K(t)为t-1时刻的卡尔曼增益，P(t|t-1)为根据t-1时刻的观测噪声协方差矩阵预测的t时刻的观测噪声协方差矩阵，S(t)为中间变量，H(t)为t时刻的观测矩阵；

(7)更新协方差矩阵P(t)＝[I-K(t)H(t)]P(t|t-1)

(8)更新k时刻的状态X(t)＝X(t|t-1)+K(t)v(t)

(9)求出海态滤波数据Y(t)；

所述遗忘因子b取值为：0.95≤b≤0.99；

所述(9)中海态滤波数据Y(t)＝X(t)H(t)。