基于改进支持向量机的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.