[發明專利]基于時間序列的高頻地波雷達海態數據融合方法

有效

• 申請號: 202110649446.7
• 申請日: 2021-06-10
• 公開（公告）號: CN113391287B
• 公開（公告）日: 2023-09-01
• 發明（設計）人: 董英凝;鄧正鑫;鄧維波
• 申請（專利權）人: 哈爾濱工業大學
• 主分類號: G01S7/41
• 分類號: G01S7/41
• 代理公司: 哈爾濱市松花江聯合專利商標代理有限公司 23213
• 代理人: 岳昕

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)。