基于宽频谱激励信号的系统模型参数快速辨识方法与流程

本发明涉及一种基于宽频谱激励信号的系统模型参数快速辨识方法，属于惯性测试技术领域。

背景技术：

系统模型是对系统本质属性的一种描述，建立系统数学模型的方法分为理论分析法和测试法，若想设计一个性能满足要求或改善现在系统性能的方式就是清楚的了解其动态特性。然而，在实际控制工程中，有纷繁复杂的各式被控对象，由于其构造方式复杂、运行机理模糊等限制往往使通过理论分析法建模获取其动态特性变得很困难。因此，在实际控制工程中往往采用理论建模法与测试法相结合的方法。首先，通过先验理论知识对系统特征进行判断并获得机理已知参数，如阶数、阈值、带宽等；接着，根据上述约束条件设计合理的试验辨识机理未知的参数；最后，将模型参数进行综合并对已获得参数可靠性进行检验。其中，测试法通过获取系统的输入、输出信号，并对其进行特殊处理可获得系统的动态信息，接着根据相应的辨识方法来确定模型参数。目前典型的且常用的输入信号为脉冲信号、阶跃信号、斜坡信号、正弦信号、方波信号等。在试验过程中，通过将上述信号作为激励信号输入系统中，在输出端进行实时记录并绘图即可得到对应的系统响应。以阶跃响应为例，阶跃响应是指系统在阶跃输入信号的作用下所产生的零状态输出响应，由于其信号产生方式容易，因此在工程上得到广泛的应用。阶跃响应能够反映系统的动态特性，可用于辨识系统模型的参数。在辨识系统模型参数方面，目前采用的方法包括：一、选取阶跃响应曲线上某些时间点tk上的观测数据y(tk)，利用观测数据构造关于参数的方程组，研究代数求解的两点法/三点法等，求出待辨识的传递函数的参数(比如一阶系统和二阶系统的时间常数t和增益k等)；二、设定模型阶数、辨识准则后对输入输出数据进行数据拟合，选择合适的估计方法估计模型中的未知参数。

丁峰教授在《传递函数辨识(1):阶跃响应两点法和三点法》、《传递函数辨识(2):脉冲响应两点法和三点法》等连载论文中论述了通过给系统施加某种特定的输入信号，记录系统的响应，获取系统的动态信息。在输入信号的选择上，采用的均是典型的输入信号、如脉冲信号、阶跃信号、斜坡信号以及周期信号等，进一步利用这些动态信息，提出相应的辨识方法来确定系统的传递函数参数模型。

金宇晖在专利《一种结合阶跃与正弦激励的二阶系统辨识方法》中提出通过对二阶系统进行经典信号结合的分步激励的方法，通过制图、取点分析来获得传递函数的特征表达式。但在每次激励信号施加的过程中，其频域特性并未得到改变，信号能容的丰富程度并不能通过分布激励而实现频谱融合，因此信号内容仍较为单一。

技术实现要素：

本发明的目的是为了解决上述现有技术存在的问题，进而提供一种基于宽频谱激励信号的系统模型参数快速辨识方法。

本发明的目的是通过以下技术方案实现的：

一种基于宽频谱激励信号的系统模型参数快速辨识方法，包括以下步骤：

步骤1：激励信号生成

1)设计激励信号y，使激励信号产生的过载在允许范围之内，其表达式为：

式中，fc表示截止频率；f0表示初始频率；t表示激励信号时长；n0表示衰减系数；g表示幅值增益，t表示时间变量，n为时间变量的阶数。n0、n、f0、fc均可以根据实际需要进行参数值设置；

2)激励信号所表示的是位置信息，通过对激励信号求一阶导数和二阶导数，得到其速度和加速度；

速度表达式如下：

加速度表达式如下：

其中，y′表示激励信号的实时速度信息，y″表示实时信号的加速度信息；

3)该时域连续信号在频域的频谱表示为覆盖f0至fc的矩形区域，设置频带阈值；引入如下表达式：

公式(4.1)表示当fc>f0时，求任意时刻的频率值或求产生任意频率所在的时刻，其中ft表示t时刻频率，其中f0、fc为根据试验需求设定的初始值，t为试验时长，因此，对于ft和t，已知任意时刻t即可对应的求相应时刻的信号频率ft，反之亦然；

公式(4.2)表示当fc<f0时，会出现频率归零的这一现象，tz为到达频率为0的时间。

公式(4.1)和公式(4.2)使得对于输入信号，即激励信号y，可以获取任意频率的时间和获取任意时间的频率；

在惯性仪表的测试过程中，通过对其工作方式及工作性能的了解，可以得出其工作频率范围，针对某一频率附近的性能进行挖掘时，设置初始频率f0在有效工作范围内，

当fc>f0时，在进行持续时长为t的试验中，到达设置任意频率的时间即为t可求，在试验结束的t时刻惯性仪表的工作频率即为截止频率fc；

当fc<f0时，在进行持续时长为t的试验中，由于衰减因子的存在会使二阶倒数y″达到峰值，即频率为0，到达频率为0的时间即为tz，这种情况下取在0-tz时刻试验数据分析；

步骤2：相关参数计算方法

1)惯性测试设备传输延迟

采用所述激励信号y，并设定初始频率f0和截止频率fc，惯性测试设备设定为追踪模式，传输延迟可以通过对输入数据和输出数据进行互相关计算求得，对于离散数据，互相关计算公式如下：

rry互相关参数，参数值越大，参数值越趋近于0则越不相关；r^*(t)表示r(t)的共轭函数，其中r(t)代表输入数据，y(t+τ)代表输出数据，k代表采样数量，τ为时间延迟；

2)惯性测试设备的比例因子k和测量偏差ε的确定

通过对公式(5)求解可得惯性测试设备传输延迟，对原有输出数据进行消除延迟处理后，对惯性测试设备测量误差进行自相关计算可得惯性测试设备测量偏差ε，

其中，代表自相关参数，参数值越趋近于0则越不相关；

对公式(7)计算结果取峰值，对峰值进行开方计算即可得惯性测试设备测量偏差ε；

惯性测试设备的比例因子k可由下式求得：

3)惯性测试设备噪声

惯性测试设备噪声可以用自相关的均方值和均值平方的差值来确定，表达式如下：

其中与表示的最大值和最小值；

4)惯性测试设备动态特性的两种计算方法

①首先，对输入数据和输出数据分别进行快速傅里叶变换，可得对应的频域离散表达式；接着，计算对应点输出数据与输入数据的比值，该结果可等效为惯性测试设备的脉冲频率响应；然后，对其进行反傅里叶变换即可得脉冲响应；最后，采用传统脉冲响应辨识方法即可得系统传递函数；

②由于频域内传递函数可表示为实部与虚部加和或幅值和相角乘积，因此应用输入数据和输出数据绘制bode图，由于惯性测试设备是二阶系统为主导的系统，因此通过对相频特性图滞后90°处进行数据读取即可获得惯性测试设备的固有频率，由固有频率可求系统的阻尼比；

如上所述，通过施加一次激励信号，即可通过计算得到一系列惯性测试设备有关参数。

所述步骤1中，出于对惯性仪表测试设备性能的考虑，令n＝2或令n＝3。

所述步骤1中，根据实际惯性仪表工作特征可以自由设置频带阈值。

本发明的有益效果为：

本发明所设计的激励信号形式新颖、内容丰富、可控性强、灵活度高在测试领域创新性强；相较于现有技术在信号内容丰富程度上有了极大的提升，从信号在频域的频谱特性可以看出，将传统信号的单频峰值或窄频波动幅值激励改善为宽频且幅值可控的激励信号，大大改善且提高了信号的频域特征。

本发明所设计的激励信号相较于现有技术在对系统参数的激励机理上有了更全面的改善，考虑到系统在实际工程中所面临的工作状态、工作环境绝非单纯的阶跃信号激励等传统信号。因此，为了更真实的获得系统的动态特性，更逼真的还原系统的实际工作状态，在对系统的激励上选择了可控频域宽度的激励信号，更全面的激励系统内的各项参数，使测试结果可靠性提升。

本发明所设计的激励信号相较于现有技术在测试试验设计冗杂度，测试试验时间成本上有了极大的改善。现有技术往往需要传统信号相结合的多次试验来确定系统的参数，对于每一步骤的零初始状态皆会存在频差，从而影响试验结果。而本发明通过一次短时激励信号即可达到目的，实现了高效测试、快速测试的目的。

本发明所设计的激励信号相较于现有技术在系统模型参数获得上具备更全面的优势，在一次信号激励中不仅可获得动态模型参数，同时可以获得系统噪声参数、测试设备偏差、测试设备比例增益以及测试系统传输延迟等参数。

本发明所设计的激励信号相较于现有技术在系统模型参数辨识上通过对所有采样点进行数据处理试验最优化拟合，而不是通过取特殊点进行近似拟合或通过响应曲线作图取点。在拟合精度以及模型参数的获取上可靠性更高。

本发明采用特征全面的信号大大减少了测试设计复杂度、测试时间成本以及规避了多次调零的初始误差；并且对输入、输出数据的处理上具有快速性、全面性，通过短时激励即可获得系统的重要参数，改善了现有技术。

附图说明

图1为脉冲信号频域特性图。

图2为斜坡信号频域特性图。

图3为正弦信号频域特性图。

图4为扫频信号(无衰减因子)频域特性图。

图5为本发明改进后扫频信号频域特性图。

图6为二阶位置信号波形曲线图。

图7为三阶位置信号波形曲线图。

图8为二阶速度信号波形曲线图。

图9为三阶速度信号波形曲线图。

图10为二阶加速度信号波形曲线图。

图11为三阶加速度信号波形曲线图。

具体实施方式

下面将结合附图对本发明做进一步的详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式，但本发明的保护范围不限于下述实施例。

如图1至图11所示，本实施例所涉及的一种基于宽频谱激励信号的系统模型参数快速辨识方法，包括以下步骤：

步骤1：激励信号生成

目前现有测试方法的输入信号大多为定值、匀速信号、匀速斜坡信号或单频正弦信号，频谱内容单一，不能满足惯性仪表测试设备工作频带内的全覆盖动态激励。基于此，本实施例在激励信号的设计中首先考虑信号频谱内容的丰富程度，其次，通过对惯性仪表测试设备工作环境的了解，使激励信号产生的过载在允许范围之内，其表达式为：

式中，fc表示截止频率；f0表示初始频率；t表示激励信号时长；n0表示衰减系数；g表示幅值增益，t表示时间变量，n为时间变量的阶数，a和b的引入是为了简化公式表示形式。n0、n、f0、fc均可以根据实际需要(设备承受范围以及元器件工作环境)进行参数值设置。图1至图5所示为各输入信号的频域特性。

考虑到惯性仪表测试设备在设计时的性能上限，如可承受最大过载、载荷等因素的限制，对于激励波形，引入n0和g的分别用来限制加速度和速度的幅值裕度。这些调节参数通过调节波形使惯性仪表和惯性仪表测试设备在输入速度和加速度值较大时不会损坏，并可以有效限制加速度、速度位置指令限制在硬件的测试能力范围内。激励信号所表示的是位置信息，对于其速度和加速度可分别对其求一阶导数和二阶导数。

速度信息表达式如公式(2)：

加速度信息表达式如公式(3)：

y′表示激励信号的实时速度信息，y″表示实时信号的加速度信息，出于对惯性仪表测试设备性能的考虑，令n＝2或令n＝3，对于激励信号的位置、速度和加速度的波形曲线如图6至图11所示。

通过波形变换速率可分析随着表示时间阶数量值的增加，动态过程开始更缓慢，但动态过程的加剧速度很快，所以针对需要确定时间阶数即可。

该时域连续信号(时域连续信号的定义是在一段时间上连续的信号，因为输入y激励信号是连续的信号，且激励时间可控，因此其满足时域连续信号的定义，即激励时间内为连续信号)，在频域的频谱表示为覆盖f0至fc的矩形区域；根据实际惯性仪表工作特征可以自由设置频带阈值。为了更全面的完善试验，引入如下表达式：

公式(4.1)表示当fc>f0时，求任意时刻的频率值或求产生任意频率所在的时刻，其中ft表示t时刻频率，其中f0、fc为根据试验需求设定的初始值，t为试验时长，因此，对于ft和t，已知任意时刻t即可对应的求相应时刻的信号频率ft，反之亦然；

公式(4.2)表示当fc<f0时，会出现频率归零的这一现象，tz为到达频率为0的时间。

公式(4.1)和公式(4.2)使得对于输入信号，即激励信号y，可以获取任意频率的时间和获取任意时间的频率；对于实际的试验操作及数据分析提供方便。下面针对一次完整的试验设计来阐述信号中各个参数的实际意义。在惯性仪表的测试过程中，通过对其工作方式及工作性能的了解，可以得出其工作频率范围，若想针对某一频率附近的性能进行挖掘，则可根据上式进行求解，初始频率可在有效工作范围内任意设置，其中分为fc<f0和fc>f0两种情况。当fc>f0时在进行持续时长为t的试验中，到达设置任意频率的时间即为t可求，在试验结束的t时刻惯性仪表的工作频率即为截止频率fc。当fc<f0时在进行持续时长为t的试验中，由于衰减因子即衰减系数n0的存在会使二阶倒数达到峰值，即频率为0，到达频率为0的时间即为tz，这种情况下取在0-tz时刻试验数据分析即可。

步骤2：相关参数计算方法

(1)惯性测试设备传输延迟

采用前文所述激励信号，并设定初始频率和截止频率。惯性测试设备设定为追踪模式，传输延迟可以通过对输入数据和输出数据进行互相关计算求得，此处的输入数据是将激励信号进行仿真，并设置采样频率对其进行采样，如采样频率为500hz，时长为20秒，得到的离散数据即为输入数据，即为10000*1的数组；输出数据是指对在上述输入数据的激励下，对惯性测试设备的输出进行相同时长相同采样频率的采样，则同样为500hz，时长为20秒，也是一个10000*1的数组，为后文计算用。

对于离散数据，互相关计算公式如下：

rry互相关参数，参数值越大，参数值越趋近于0则越不相关。中为若r(t)和y(t+τ)均为连续函数的情况下互相关的计算公式，其中r^*(t)表示r(t)的共轭函数。考虑到r(t)和y(t+τ)在实际计算中均为离散的，因此采用的表示方法，其中r(t)代表输入数据，y(t+τ)代表输出数据，k代表采样数量。对互相关进行作图后可以通过峰值确定时间延迟τ。上述公式中，第一个等号代表连续函数互相关表达式，第二个等号代表离散形式下的互相关表达式，对于本实施例，可用第二个等式近似取值近似。

(2)惯性测试设备的比例因子k和测量偏差ε的确定

通过对公式(5)求解可得惯性测试设备传输延迟，对原有输出数据进行消除延迟处理后，对惯性测试设备测量误差进行自相关计算可得惯性测试设备测量偏差ε。

公式(7)中代表自相关参数，参数值越趋近于0则越不相关。表示当和和均为连续函数的情况下自相关的计算公式，其中代表的共轭函数；考虑到和在实际计算中均为离散的，因此采用来表示离散情况下的计算方法，在步骤2的(1)中提及到惯性测试设备的传输延迟的计算方法中可以的到时间延迟τ的参数值，得到参数值后对输出数据进行延迟调整后，并对输入数据和输出数据做差值运算，结果即为其物理含义为t时刻的偏差，为间隔为σ时的函数，自相关函数为偶函数，且一般当σ＝0时，自相关参数最大。公式(7)中，第一个等号代表连续函数自相关表达式，第二个等号代表离散形式下的自相关表达式，对于本实施例，可用第二个等式取值近似。

对公式(7)计算结果取峰值，对峰值进行开方计算即可得惯性测试设备测量偏差ε。

惯性测试设备的比例因子k可由下式求得：

(3)惯性测试设备噪声

惯性测试设备噪声可以用自相关的均方值和均值平方的差值来确定。表达式如下：

其中与表示的最大值和最小值；

(4)惯性测试设备动态特性的两种计算方法

①首先，对输入数据和输出数据分别进行快速傅里叶变换，可得对应的频域离散表达式；接着，计算对应点输出数据与输入数据的比值，该结果可等效为惯性测设设备的脉冲频率响应；然后，对其进行反傅里叶变换即可得脉冲响应；最后，采用传统脉冲响应辨识方法即可得系统传递函数。

②由于频域内传递函数可表示为实部与虚部加和或幅值和相角乘积，因此应用输入数据和输出数据绘制bode图，由于惯性测试设备是二阶系统为主导的系统，因此通过对相频特性图滞后90°处进行数据读取即可获得惯性测试设备的固有频率，由固有频率可求系统的阻尼比。

如上所述，通过施加一次激励信号，即可通过计算得到一系列惯性测试设备有关参数，避免了繁杂的试验操作。在更逼真的模拟系统所处工作环境的同时大大缩短了测试时间成本。

以上所述，仅为本发明较佳的具体实施方式，这些具体实施方式都是基于本发明整体构思下的不同实现方式，而且本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。