一种基于阶跃响应相关分析的非参数系统辨识方法

摘要

本发明公开了一种基于阶跃响应相关分析的非参数系统辨识方法，用于利用电信号进行系统辨识的技术领域。本方法包括：对m序列激励系统后产生的系统输入信号与输出信号求互相关；确定系统的输入信号与输出信号的互相关除去直流分量的部分与系统的冲激响应的关系；建立系统的阶跃响应与系统输入信号与输出信号互相关除去直流分量的部分之间的时域对应关系；对工程仿真及应用中的具体离散系统或连续系统，在获得系统的阶跃响应后，根据系统的阶跃响应与冲激响应的关系求得系统的冲激响应。本发明消除了传统方法采用采用泰勒级数展开舍去冲激响应的高阶分量求近似所带来的系统响应的分量丢失的问题，极大地改善了相关分析法的辨识精度。

说明书

技术领域

本发明属于利用电信号进行系统辨识的技术领域，尤其涉及一种利用互相关来计算未 知系统的阶跃响应的方法，主要应用于航空航天、地球物理勘探、电力系统检测、环境系 统分析等研究领域。

背景技术

系统辨识(SystemIdentification)是根据系统的输入输出时间函数来确定描述系统行为的 数学模型，是现代控制理论中的一个分支。系统辨识所建立的辨识模型是对象输入输出特性 在某种准则意义下的一种近似，近似的程度取决于人们对系统先验知识的认识和对数据集性 质的了解程度，以及所选用的辨识方法是否合理。

目前非参数系统辨识方法主要采用基于系统输入和输出的互相关而进行系统响应分析的 相关分析法。现在常用的相关分析法系统辨识方法主要是李白南教授的处理方法，出自李白 南教授的《伪随机信号及相关辨识》P29-34，科学出版社出版，出版时间为1987-04-01。用 户利用自己掌握的系统的输入与输出的互相关算得系统的冲激响应如下：

$h\left(\tau \right)=\left(\begin{array}{c}[R_{xy}\left(\tau \right)-R_{yx}\left(\Delta t\right)]/S,\tau >\Delta t\\ [R_{xy}\left(\tau \right)-R_{yx}\left(\Delta t\right)]/\left[\left(\frac{1}{2}+\frac{\tau }{\Delta t}\left(1-\frac{\tau }{2\Delta t}\right)\right)S\right],0<\tau <\Delta t\end{array}\right)$

h(0)＝2[R_xy(0)-R_yx(Δt)]/S

其中中间参数$S=(1+\frac{1}{N})a^2\Delta t;$

式中：

τ为系统辨识过程中的时间变量；

a为m序列的幅值；m序列是指伪随机序列；

N为n阶m序列的单周期长度；

Δt为m序列的位宽；

R_xy(τ)为系统输入与输出信号的互相关；

R_yx(τ)为系统输出与输入信号的互相关；

h(τ)为系统的冲激响应；

常用的互相关分析法是针对输入与输出互相关进行近似运算的，这种近似是通过对系统 的冲激响应进行泰勒级数展开，在省去冲激响应的高阶导数分量的基础上建立的。目前相关 分析法广泛应用于地球物理勘探、电力系统检测等经济国防等领域。但是常用的相关分析法 存在缺陷，主要有如下三个方面的不足：

1)、对系统冲激响应的起始区间辨识失真度较大，尤其是初值及主要起始点的辨识；

2)、传统的相关分析法采用的冲激响应不变法，其频率响应中频谱混叠现象较严重；

3)、系统的幅频响应与相频响应与系统实际的频率响应相差较大，不利于系统辨识的频 域分析。

其中，前1)和2)两个不足虽然在一定程度上改善了直接用系统输入与输出的互相关近 似系统冲激响应所带来的失真，但是随着系统辨识技术的发展，辨识精度要求不断提高，创 新研发更高精度的系统辨识新方法就显得尤为重要。

即便前1)和2)瓶颈解决后，第三个不足的限制仍然使得目前的非参数辨识在频域的处 理上，尤其是采用滤波器抗噪上还会面临很多处理上的难题。

发明内容

本发明的目的是提供了一种利用系统的输入与输出数据的互相关来直接计算系统阶跃响 应的方法，能解决现有相关分析法直接估计待测系统的冲激响应所带来的误差较大与系统主 要信息丢失等问题。

本发明提供了一种基于阶跃响应相关分析的非参数系统辨识方法，实现步骤如下：

第一步，对m序列激励系统后产生的系统输入信号与输出信号求互相关；

$\begin{array}{c}{R}_{xy}\left(\tau \right)=\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)y\left(t+\tau \right)dt=\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)x\left(t+\tau -t_1\right){\mathrm{dt}}_{1}dt\\ =\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)x\left(t+\tau -t_1\right){\mathrm{dtdt}}_1=\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)R_{xx}\left(\tau -t_1\right){\mathrm{dt}}_1\\ =R_{xx}\left(\tau \right)\otimes h\left(\tau \right)\end{array}---\left(1\right)$

其中，x(t)为系统的输入信号；y(t)为系统的输出信号；τ为系统辨识过程中的时间变量； R_xy(τ)为系统的输入与输出信号的互相关；R_xx(τ)为系统输入信号的自相关；T₀为m序列的 相关长度；h(τ)为系统的冲激响应。

第二步，确定系统的输入信号与输出信号的互相关除去直流分量的部分与系统的 冲激响应h(τ)的关系，如下：

${\tilde{R}}_{xy}\left(\tau \right)={\tilde{R}}_{xx}\left(\tau \right)\otimes h\left(\tau \right)---\left(2\right)$

其中，为系统的输入信号的自相关除去直流分量的部分；

将输入信号的自相关函数分解成两个长度为Δt的方波w(τ)的卷积形式，

${\tilde{R}}_{xx}\left(\tau \right)=w\left(\tau \right)\otimes w\left(\tau \right)---\left(3\right)$

将式(3)代入式(2)得到：

${\tilde{R}}_{xy}\left(\tau \right)={\tilde{R}}_{xx}\left(\tau \right)\otimes h\left(\tau \right)=w\left(\tau \right)\otimes w\left(\tau \right)\otimes h\left(\tau \right)---\left(4\right)$

第三步，建立系统的阶跃响应与系统输入信号与输出信号互相关除去直流分量的部分之 间的频域对应关系：

对式(4)求傅里叶变换，得到如下：

${\tilde{R}}_{xy}\left(j\Omega \right)={\tilde{R}}_{xx}\left(j\Omega \right)H\left(j\Omega \right)={\left[W\left(j\Omega \right)\right]}^{2}H\left(j\Omega \right)---\left(5\right)$

其中，Ω是连续信号与系统中的模拟角频率，W(jΩ)为方波信号w(τ)经过傅里叶变换得 到的；H(jΩ)为频率响应；方波的信号傅里叶变换W(jΩ)表示为下式：

$W\left(j\Omega \right)=\frac{1-e^{-j\Omega \Delta t}}{j\Omega }---\left(6\right)$

将式(6)代入式(5)，得到下式：

${\tilde{R}}_{xy}\left(j\Omega \right)=W\left(j\Omega \right)\frac{1-e^{-j\Omega \Delta t}}{j\Omega }H\left(j\Omega \right)---\left(7\right)$

由于系统的频率响应H(jΩ)与系统阶跃响应g(t)的傅里叶变换G(jΩ)具备以下对应关 系：

$G\left(j\Omega \right)=\frac{1}{j\Omega }H\left(j\Omega \right)---\left(8\right)$

将式(8)代入式(7)得到下式：

${\tilde{R}}_{xy}\left(j\Omega \right)=W\left(j\Omega \right)G\left(j\Omega \right)(1-e^{-j\Omega \Delta t})---\left(9\right)$

所以：

$G\left(j\Omega \right)=\frac{{\tilde{R}}_{xy}\left(j\Omega \right)}{W\left(j\Omega \right)(1-e^{-j\Omega \Delta t})}---\left(10\right)$

根据上式能进一步获得系统的阶跃响应与系统输入信号与输出信号互相关之间在频域上 的对应关系；

其中，对上式中G(jΩ)求傅里叶反变换得到阶跃响应g(t)＝F^-1[G(jΩ)]；t表示时间变量；

对阶跃响应求微分得到系统的冲激响应h(t)：

$h\left(t\right)=\frac{d}{dt}g\left(t\right)---\left(11\right)$

对冲激响应h(t)求傅里叶变换得系统的频率响应H(jΩ)：

$H\left(j\Omega \right)=F\left[h\right(t\left)\right]={\int }_{-\infty }^{+\infty }h\left(t\right)e^{-j\Omega t}dt---\left(12\right)$

第四步，在实际实验仿真与工程应用中，处理的实际系统往往都为离散系统，因此需要 将估计出来的阶跃响应进行离散化处理，得到采样后的阶跃响应g(n)；

g(n)＝g_a(t)|_t＝nT＝g_a(nT)(13)

上式中T是抽样周期，n是采样个数，g_a(t)指的是连续系统的阶跃响应；

对阶跃响应g(n)求差分得到离散的冲激响应h(n)：

h(n)＝g(n)-g(n-1)(14)

再对冲激响应h(n)求离散傅里叶变换得到系统的频率响应H(jω)，如下：

$H\left(j\omega \right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}h\left(n\right)e^{-j\omega n}---\left(15\right)$

其中离散系统与连续系统的对应关系是：

H(jω)＝H(jΩT)(16)

ω是离散信号与系统中的数字角频率。

相对于现有技术，本发明的优点和积极效果在于：本发明方法在系统的输入与输出信号 的互相关计算系统阶跃响应上，创新地提出并建立了系统阶跃响应与互相关的直接对应，克 服了传统方法利用系统的互相关计算系统的冲激响应的不足，保留了最完整的系统响应信息， 较大地提高了系统的辨识精度。

附图说明

图1为输入信号的自相关的示意图；

图2为将m序列自相关函数处理成两个方波信号卷积的示意图；

图3为本发明非参数系统辨识方法的流程示意图；

图4为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统阶跃响应辨识 的比较示意图；

图5为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统冲激响应辨识 的比较示意图；

图6为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统冲激响应辨识 的冲激分量细节比较示意图；

图7为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统幅频响应辨识 的比较示意图；

图8为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统相频响应辨识 的全频段比较示意图；

图9为在连续系统中本发明方法与目前常用的非参数系统辨识方法对系统相频响应辨识 的低频部分比较示意图。

具体实施方式

下面结合附图和实施例对本发明的技术方案进行说明。在以下的描述中，将描述本发明 方法在不同系统上的应用，然而，对于本领域内的普通技术人员而言，可以仅仅利用本发明 的部分或者全部方法来实施本发明。为了解释的明确性而言，阐述了特定的阶数、序列重复 周期数和过采样倍数，但是很明显，在没有这些特定细节的情况下也可以实施本发明。在其 他情况下，为了不混淆本发明，对于一些众所周知的特征将不再进行详细阐述。

m序列激励系统后输出信号与输入信号的互相关与输入信号的自相关的之间存在着如图 1所示的对应关系。图中，R_xx(t)是m序列的自相关；Δt为m序列的位宽；a为m序列的幅 度。

如图2中，根据三角波可以表示成两个方波卷积的原理，将m序列的自相关函数R_xx(t)中 除去直流分量的部分分解成两个长度为Δt的方波w(t)的卷积形式，如下式：

${\tilde{R}}_{xx}\left(t\right)=w\left(t\right)\otimes w\left(t\right)---\left(1\right)$

本发明方法利用这种分解方式可进行公式化简。

下面结合图3对本发明的基于阶跃响应相关分析的非参数系统辨识方法的实现步骤进行 说明。

第一步，对m序列激励系统后产生的系统输入信号与输出信号求互相关。

根据互相关的性质与系统的卷积关系有如下化简：

$\begin{array}{c}{R}_{xy}\left(\tau \right)=\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)y\left(t+\tau \right)dt=\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)x\left(t+\tau -t_1\right){\mathrm{dt}}_{1}dt\\ =\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)\frac{1}{T_0}\underset{0}{\overset{T_0}{\int }}x\left(t\right)x\left(t+\tau -t_1\right){\mathrm{dtdt}}_1=\underset{0}{\overset{\infty }{\int }}h\left(t_1\right)R_{xx}\left(\tau -t_1\right){\mathrm{dt}}_1\\ =R_{xx}\left(\tau \right)\otimes h\left(\tau \right)\end{array}---\left(2\right)$

其中，R_xy(τ)为系统的输入与输出信号的互相关；R_xx(τ)为系统输入信号的自相关；T₀为m序列的相关长度；h(τ)为系统的冲激响应；x(t)为系统的输入信号；y(t)为系统的输出 信号。

第二步，将输入信号m序列的自相关函数分解成两个方波的卷积形式，建立输入与输出 的互相关、输入自相关、系统冲激响应的卷积对应关系。本步骤确定了系统的输入信号与输 出信号的互相关除去直流分量的部分与系统冲激响应h(τ)的关系。

由公式(2)可得，在利用m序列激励系统时，有如下对应关系：

${\tilde{R}}_{xy}\left(\tau \right)={\tilde{R}}_{xx}\left(\tau \right)\otimes h\left(\tau \right)=w\left(\tau \right)\otimes w\left(\tau \right)\otimes h\left(\tau \right)---\left(3\right)$

第三步，通过对上一步中卷积关系的在频率域中的分解，求反傅里叶变换，确定时间域 中系统的输入信号与输出信号的互相关与系统阶跃响应的关系。本步骤建立系统的阶跃响应 与系统输入信号与输出信号互相关除去直流分量的部分之间的频域对应关系以及相应的时域 对应关系。

下面对公式(3)进行傅里叶变换，然后进行公式推导以获取系统的阶跃响应与输入信号与 输出信号互相关除去直流分量的部分的时域对应关系。

对连续系统做傅里叶变换，对于离散系统做离散傅里叶变换。以连续系统为例，对式(2) 求傅里叶变换：

${\tilde{R}}_{xy}\left(j\Omega \right)={\tilde{R}}_{xx}\left(j\Omega \right)H\left(j\Omega \right)={\left[W\left(j\Omega \right)\right]}^{2}H\left(j\Omega \right)---\left(4\right)$

根据上式可以得到系统输入与输出互相关的傅里叶变换、方波信号的傅里叶变换和系统 频率响应之间的对应关系，

方波的信号傅里叶变换为：

$W\left(j\Omega \right)=\frac{1-e^{-j\Omega \Delta t}}{j\Omega }---\left(5\right)$

将式(5)代入式(4)可得：

${\tilde{R}}_{xy}\left(j\Omega \right)=W\left(j\Omega \right)\frac{1-e^{-j\Omega \Delta t}}{j\Omega }H\left(j\Omega \right)---\left(6\right)$

由于系统的频率响应H(jΩ)与系统阶跃响应g(t)的傅里叶变换G(jΩ)具备以下对应关 系：

$G\left(j\Omega \right)=\frac{1}{j\Omega }H\left(j\Omega \right)---\left(7\right)$

所以:

$G\left(j\Omega \right)=\frac{{\tilde{R}}_{xy}\left(j\Omega \right)}{W\left(j\Omega \right)(1-e^{-j\Omega \Delta t})}---\left(8\right)$

根据上式能进一步获得系统的阶跃响应与系统输入信号与输出信号互相关之间在频域上 的对应关系；

其中，对上式中G(jΩ)求傅里叶反变换可以得到阶跃响应g(t)＝F^-1[G(jΩ)]；

对阶跃响应求微分可以得到系统的冲激响应

$h\left(t\right)=\frac{d}{dt}g\left(t\right)---\left(9\right)$

对冲激响应h(t)求傅里叶变换得系统的响应：

$H\left(j\Omega \right)=F\left[h\right(t\left)\right]={\int }_{-\infty }^{+\infty }h\left(t\right)e^{-j\Omega t}dt---\left(10\right)$

第四步，在实际实验仿真与工程应用中，处理的实际系统往往都为离散系统，因此需要 将估计出来的阶跃响应进行离散化处理，得到采样后的阶跃响应，

g(n)＝g_a(t)|_t＝nT＝g_a(nT)(11)

上式中T是抽样周期，n是采样个数，g_a(t)指的是连续系统的阶跃响应；

对阶跃响应g(n)求差分得到离散的冲激响应h(n)：

h(n)＝g(n)-g(n-1)(12)

求得系统的冲激响应后，再对冲激响应h(n)求离散傅里叶变换可以得到系统的频率响应 H(jω)，ω是离散信号与系统中的数字角频率即：

$H\left(j\omega \right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}h\left(n\right)e^{-j\omega n}---\left(13\right)$

其中离散系统与连续系统的对应关系是：

H(jω)＝H(jΩT)(14)

实施例

为了解决现有相关分析法非参数系统辨识在系统响应初始区间内精度不高的问题，本发 明实施例提供了一种基于m序列激励未知系统后的输入与输出信号的互相关来高精度地计算 系统阶跃响应的方法，如图1所示，本实施例所述方法包括以下步骤：

步骤S101：设定m序列的初始参数，如阶数、采样率、过采样倍数、幅值、周期重复次 数等；

步骤S102：对设定好的m序列激励系统后产生的系统输入与输出信号求互相关，如公式 (2)所示；

步骤S103：建立输入信号的自相关的交流分量、系统冲激响应与系统输入与输出信号的 互相关的交流分量之间的卷积对应关系，如公式(3)所示。

步骤S104：对所述步骤S103中所述关系式做傅里叶变换，推导出系统的阶跃响应与互 相关除去直流分量部分的时域对应关系。

步骤S105：根据步骤S104获取系统的冲激响应。

通常系统非参数辨识需要得系统的冲激响应，现有的相关分析法由于广泛应用泰勒级数 的展开来建立系统输入、输出信号的互相关与系统冲激响应的对应关系，现有方法无法回避 利用泰勒级数近似计算所带来的系统响应的主要信息(尤其是冲激分量)丢失等问题。而本 实施例方法未涉及泰勒级数展开来近似计算系统的响应、也未涉及系统响应分量的舍去。相 应地，本实施例方法保留了待测系统响应全部分量，从而提高了利用互相关来辨识系统的精 度。

本实施例针对利用系统输入与输出信号的互相关的传统系统辨识方法中误差较大的问 题，采用如下两点进行改进：

1)直接建立系统输入与输出信号的互相关与系统的阶跃响应的对应关系，采用阶跃响应 不变法进而求得系统的冲激响应、幅频响应与相频响应。

2)将时域上的对应关系通过傅里叶变换转移到频域上，在频域上完成等式的化简与变形， 保留了系统响应最完整的信息量。

以下将现有清华大学李白南教授的相关分析法系统辨识和本发明的系统辨识法进行对比 描述。图4～图9中提到的“常用方法”出自《伪随机信号及相关辨识》李白南著，P29-34。 现有李白南教授的相关分析法和本发明的系统辨识方法计算的系统的阶跃响应如图4所示。 其中李白南教授的相关分析法系统辨识方法算得冲激响应再积分。从图4可以看出，本发明 的系统辨识方法辨识的精度明显高于现有相关分析法辨识的精度。

将现有辨识方法估计的系统频率响应的曲线、本发明的系统辨识方法计算的系统频率响 应的曲线、待测系统实际的频率响应曲线如图5所示；可以看出，本发明方法的辨识效果明 显优于现有的相关辨识方法，辨识精度有很大的提高。通过图6～图9的效果对比，也可以 明显看出采用本发明方法比采用现有的相关辨识方法，辨识精度有很大的提高。因此可以看 出，本发明通过建立输入、输出的互相关及其高阶导数时移分量的累加与系统阶跃响应的直 接对应关系，消除了传统方法采用采用泰勒级数展开舍去冲激响应的高阶分量求近似所带来 的系统响应的分量丢失的问题，极大地改善了相关分析法的辨识精度。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何 熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应 涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围 为准。