CN107688554B - 基于自适应傅里叶分解的频域辨识方法 - Google Patents

Abstract

本发明公开了一种基于自适应傅里叶分解的频域辨识方法，主要包含两阶段，第一阶段是估计频率响应：首先采集系统的输入和输出数据，然后对数据进行预处理和谱分析，得到系统的频率响应估计值；第二阶段是辨识传递函数模型：主要运用AFD算法，根据第一阶段得到的频率响应估计值将模型表示为加权Blaschke乘积，从而得到离散传递函数。本发明的频域辨识方法通用性强，计算速度快，辨识精度高。

Description

基于自适应傅里叶分解的频域辨识方法

技术领域

本发明属于系统辨识领域，特别是一种基于自适应傅里叶分解的频域辨识方法。

背景技术

建模的方法主要有两种：机理建模与系统辨识方法建模。理论上，基于基本原理建立的数学模型物理意义明确，能够充分反映系统的内在运动机理。然而对于部分复杂的系统，无法对其内部进行精确分析，而系统辨识方法建模是利用系统的输入输出数据，通过辨识算法和误差准则去得到一个具有等价输入输出特性的模型，其优势在于不需要进行详细受力分析和精密测量。

系统辨识可以分为时域辨识和频域辨识。与时域辨识相比，频域辨识无需对模型的结构做任何假设即可得出非参数模型，可以提供系统的初步信息，直观了解系统的特性。通常，在外界干扰比较大的情况下如果采用时域辨识则不能通过一个简单的滤波进行直接处理，还要通过复杂的计算在原有模型基础上再建立相对繁琐的噪声模型。而在频域辨识中采用合适的方式，可以去除指定频段内的噪声，使得最终辨识结果集中于所关注的频段，使得辨识结果准确性得到很大提高。

传统的辨识方法经过多年发展已经变得比较成熟，但是对于实际应用中的一些复杂系统，难以得到精确的模型或令人满意的结果，其主要存在以下局限和不足：1)部分辨识方法存在着不能同时确定系统的结构与参数的问题，需根据假设与尝试来分别进行辨识；2)部分辨识方法要求输入信号已知，且必须具有较丰富的变化，这一条件在许多普通闭环控制系统是可以满足的，但在某些系统中对输入有严格的限制，因此这些方法不便直接应用；3)部分辨识方法运行过程复杂，计算时间耗费大、效率低，而且可能得到的是局部最优解，因此辨识精度相对较低。

发明内容

本发明的目的在于提供一种基于自适应傅里叶分解的频域辨识方法，能够更加精确估计频率响应以及提高辨识模型的精度。

实现本发明目的的技术解决方案为：一种基于自适应傅里叶分解的频域辨识方法，包括以下步骤：

步骤1、采集实验数据：确定采样率与数据长度、选择试验信号、选择采集输入输出数据的方式；

步骤2、数据预处理：对输入输出数据剔除野值，去除均值和趋势项，进行数字滤波以及数据串接；

步骤3、谱分析：采用层叠分窗法得到输入功率谱密度、输出功率谱密度、互谱密度；

步骤4、频率响应辨识：计算频率响应估计值；

步骤5、传递函数辨识：将频率响应估计值作为AFD算法的输入，选择需要拟合的阶数，经过算法的迭代与分解输出一系列参数，得到离散传递函数；

步骤6、传递函数转化：将离散传递函数经过双线性变换得出连续传递函数；

步骤7、使用评估拟合精度高低的代价函数和反映模型预测能力和可靠性的时域验证，对步骤6中连续传递函数进行评价。

本发明与现有技术相比，其显著优点为：

(1)本发明使用层叠分窗法，通过对多段数据上得到的粗略谱估计求平均，可以大幅减小谱估计中的随机噪声，进而得到平滑的谱估计结果；而且在对频率响应进行估计时，同时考虑了输入输出测量噪声的影响，使之更接近于真实值；

(2)本发明运用了一种高效的有理逼近算法-AFD，由于其自适应特性和在能量角度的快速收敛性，辨识得出的模型具有更高的精度，并可以满足后续对系统设计控制器的需求；

(3)本发明的辨识方法通用性强，计算效率快，辨识精度高。

附图说明

图1是本发明基于自适应傅里叶分解的频域辨识方法流程图。

图2是本发明中的频率响应辨识的总流程图。

图3(a)和图3(b)是输入时间历程数据和输出时间历程数据的图像。

图4(a)和图4(b)是本发明中根据四种频率响应估计公式计算得出的结果的幅频特性曲线和相频特性曲线波特图。

图5(a)和图5(b)是本发明中使用AFD算法拟合得出的传递函数与频率响应估计值的对比图，图5(a)为幅频特性曲线图，图5(b)为相频特性曲线图。

图6(a)和图6(b)是本发明中使用CIFER软件得出的连续传递函数与频率响应数据的对比图，其中图6(a)为幅频特性曲线图，图6(b)为相频特性曲线图。

图7是本发明中的AFD辨识方法与CIFER软件进行比较的方案示意图。

图8是本发明中使用测试信号1来进行两种辨识方法时域验证的结果对比图。

具体实施方式

结合图1、图2，一种基于自适应傅里叶分解的频域辨识方法，包括以下步骤：

步骤1、采集实验数据：确定采样率与数据长度、选择试验信号、选择采集输入输出数据的方式；

步骤2、数据预处理：对输入输出数据剔除野值，去除均值和趋势项，进行数字滤波以及数据串接；

步骤3、谱分析：采用层叠分窗法得到输入功率谱密度、输出功率谱密度、互谱密度；

步骤4、频率响应辨识：计算频率响应估计值；

步骤5、传递函数辨识：将频率响应估计值作为自适应傅里叶分解算法(AFD)的输入，选择需要拟合的阶数，经过算法的迭代与分解输出一系列参数，得到离散传递函数；

步骤6、传递函数转化：将离散传递函数经过双线性变换得出连续传递函数；

步骤7、使用评估拟合精度高低的代价函数和反映模型预测能力和可靠性的时域验证，对步骤6中连续传递函数进行评价。

进一步的，步骤2中的剔除野值运用的是拉依达准则，首先计算测量数据列的均值以及标准偏差，以置信概率99.7％为标准，经判断若序列值大于三倍标准偏差，则属于粗大误差然后剔除，若小于则予以保留；

具体公式为：

式中，x_i是一组数据序列，i＝1,2…n，n为数据序列个数，

是该序列的均值，σ是序列的标准偏差。

进一步的，步骤2中的去除趋势项具体是指先用最小二乘拟合出数据的趋势线，然后减去该趋势线；其中的数字滤波是指先根据实验对象的要求设计IIR滤波器或者FIR滤波器，然后将数据通过滤波器进行滤波。

进一步的，步骤3中的层叠分窗法是一种谱估计方法，其具体方法包括以下步骤：

步骤3-1、计算时域数据的傅里叶变换；

其中，x(n),y(n)是时域数据记录，X(f_k),Y(f_k)是第k个频率点的傅里叶系数，w(n)是离散的窗函数，Δt是采样时间间隔，

是频率分辨率，N是离散频率点的个数；

步骤3-2、计算谱函数的粗略估计值；

其中，

为第k个窗口计算得出的输入自谱，

为第k个窗口计算得出的输出自谱，

为第k个窗口计算得出的互谱，T_win为窗口长度；

步骤3-3、计算谱函数的平滑估计；

输入自谱的平滑估计：

输出自谱的平滑估计：

互谱的平滑估计：

其中，U是窗函数功率恢复系数，用于修正谱密度的幅值；m是窗口片段的数量；

根据得到的谱函数，引出相干函数估计值

它在各频率点上的定义为：

相干函数表示了在频率点f，由输入自谱线性地引起的输出自谱的部分占整体输出自谱的比例，其取值范围是0～1；若相干函数值满足大于等于0.6且不振荡，则认为辨识精度可接受。

进一步的，步骤4中的频率响应估计公式需要根据具体情况来选择，具体分为以下几种：

(1)H1估计

若只考虑输出测量噪声，由H1估计得到的频率响应值能保证对所有非相关的输出测量噪声情况是无偏的,但结果会导致估计值相对于真实响应偏小；

(2)H2估计

如果假定噪声仅存在于输入测量中，而输出的测量中不含噪声，则H2估计可获得无偏的频率响应估计值。但在实际中，会导致估计值相对于真实响应偏大；

在一般的实际测量中，输入输出噪声是同时存在的，考虑到H1估计和H2估计分别为欠估计和过估计，为了减小误差影响，提出了以下折中算法；

(3)H3估计

(4)H4估计

进一步的，步骤5中的AFD算法是一种有理逼近算法，其具体方法包括以下步骤：

步骤5-1、根据先验知识，选择需要拟合的传递函数的阶数n′，并把步骤4得到的频率响应估计值取共轭后，作为算法的初始值，

步骤5-2、根据极大选择原理，计算出系数a_j，j＝1,2…n′；

其中涉及到的公式有：

式中，

表示单位圆盘，e_{a}(z)是算法中定义的一种字典，<>符号表示内积；

步骤5-3、根据已有的系数a_q和g_q，更新下一项g_q+1，并计算部分分量的和；

其中涉及到的公式有：

式中，

是算法中选用的一种有理正交系，也叫Takenaka-Malmquist system；

步骤5-4、将g_q+1返回到步骤5-2，求出下一个系数a_q+1，如此循环，直到求出n′个系数a_j，最后得到离散传递函数公式；

式中，m_j表示第j项前面的系数；

步骤5-5、最后经过转换

得到最终的离散传递函数，公式即为：

进一步的，步骤7中的时域验证是一种通过比较在同样激励下模型的输出和测量的系统输出来实现对模型的评价；具体使用了均方误差(Mean Squared Error,MSE)这一指标，公式如下：

式中，x_i为系统真实的输出序列，

为模型预测的输出序列；i＝1,2,…,n；

而代价函数是用来描述所得到的传递函数的拟合精度，当其值小于100时，大体上反映了飞行动态建模可接受的精度水平，公式如下：

其中，||符号表示每个频率ω处的幅值(dB)，∠符号表示每个频率ω处的相位(°)，

表示频率响应估计值，T表示传递函数模型，n_ω表示频率采样点的数量，ω₁和

表示拟合的起始频率值和结束频率值，W_g和W_p分别是幅值和相位平方误差的相对权重，通常取W_g＝1.0，W_p＝0.01745；

W_γ是加权函数，它依赖于每个频率ω处的相干函数值r_xy，公式如下：

下面结合附图和实施例对本发明进行详细说明。

实施例

四旋翼无人飞行器具有造价低、飞行特性良好、独特的位姿控制方式等多方面优点。然而，作为一个多输入多输出、非线性时变系统，四旋翼飞行器的输入量和输出量之间存在耦合，再加上系统本身的不确定性和外部干扰的引入，使得系统的控制问题变得十分复杂。若要实现对四旋翼飞行器的精确控制，必须先建立精准的模型，而系统辨识就是利用系统的输入输出数据，通过辨识算法和误差准则去得到一个可以描述系统等价输入输出特性的模型。

四旋翼飞行器主要包括四个输入通道，分别是横滚通道、俯仰通道、油门通道、偏航通道。而输出响应包括轴向的速度分量、角速度分量(横滚角速率、俯仰角速率、偏航角速率)、机体欧拉角(横滚角、俯仰角、偏航角)等。而本实施例将需要辨识的对象简化为四旋翼飞行器悬停状态点附近的线性化模型，然后运用本发明提出的频域辨识方法来辨识。具体辨识的是以遥控器横滚通道c1为输入，飞行器横滚角roll为输出的单输入单输出传递函数模型。

结合图1、图2，下面将按照步骤进行详细说明：

步骤1，首先进行试验获取四旋翼飞行器的输入输出数据。采集和记录飞行数据的平台是四旋翼飞行器的飞控，飞控中含有惯性测量单元，包括三轴加速度计、三轴陀螺和磁力计等。频域辨识的推荐输入信号是扫频输入，因为扫频信号可在感兴趣的频率范围内提供很均匀的谱激励，从而利于得到好的辨识结果；另外对系统动态中先验知识的不确定性也具有鲁棒性。

本实验扫频输入的频率覆盖范围在0.1Hz到3Hz，采样率设置为10Hz。具体操作如下：四旋翼飞行器首先在高度1.5-2米处实现悬停，同时保持其他通道相对稳定，悬停3s后在横滚通道施加扫频输入，单次采集的数据长度约100s左右，最后回到平衡点悬停一段时间后降落；重复采集5组数据，从通道记录的数据中挑选几组较好的数据作为待辨识数据，并将输入信号做归一化处理。最后得到的一组数据图像如图3(a)和图3(b)所示。

步骤2，在Matlab软件中，编写数据预处理程序，对采集得到的数据进行初步处理。首先运用根据拉依达准则编写的程序，把由于测试仪器工作不稳定而使信号发生偶然跳动或因外界异常干扰的影响导致测量数据中出现的不合理跳点去除；接着运用detrend函数，把数据测量中引入的不必要的直流分量和由于传感器测量时产生的零点漂移、低频性能的不稳定以及周围的环境干扰等导致的趋势项成分去除；然后为了增加信息量，将几组处理后的实验数据串接成一个数据记录；最后根据需要关注的频率范围，选择合适的截止频率和窗函数种类设计FIR低通滤波器，来滤除信号中的高频噪声。

步骤3，对预处理后的信号进行谱分析。Welch法也叫层叠分窗法，是一种在实际谱分析中用于大幅减小谱估计中随机噪声的关键技术。其基本原理是：将原始的时间历程数据记录分割为短的、相互重叠的片段，每个数据片段用窗口函数进行加权，然后对每一段进行傅里叶变换得到谱密度，最后取平均。

在分窗法中，有几个关键点：1、数据重叠在很大程度上能抵消窗口边缘缩减所造成的误差增大，增加重叠比可以减小随机误差，采用50％到80％的重叠比。2、所用的窗函数一般选择汉宁窗，因为渐变的拱形窗口可以减小旁瓣泄露引起的谱估计误差。其中窗口大小的选择需要特别注意，增大窗口尺寸，频率分辨率会提高，能改善低频段辨识结果，但同时会降低高频段的辨识精度，也会减少数据段的数量，从而使随机误差增加；选择小的窗口则正好与此相反。

步骤4，根据求得的谱函数来计算频率响应估计值。在实际辨识应用中，输入输出的测量值中通常包含噪声。输入端噪声来源，如对驾驶杆位移，作动器位移等测量时引入的扰动；输出端噪声来源，如飞行器的动态响应测量量，如陀螺仪、加速度计，受到来自机械、电子、大气扰动等方面的干扰。

选择何种频率响应估计公式(有

)更适合于当前应用，取决于测量噪声的主要来源有关。在本实施例中，考虑到输入输出噪声都存在，所以选择了公式

得出的频率响应估计值一般以波特图的形式直观显示，如图4(a)和图4(b)所示，分别为四种公式对应的幅频特性曲线和相频特性曲线。

步骤5，运用AFD算法对得到的频率响应估计值进行拟合逼近，得出离散传递函数。首先将频率响应估计值作为算法的初值，然后根据先验知识选择需要拟合的阶数，本实施例选择3阶或者4阶皆可。AFD算法拟合的结果同样以波特图来显示，如图5(a)和图5(b)所示，同时列举得出的4阶离散传递函数如下：

步骤6，由于后续对比的需要，在Matlab中，将离散传递函数经过双线性变换得到连续传递函数。列举得出的4阶连续传递函数如下：

步骤7，得出传递函数模型后，需要进行时域验证来对模型进行评价。模型的辨识是用频域法完成的，但是对于模型进行检验需要在时域中进行，以验证其是否具有较好的预测能力和对输入形式的鲁棒性。具体方法：先用某种测试输入信号激励系统，再测量系统在这种信号作用下的真实输出，最后通过比较在同样激励下模型的输出与测量的真实输出来实现。本实施例将双极方波作为测试信号，因为这种信号可以很好地显示关键的飞行器动态以及模型性能，而且该机动的瞬态特性是常见机动飞行中有代表性的瞬态阶段，是一种确保模型不被过调的很好的检测信号。

最后，为了验证本发明的辨识精度高，在频域辨识领域选取了一款专门针对飞行器的软件CIFER来进行比较。CIFER软件是由艾姆斯研究中心开发的，在众多的系统辨识软件包中，它是第一个完成对端到端的频率响应辨识的集成软件包，实践证明它是解决旋翼飞行器系统辨识中困难问题的有效工具。CIFER在它的用户界面里引入了较多非常实用的有关于飞行器系统辨识的准则，并将许多高效的辨识功能进行了高度集成，目前已被广泛应用于飞行器项目中。

选取由CIFER软件辨识得出的4阶连续传递函数模型，图6(a)和图6(b)用波特图的形式展示了其拟合的效果，而传递函数公式如下所示：

两种方法的比较方案如图7所示，图中的“FRESPID模块”指CIFER中对时域数据进行处理，最终得到频率响应数据的一个模块；而“NAVFIT模块”指CIFER中通过拟合频率响应数据，然后得到零极点式传递函数模型的一个模块。

在得到本发明辨识方法和CIFER软件计算出的传递函数模型后，分别选取多组试验信号对其进行时域验证，然后计算均方误差，结果列表如下：

表1 4阶情况下两种方法的均方误差比较结果

表2 3阶情况下两种方法的均方误差比较结果

从两个表格的数据上来看，本发明方法的均方误差指标比CIFER软件的更小，表明本方法辨识得出的模型精度更高。同时，仅选取4阶模型下的时域验证结果，如图8所示，两种方法所得模型的输出很接近且与真实输出相差不大。

综上所述，本发明首先通过一系列预处理将辨识所需的数据变得优质，接着使用层叠分窗法减小了谱估计的随机误差，然后根据实际情况选择频率响应估计公式来减小测量噪声的影响，得到更加接近真实的系统频率响应。之后运用高效的AFD算法来逼近频率响应估计值，得到离散传递函数模型，最后通过时域验证等方式进行评价，还与专业的CIFER软件进行比较，证明了本发明辨识得出的模型具有很高的精度，可以在一定程度上真实反映系统的输入输出情况，具有可靠的实用性。

Claims (5)

1.一种基于自适应傅里叶分解的频域辨识方法，其特征如下，包括以下步骤：

步骤1、采集实验数据：确定采样率与数据长度、选择试验信号、选择采集输入输出数据的方式；

步骤2、数据预处理：对输入输出数据剔除野值，去除均值和趋势项，进行数字滤波以及数据串接；

步骤3、谱分析：采用层叠分窗法得到输入功率谱密度、输出功率谱密度、互谱密度；

步骤4、频率响应辨识：计算频率响应估计值；

步骤5、传递函数辨识：将频率响应估计值作为自适应傅里叶分解算法的输入，选择需要拟合的阶数，经过算法的迭代与分解输出一系列参数，得到离散传递函数；自适应傅里叶分解算法是一种有理逼近算法，其具体方法包括以下步骤：

步骤5-1、根据先验知识，选择需要拟合的传递函数的阶数n′，并把步骤4得到的频率响应估计值取共轭后，作为算法的初始值，

步骤5-2、根据极大选择原理，计算出系数a_j，j＝1,2…n′；

其中涉及到的公式有：

a_j＝arg max{|<g_i,e_{a}>|²:a∈D} (18)

式中，D＝{z∈C:|z|＜1}表示单位圆盘，e_{a}(z)是算法中定义的一种字典，<>符号表示内积；

步骤5-3、根据已有的系数a_q和g_q，更新下一项g_q+1，并计算部分分量的和；

其中涉及到的公式有：

式中，B_n′(z)是算法中选用的一种有理正交系，也叫Takenaka-Malmquist system；

步骤5-4、将g_q+1返回到步骤5-2，求出下一个系数a_q+1，如此循环，直到求出n′个系数a_j，最后得到离散传递函数公式；

式中，m_j表示第j项前面的系数；

步骤5-5、最后经过转换

得到最终的离散传递函数，公式即为：

步骤6、传递函数转化：将离散传递函数经过双线性变换得出连续传递函数；

步骤7、使用评估拟合精度高低的代价函数和反映模型预测能力和可靠性的时域验证，对步骤6中连续传递函数进行评价。

2.根据权利要求1所述的基于自适应傅里叶分解的频域辨识方法，其特征在于，步骤2中的剔除野值运用的是拉依达准则，首先计算测量数据列的均值以及标准偏差，以置信概率99.7％为标准，经判断若序列值大于三倍标准偏差，则属于粗大误差然后剔除，若小于则予以保留；

具体公式为：

式中，x_i是一组数据序列，i＝1,2…n，n为数据序列个数，

是该序列的均值，σ是序列的标准偏差。

3.根据权利要求1所述的基于自适应傅里叶分解的频域辨识方法，其特征在于：步骤2中的去除趋势项具体是指先用最小二乘拟合出数据的趋势线，然后减去该趋势线；其中的数字滤波是指先根据实验对象的要求设计IIR滤波器或者FIR滤波器，然后将数据通过滤波器进行滤波。

4.根据权利要求1所述的基于自适应傅里叶分解的频域辨识方法，其特征在于，步骤3中的层叠分窗法是一种谱估计方法，其具体方法包括以下步骤：

步骤3-1、计算时域数据的傅里叶变换；

其中，x(n),y(n)是时域数据记录，X(f_k),Y(f_k)是第k个频率点的傅里叶系数，w(n)是离散的窗函数，Δt是采样时间间隔，

是频率分辨率，N是离散频率点的个数；

步骤3-2、计算谱函数的粗略估计值；

其中，

为第k个窗口计算得出的输入自谱，

为第k个窗口计算得出的输出自谱，

为第k个窗口计算得出的互谱，T_win为窗口长度；

步骤3-3、计算谱函数的平滑估计；

输入自谱的平滑估计：

输出自谱的平滑估计：

互谱的平滑估计：

其中，U是窗函数功率恢复系数，用于修正谱密度的幅值；m是窗口片段的数量；

根据得到的谱函数，引出相干函数估计值

它在各频率点上的定义为：

相干函数表示了在频率点f，由输入自谱线性地引起的输出自谱的部分占整体输出自谱的比例，其取值范围是0～1；若相干函数值满足大于等于0.6且不振荡，则认为辨识精度可接受。

5.根据权利要求1所述的基于自适应傅里叶分解的频域辨识方法，其特征在于，步骤7中的时域验证是通过比较在同样激励下模型的输出和测量的系统输出来实现对模型的评价；具体使用均方误差MSE这一指标，公式如下：

式中，x_i为系统真实的输出序列，

为模型预测的输出序列；i＝1,2,…,n；

代价函数是用来描述所得到的传递函数的拟合精度，公式如下：

其中，||符号表示每个频率ω处的幅值，∠符号表示每个频率ω处的相位，

表示频率响应估计值，T表示传递函数模型，n_ω表示频率采样点的数量，ω₁和

表示拟合的起始频率值和结束频率值，W_g和W_p分别是幅值和相位平方误差的相对权重；

W_γ是加权函数，公式如下：

式中，r_xy是每个频率ω处的相干函数值。

Images