基于改進支持向量機的TVARMA模型辨識*

王躍鋼，鄧衛強，單 斌

(第二炮兵工程學院304實驗室，西安 710025)

TVARMA模型是ARMA模型眾多改進結構中的一種。由于引入變系數的概念，TVARMA模型徹底解決了ARMA模型對非平穩信號建模效果不佳的問題，其理論一經提出便獲得了迅速的發展，并廣泛應用于以建模與預測、目標識別、故障診斷等為目的的非平穩信號處理領域。

模型參數估計是TVARMA模型研究的重點和難點。目前，關于TVARMA模型參數估計的研究成果并不多。其中，以時間基函數展開思想為基礎的模型參數估計方法得到普遍重視和應用。文獻［1］提出了時變ARMA模型的分步估計方法;文獻［2－3］提出了時變ARMA模型參數估計的最小二乘方法;文獻［4］提出了一種基于高階統計量的模型參數估計方法;文獻［5］針對AR參數時變而MA部分參數不變的一類特殊ARMA模型，提出了一種反饋線性估計法。文獻［6］通過樣本周期圖和多點平均方法得到時變參數的函數形式，再分別采用最小二乘法和極大似然法確定其中的待定參數，從而提出一種確定廣義時變ARMA模型參數函數的方法;文獻［7］分別將AR、MA和ARMA模型的時變系數展開為小波基等一系列時頻基函數的加權和，并推導了相應參數估計方法。

支持向量機是Vapnik于1995年提出的一種機器學習算法［8］。與神經網絡相比，支持向量機由于同時最小化經驗風險和VC維的界而對預測樣本有較好的泛化能力。最小二乘支持向量機［9］在保留支持向量機結構風險最小化、小樣本等特性的前提下作了如下改進:一是將支持向量機優化模型中的損失函數設定成最小二乘損失函數;一是將不等式約束轉化為等式約束。這樣，支持向量機求解過程的二次尋優問題即轉化為線性KKT(Karush-Kuhn-Tucker)方程組的解，從而降低了求解復雜性。2004年，Rojo-Alvarez將支持向量機應用于常參數ARMA模型的辨識，取得了較好的效果［10］。

本文在前人研究成果基礎上，提出一種基于改進最小二乘支持向量機的TVARMA模型辨識方法。實驗結果證明了方法的有效性。

1 問題提出

假設非平穩時間序列x(1)，x(2)，…，x(N)滿足 TVARMA 模型存在條件［11］，則其 TVARMA(p，q)模型可表示為

其中，a1(t)，…，ap(t)，b1(t)，…，bq(t)為時變系數，{w(t)}為激勵白噪聲序列且{w(t)}～N(0，σ2)。

進一步假設模型的時變系數可展開為一系列時間基函數的線性加權和，即

fj(t)、gj(t)為時間基函數，可供選擇的時間基函數有Fourier基、多項式基、切比雪夫基等［11］。

將式(2)代入式(1)，則式(1)的變參數方程轉化為如下常方程

TVARMA模型辨識的目的就是通過某種合理的參數估計方法獲得未知參數aij和bij的估計，進而按下式估計其激勵白噪聲的方差并最終獲得完整的模型表達式。

根據下式則可以獲得模型的進化譜

2 基于改進的最小二乘支持向量機的TVARMA模型辨識

與常參數ARMA模型類似，經過基函數展開后的TVARMA模型的參數估計也是一非線性回歸過程，直接求解非常困難。然而，根據柯爾莫哥洛夫定理，基于同一觀測數據建立的ARMA模型均可等價表示為一高階AR模型或MA模型［12］。該結論對式(3)同樣成立。因此，可首先采用長自回歸法估計出TVARMA模型的白噪聲序列，從而可將非線性估計問題轉化為線性估計問題。

在以長自回歸方法獲得TVARMA模型的激勵白噪聲序列后［13］，分別將原始觀測時間序列及其激勵序列{w(t)}代入式(3)，有

進一步，令M=［X－W］，則式(6)可改寫為如下矩陣方程

對上述方程，輸入變量為矩陣M的各列所對應的變量，輸出向量為觀測數據x(1)，x(2)，…，x(N)，二者構成方程的訓練樣本集。依據結構風險最小化準則，求解上述問題的最小二乘支持向量機模型如下［14、15］:

式中，γ＞0為懲罰系數即模型平衡系數;wi為回歸函數值與實際值的誤差。

由式(8)可以看出，最小二乘支持向量機模型對未知向量θ元素所對應的不同變量賦予相同的結構風險系數1，而對不同歷史數據對應的經驗風險部分的加權系數也均為1。為此，對上述目標函數，同時引入結構風險權值矩陣Q和經驗風險權值因子vk，即考慮如下改進的最小二乘支持向量機模型

其中，矩陣 Q=diag(Q11，Q22，…，Qnn)為對角矩陣且0＜Qii≤1(i=1，2，…n;n為變量數目且n=m×p+m×q)，以控制不同變量在系統結構風險中的權重大小;vk(0≤vk≤1)用以控制不同歷史觀測樣本在經驗風險中的權重［16］。

引入拉格朗日函數，將上述約束優化問題轉化為無約束優化問題，在對偶空間可得到下式

式中，α=［α1α2… αN］為拉格朗日乘子。由KKT優化條件可得到如下等式約束條件

證明:

(2)式(10)與下式等價

對式(11)消去θ和wi可得到如下線性方程組

求解上式得到αT

代入下式可得θ的估計

由上述過程可知，對于上述改進的最小

二乘支持向量機，γ、Q和vk均為預置變量，且只需求解線性方程組就可得到α，具有簡單、快速、穩定等優點。

3 算法實驗

(1)仿真實驗

采用如下雙線性調頻信號［11］:

進行仿真實驗，采樣頻率1 000 Hz，數據長度為128，且{w(t)}～N(0，0.01)，信號時域波形如圖1所示。

采用SPWVD獲得雙線性調頻信號的功率譜等高線圖如圖2所示。

圖1 含白噪聲的雙線性調頻信號

圖2 仿真信號的SPWVD時頻譜等高線圖

首先，對矩陣MTM作奇異值分解，并將其奇異值按由小到大順序排列σ1≤σ2≤…≤σn。這時，將σn所對應的變量的結構風險權重矩陣元素設置為1，進而按σi/σn確定其他變量的結構風險權重。其次，確定經驗風險權值因子vk和γ，實驗中按指數式e－(i－N)/5確定各觀測數據對應的經驗風險權值因子，γ=500。最后，按式(12)和式(13)獲得 TVARMA模型參數估計，代入式(4)得模型激勵白噪聲序列的方差為0.033 7。此時，由式(5)可知模型參數譜等高線示意圖如圖3所示。

圖3 TVARMA模型進化譜等高線圖

模型白噪聲方差反映了模型的時域精度，TVARMA模型進化譜等高線圖與SPWVD等高線圖的匹配程度則反映模型的時頻域精度。實驗結果表明，本文方法辨識后的TVARMA模型在時域和時頻域均具有較好的精度。

(2)工程應用實驗

飛行器結構響應序列是導彈彈體結構或彈載儀器在多源外載荷作用下的結構響應，具有正態、非平穩等特點。該信號的處理結果可以為改進彈體結構設計、優化儀器艙儀器設計和配置等提供重要依據。下面以某飛行器結構響應序列對本文方法作進一步驗證，采樣長度64，采樣頻率5 120 Hz，如圖4所示。

圖4 飛行器結構響應序列

為便于對比，首先采用時頻分析方法獲得上述結構響應序列的WVD和SPWVD，如圖5和圖6所示。進而，分別采用最小二乘方法和文中方法獲得TVARMA模型的參數估計結果并分別由式(4)和式(5)估計模型白噪聲方差(0.012和0.0233)和模型進化譜(圖7和圖8)。

圖5 結構響應序列的WVD

圖6 結構響應序列的SPWVD

圖7 基于LS的TVARMA模型進化譜

圖8 基于改進LS-SVM的TVARMA進化譜

實驗結果分析。由圖5～圖7可以直觀地看出，飛行器結構響應序列是一個多分量信號，故而其WVD出現了許多由交叉項引起的虛假分量，而SPWVD和進化譜顯然受交叉項影響極小。提取三種不同方法中譜峰所對應的時間和頻率分別為(34，1480)、(34，1200)、(22，1474)和(34，1489)，提取結果表明，盡管SPWVD和基于LS的TVARMA模型進化譜可以較好地抑制交叉項干擾，但其譜峰在此過程中也發生了平滑和漂移。相比之下，基于改進的LS-SVM的TVARMA進化譜既可以取得較好的時域精度，其進化譜也保持了較高的時頻聚集性且較好地抑制了交叉項和譜峰漂移。

4 結論

本文提出了一種基于改進的最小二乘支持向量機，并將之應用于TVARMA模型參數的辨識過程。結合仿真實驗和實際工程應用實驗，給出了一種結構風險權值矩陣確定方法和經驗風險加權因子的經驗確定方法，實驗結果表明本文方法的有效性。此外，預置變量γ對參數估計結果也有較大影響，一般取γ為較大的數，約為103量級。

［1］Maiwald D，Dalle Molle J W，Bohme J F.Model Identification and Validation of Nonstationary Seismic Signals［C］//Proceedings of IEEE Signal Processing Workshop on Higher-Order Statistics，1993:319－322.

［2］Emin Yuksel M，Sadik Kara，Necmi Taspinar.Time-Dependent ARMA Modeling of Continuous Wave Ultrasonic Doppler Signals［C］//Proceedings of 8th IEEE Signal Processing Workshop on Statistical Signal and Array Processing，1996:260－263.

［3］Emin Yuksel M，Sadik Kara，Necmi Taspinar.Time-Varying Parametric Modeling and Simulation of Ultrasonic Doppler Signals［C］//TFTS’96:157－160.

［4］Donghae Kim，Paul R White.Nonstationary Parametric System Identification Using Higher-Order Statistics［C］//Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis，1998:457－460.

［5］王文華，王宏禹.非平穩信號的一種ARMA模型參數估計法［J］.信號處理，1998，14(1):33－38.

［6］傅惠民，王治華.廣義時變ARMA模型參數函數的確定方法［J］.機械強度，2004，26(6):636－641.

［7］Jachan M，Matz G，Hlawatsch F.Time-Frequency ARMA Models and Parameter Estimators for Underspread Nonsationary Random Processes［J］.IEEE Trans.on Signal Processing，2007，55(9):4366－4381.

［8］Vapnik V N，著.張學工，譯.統計學習理論的本質［M］.北京:清華大學出版社，2000.

［9］Suykens J A K，Andewalle J V.Least Squares Support Vector Machine Classifiers［J］.Neural Processing Letters，1999，9(3):293－300.

［10］Rojo-Alvarez J L，Martinez-Ramon M，Prado-Cumplido M D，et al.Support Vector Method for Robust ARMA System Identification［J］.IEEE Transactions on Signal Processing，2004，52(1):155－164.

［11］王宏禹.非平穩隨機信號分析與處理［M］.北京:國防工業出版社，1999.

［12］楊叔子.時間序列分析的工程與應用［M］.武漢:華中理工大學出版社，1991.

［13］何書元.應用時間序列分析［M］.北京:北京大學出版社，2005.

［14］丁蕾，廖同慶，陶亮.基于SVR的多傳感器數據融合處理方法［J］.傳感技術學報，2011，24(5):710－713.

［15］王曉紅，吳德會.基于SVR的傳感器Hammerstein模型辨識［J］.傳感技術學報，2007，20(5):1042－1046.

［16］付志超，程偉，徐成.基于LS-SVM的模態參數識別方法［J］.航空學報，2009，30(11):2087－2092.