Volterra级数模型的非线性压缩测量辨识算法*

邱 棚，姚旭日，李鸣谦，翟光杰

(1. 中国科学院国家空间科学中心， 北京 100190； 2. 中国科学院大学， 北京 100049)

一直以来，关于非线性系统的研究吸引了来自各个领域学者的关注。特别是在工程领域中，理想的线性系统是几乎不存在的，大部分需要面对的都是非线性系统。一般认为，输出与输入不成比例的系统均是非线性的，如失真、自激振荡、分叉、混沌等。本文重点研究弱非线性系统，比如功率放大器[1]、扬声器[2]、生物过程[3]均属此类。对于此类非线性系统，其模型是很难通过理论推导得到的。因而，需要通过对输入输出信号的辨识得到系统模型。

在系统辨识中，首要考虑的问题是如何表示非线性系统。到目前为止，还没有一个通用的非线性系统表示方法。非线性系统的模型可以从以下几个方式来获得：获得系统的频域等价模型；将非线性环节表示为多个线性环节的卷积；将非线性系统展开为线性形式，再利用线性的辨识方法获得非线性系统的近似等。本文采用后一种方法，利用多项式展开的方法来表示非线性对象。根据泰勒公式可知，如果系统的非线性环节足够平滑并且没有历史记忆，那么就可将系统的输出表示为输入信号的多项式展开。但是，如果系统是有记忆的，就需要用Volterra级数来描述[4]。本文选用p阶、记忆长度为L的Volterra核作为基础，并利用多项式回归的方法得到核的参数，从而使得展开项的和尽可能地接近原函数。近年来，Volterra级数以及其他多项式展开的方式在人脸识别[5]、语音识别[6]、文字识别[7]等非线性问题上都取得了成功。

系统辨识的第二个问题是如何准确估计出模型中的参数。鉴于Volterra级数中的未知参数可以表示为输入输出的线性组合，经典的最小二乘(Least-Square，LS)方法就可辨识得到[8]。但是，此方法会碰到“维数灾难”的问题，即待辨识的参数数量将随着模型阶数的增加而呈指数增长。对于LS方法而言，就需要远大于Θ(Lp)次的数据，辨识结果才能收敛。这不仅增加了辨识的成本，还会增大计算复杂度。一个可能的方法是，将该问题视为核回归问题[9-10]。然而，Volterra级数的高阶参数是稀疏度的，也就是说待辨识参数中只有个别位置是非零的。不过，多项式的阶数、记忆长度以及非零项的位置和个数都是未知的，所以无法在辨识之前选定合适的数值。此外，在实际系统中，核回归的方法由于输入输出数据存在误差将导致过拟合，使得辨识结果变差。考虑到参数的稀疏性，有学者提出采用稀疏估计器来解决过拟合的问题，典型的如非负参数推断(Non-Negative Garrote, NNG)[11]和递归最小二乘(Recursive Least Squares,RLS)估计器。文献[18]给出了针对稀疏Volterra级数的EM-RLS算法，对参数的期望和分布进行计算，从而得到准确的辨识结果。另一个常见的稀疏估计器叫作最小绝对收缩和选择算子(Least-Absolute Shrinkage and Selection Operator, LASSO)[13]。文献[14]基于该算法，提出了针对Volterra级数的wLASSO(weight LASSO)估计器，并取得了良好的效果，而且wLASSO的性能已经在文献[15]之中分析过了。除了利用稀疏性解决“维数灾难”问题的方法以外，还有利用核的对称性以减少待辨识参数的个数[8]。另外，还可以利用正交三角分解的方法来减少辨识计算量[16]。

近年来，Donoho[17]、Candès等[18]提出的压缩感知理论是处理信号稀疏性的一种高效方法。根据采样定理可知，对任何信号的测量都要求采样频率是信号带宽的2倍以上。然而，压缩感知理论指出，对稀疏信号的测量需要关注的是稀疏度，而不是信号带宽。该定理目前已经被大量应用于有关成像的各个领域中，如高速流显微镜成像[19]、核磁共振成像[20]、合成孔径雷达[21]、毫米波雷达信号处理[22]等。Volterra级数中的一阶核可以视为一维的信号，其二阶核也可视为二维的信号。依此类推，级数的各阶核参数均可以视作高维信号，而且这些信号均是极度稀疏的。因而，利用压缩感知理论，可以极大地提高测量的效率。

1 多项式级数展开

本节主要介绍非线性系统的Volterra级数展开方式。假设一个离散非线性时不变系统的输出可表示为：

y(n)=f(x(n),x(n-1),…,x(n-l))

(1)

其中，x(n)代表系统在时刻n的输入，而f(x)则是一个非线性函数。上式表示，系统的输出与当前时刻之前的所有时刻都相关，即系统是有无限记忆的。但是，对于大部分系统而言，系统的响应都会随着时间的推移而衰减。所以，在此假设系统的记忆是有限且合理的，也就是系统的输出只与前L个时刻的输入有关。本文只考虑由死区、饱和、谐波等导致的信号失真等弱非线性系统。理论上，Volterra级数是一个无穷级数，而本文允许展开级数与原系统存在截断误差，所以系统的输出就可以利用P阶的Volterra级数来表示：

(2)

式中：v(n)包含截断误差以及测量噪声，而且该噪声是一个零均值并且独立于输入信号的白噪声；Hp[Χ1(n)]表示第p阶的Volterra核，其具体形式如下:

(3)

式中，L代表Volterra级数的记忆长度，而且各阶核的记忆长度是一致的。针对Volterra级数展开的收敛域以及使用范围已经被广泛研究，可以参考文献[8]。本文的目的是利用系统的输入输出数据，以辨识得到可以准确描述系统的各阶核参数hp(k1,…,kp)。

2 非线性压缩测量辨识

非线性压缩测量辨识算法将压缩感知理论引入对非线性系统的辨识之中，并利用测量的方法对模型参数进行观测。在测量的框架下，首先要确定被测信号的形式。在系统辨识中，被测信号包含了系统的全部特征，每一个信号都可以唯一地表示一个系统。对于线性系统来说，可以选择系统的脉冲响应序列来表示，该信号包含了线性系统的全部动态和静态特征。那么，对于非线性系统来说，是否存在一个信号是可以唯一地表达系统的呢？本文采用Volterra级数的各阶核来表示。另外，在测量的角度下，可以认为表示系统的信号是客观存在的，只不过需要特殊的观测方法才能够得到，或者是需要通过间接的方法来对其测量。本文的重点在于介绍测量在辨识非线性系统中的思想，因而只考虑二阶Volterra级数的辨识问题，而更高阶系统可以依此推广。由于忽略了高阶核，所以测量结果与原信号存在一定偏差，而这个偏差从测量的角度来看，是由于测量工具的局限性导致的。这就如同用米尺去测量毫米级的物体一般，一定会存在误差。

假设非线性系统展开为二阶Volterra级数，该系统输出可以表示为：

(4)

根据上文对Volterra级数的介绍可知，一个非线性系统的输出可以写为多维卷积的形式。对于一阶Volterra核，可以将一维卷积改写为向量相乘的形式，即：

(5)

(6)

根据式(2)可知，二阶Volterra核的输出可以写为线性二次型的形式，即：

(7)

(8)

此时，将一阶核和二阶核先后排列成一个一维列向量，那么非线性系统输出可以表示为：

Y=Y(0)+Y(1)+Y(2)=[U(0),U(1),U(2)][Hc,H1,H2]T

(9)

式中，U(0)是一个M×1的常数向量，Hc是一个1×1 的常数向量。H=[Hc,H1,H2]T是待测量信号，该信号长度N=L2+L。不过，由于二阶核是关于对角线对称的，因此可以将相同项合并，即信号长度缩减为N=(L2+3L+2)/2；U=[U(0),U(1),U(2)]是测量矩阵，该矩阵中的每一行都是一个测量向量，用于对信号H进行一次观测。因而，矩阵U的性质直接决定了测量的准确性。压缩感知理论给出了准确测量的要求，即只有当测量矩阵满足有限等距约束(Restricted Isometry Property, RIP)性质时，才能保证测量是无损的。学者已经证明一般常见的独立同分布的随机矩阵都满足RIP性质。不过，对于本文中使用的复杂随机矩阵，要想得到RIP性质的充要条件是十分困难的。因而，可以通过证明得到满足RIP性质的充分条件来得到对测量矩阵的要求。文献[23]中给出了当输入信号是符合均匀分布的随机信号时的定理如下：

该引理表明，当利用符合均匀分布的测量矩阵对信号进行测量时，如果测量数是与S2lgN在同阶大小时，恢复算法是有极大的概率从观测值中重建出被测信号的。不过，此处只是给出了均匀分布的测量准则，而且该准则只是给出了一个相对宽泛的下限。此外，如果系统输入选取的是其他随机分布，就需要重新推导RIP定理。因而，本文尝试间接地给出测量矩阵需要满足的性质。根据RIP定理的定义可知，RIP性质描述的是从测量矩阵中任意抽取任意列组成的子矩阵能否完整地保存原信号的能量，即不对原信号过分放大也不过分衰减。因而，测量矩阵的子矩阵满足RIP性质是一个必要条件。利用圆盘定理来解释，可以认为RIP性质要求测量矩阵的Grammian矩阵R=UTU的特征值在[1-δS,1+δS]范围内。概括来说，要想保证测量是准确，就需要满足以下三个条件：

②E[Ri,i]=1,i=1,…,N；

③E[Ri,j]=1,i,j=1,…,N,i≠j。

其中

(10)

3 仿真

本节将测试非线性测量辨识算法在实际非线性系统中的辨识结果与真实系统的接近程度。从测量的角度来看，对测量结果影响最大的三个因素就是测量值的个数、测量噪声以及测量向量的形式。因而，本文将测试这些参数对测量结果的影响。

本文选用一个常见的非线性系统，即Wiener-Hammerstein模型作为辨识对象。该模型是由一个线性环节、一个非线性环节和另一个线性环节串联组成，如图1所示。

图1 Simulink模型Fig.1 Simulink model

(11)

(12)

图2呈现了两种辨识算法在两种测量数下的辨识结果。其中，一阶核在图2(a)中采用了折线来表示，而二阶核则采用了三维柱状图来表示。本文定义辨识过程中采集到的系统输出个数为测量数M，又定义测量数M与信号长度N的比率为测量比率(Measurement Ratio，MR)。图2(a)是在MR为40%的情况下的辨识结果，图中上半部分展现的是Volterra级数一阶核的测量结果和原信号，图2(a)下半部分展现的是二阶核的测量结果以及原信号。其中，蓝色代表原信号，黄色代表RLS的辨识结果，红色代表非线性CMI的辨识结果。图2(b)是在MR为60%的情况下的辨识结果，图2(b)中定义同上。从这两个图中可以看到，本文提出的算法相比RLS算法辨识到的一阶核和二阶核的准确性有很大程度的提高。由于原信号是稀疏的，也就是只有靠近原点附近的项是非零的，非线性CMI方法的测量结果基本符合了这一特点，而RLS算法却将原信号能量散布在了各项上。而且，随着测量率的提高，可以明显地看到非线性CMI方法的测量准确性有了进一步的提高。

(a) MR=40%

(b) MR=60%图2 NCMI和RLS辨识得到的Volterra核Fig.2 Volterra core identified by NCMI and RLS

本文仿真了测量率从10%到200%的所有情况，并将测试结果绘制在图3中。从图中可以看到，当测量率达到66%时，非线性CMI的拟合程度已经达到94.13%，而递推最小二乘算法的测量率达到177%时，其平均拟合程度才达到同样水平。另外，当测量数大小接近于信号长度时，RLS法出现了较大的偏差。这主要是由于观测存在误差，而当测量矩阵接近方阵时，RLS算法的解空间越来越小，直到被压缩到唯一解。此时，所有观测噪声都被纳入系统模型之中，所以其偏差也就很大了。虽然，RLS算法随着测量数的继续增多，其结果的准确性也继续提高，不过最终也只能达到与非线性CMI方法同样的精度。也就是说，在该条件下，非线性CMI方法只需要90次测量就可以得到原信号，相比RLS能够节约大概150次测量值。

(a) 平均拟合程度与测量比率的关系(a) AFR of identification results under different MR

(b) 均方误差与测量比率的关系(b) MSE of Identification results under different MR图3 测量比率对NCMI和RLS的辨识影响Fig.3 Influence of NCMI and RLS identification results by MR

非线性CMI的测量结果在测量数达到90个时，就达到了最佳的测量结果。那么，如何确定这个数量，从而避免测量的浪费呢？在此，定义使得测量结果达到最佳的最小所需测量数为最小测量数。根据定理可知，最小测量数应该与稀疏度和信号长度有关，而稀疏度的影响远超过信号长度。因此，本文假设Volterra核的记忆长度为30，阶数为2，即信号总长度N为496。其中，非零项随机地分布在整个信号上。图4展示了针对不同稀疏度信号的最小测量数。根据定理可知，最小测量数应该与稀疏度的二次方成比例。但是，在线性系统中，该数值与稀疏度呈线性关系。因而，本文利用直线和抛物线分别对结果进行拟合，并且计算了两种方法拟合结果与实际值的接近程度。其中，抛物线的平均接近程度为89.3%，而直线拟合的平均接近程度为77.9%。由此可见，最小测量数与稀疏度的二次方是同阶的。换句话说，非线性CMI方法在测量过程中对信号的非零项个数更敏感。

图4 稀疏度对最小测量数的影响Fig.4 Influence of sparsity on minimum measurement number

在测量中，噪声是不可避免的，因而在包含测量噪声的条件下能否稳定地测量是至关重要的。对于测量算法来说，测量噪声的绝对大小是没有意义的，相对大小才是最重要的。因此，本文采用信噪比(Signal Noise Rate，SNR)作为指标，以测试噪声对算法的影响。图5展现了两种算法在不同噪声大小下的辨识准确度。可以看到，非线性CMI方法对噪声具有一定的鲁棒性，及时噪声达到5 dB时，信号的平均拟合程度还能保持在69.3%的水平。

图5 测量噪声对辨识结果的影响Fig.5 Influence of measurement noise on identification results

除了测量数和测量噪声以外，还有一个重要因素直接影响测量结果，那就是测量向量的形式。这是因为测量向量的形式直接决定了每次测量的方式以及测量获得的信息多少，或者说压缩感知正是凭借着精心设计的测量向量，从而实现高效采样。由前文可知，选用的输入信号是随机地从一个高斯分布中取出的。由于一阶核的测量向量是系统的输入序列，因而一阶核的测量向量是服从高斯分布的。对于二阶核来说，可将核参数分为两类，其中一类是两个不同时刻的输入序列相乘。由于两个服从高斯分布的随机变量相乘后仍然服从高斯分布，所以第一类的测量向量仍然服从高斯分布。第二类是两个相同时刻的输入序列相乘，因而观测对角线信号的测量向量就不再服从高斯分布，而是服从卡方分布。图6比较了这两种分布对同一信号的测量结果。可以看出，虽然两个分布的最佳测量结果是一致的，但是高斯分布比卡方分布的最小测量数要少得多。这主要是因为卡方分布的各项均是非负的，也就是服从该分布的测量空间减少了一半，那么肯定就需要更多的观测才能达到同样的效果。

图6 测量向量的分布对测量准确性的影响Fig.6 Influence of measurement vector distribution on accuracy

既然，Volterra二阶核对角线元素的测量服从卡方分布，那么是否可以通过修改测量向量，从而提高测量效率呢？为了解决非负的问题，一个直接的方法是，将测量向量的各项减去常数，使得其均值为零。另外，如果某一列在测量中所占比例过大，那么也会影响到其他项的测量。因而，再对整个测量矩阵的每列分别做归一化。相应地，在获得测量结果后，要将这些修改复原才能得到原信号。图7展示了没有归一化测量向量之前的测量结果，以及归一化测量向量之后的测量结果。从图中可以明显地看到，在对测量向量进行归一化后，最小测量数只需要120次测量，而此时信号的稀疏度为20。也就是说，这个最小测量数远比定理中要求的S2lgN要小得多。由此可见，通过修改测量向量的形式从而减少观测次数是可行的。

图7 归一化对测量准确性的影响Fig.7 Influence of normalization on accuracy

4 结论

本文的重要贡献在于将压缩测量算法的思想从线性领域推广到非线性领域之中，即用测量的视角观察非线性系统的辨识问题。首先，仿真也证明了本文提出的测量方法在没有噪声的情况下是可以无损地获取系统模型，即使在有噪声的条件下，也可以保证稳定的测量结果。其次，在测量的视角下，辨识算法更关注于信号本身的信息量或者说信号的稀疏度，而对信号长度不敏感。这就给解决“维数灾难”问题指出了一条新的道路。随后，本文在仿真中详细分析了影响辨识结果的三个重要因素，即测量向量、测量数和测量噪声。特别是，本文提出了一种测量向量均值归零的方法，进一步缩减了最小测量数。这也就说明，通过设计合理的测量向量，是有可能进一步提高测量效率，使得测量数控制在稀疏度的常数倍范围内。在同等测量次数下，利用本文提出的方法可以获得更多的信息，即辨识模型的记忆长度和阶数规模都可以更大。当然，这两个参数越大，Volterra模型可近似的非线性系统也就越复杂。

本文还简述了Volterra级数的物理含义。随后，本文将Volterra级数重新排列为符合压缩感知测量的形式，而且将原有测量矩阵重新排列为多个循环矩阵相连。这种方法不仅给RIP的证明提供了一个基础，还简化了辨识算法的程序设计。本文只是针对二阶以内的Volterra级数模型进行了仿真，不过，根据本文的方法可以推广到更高阶核的辨识上，该方向也将是未来的研究重点。此外，本文使用的是高斯分布组成的随机测量向量，随后将继续探索其他形式的测量向量，争取进一步地缩减最小测量数，提高测量效率。另外，本文提出的方法不仅可以用于Volterra级数，还可以用于其他非线性核的线性组合表示的模型辨识问题之中。