改进实现主动吸收的IIR数字滤波器优化算法

Improved IIR Digital Filter Optimization Algorithm for Active Absorption

摘要: 在海洋工程试验水池中，进行模型试验时，造波机产生的波浪从池壁和结构模型上反射后生成反射波，传播到造波机处再次反射产生二次反射波，二次反射波干扰了水池的波浪场，从而降低了试验的精度，这对试验是不利的。为消除二次反射波常采用主动吸收系统，其中，力反馈式主动吸收系统考虑造波板上作用力，以力与造波板冲程之比为吸收传递函数，实时控制系统的输入与输出，使得生成的波逐渐接近理想波浪。为了使得造波板冲程和力保持最优平衡关系，本文基于理论传递函数，设计了IIR数字滤波器进行滤波，实现主动吸收控制。对于设计IIR数字滤波器，粒子群优化算法是一种非常有效的设计方法，该算法设计简单，运算速度快，但是该算法在优化过程中可能会局限于局部最优，而无法得到全局最优解。本文在更新粒子速度时，采用了一种新型的方法——即基于惯性权重改进学习因子，惯性权重与学习因子呈非线性函数关系，对比惯性权重线性递减、非线性递减的改进方法。仿真实验结果表明，该方法收敛速度较快，迭代次数较少，求解效率较高，优化结果最优。

Abstract: In the ocean engineering test basin, during the model test, the wave generated by the wave-makers generates reflected waves after being reflected from the basin wall and structural model, which propagates to the wave-making plates and reflects again to generate the secondary reflected wave. The secondary reflected wave interferes with the wave field of the basin, thus reducing the accuracy of the test, which is unfavorable to the test. An active absorption system is often used to eliminate the secondary reflection wave. The force-feedback type active absorption system takes the ratio of the force to the stroke of the wave-making paddles as the absorption transfer function, and controls the input and output of the system in real time, so that the generated wave gradually approaches the ideal wave. In order to maintain the optimal balance between the stroke and force of the wave-making paddles, IIR digital filter based on the theoretical transfer function is designed to realize the active absorption control. To design IIR digital filter, particle swarm optimization algorithm is a very efficient design method, the algorithm is simple and has fast calculation speed, but the algorithm in the optimization process could be limited to the local optimum, and can’t get the global optimal solution. In this paper, a new kind of method is used to update the particle velocity based on inertia weight to improve the learning factor. Inertia weight and learning factor have a nonlinear function relation. Compared with the improved method of linear decrease and nonlinear decrease of inertia weight, the simulation results show that this method has fast convergence speed, fewer iteration times, higher solving efficiency and optimal optimization results.

1. 引言

在深水海洋试验水池中，造波机是再现海洋波浪环境的主要装备，尤其是具备主动消波系统的造波机更是对试验的意义重大。但是，在造波过程中会产生一些影响试验的因素，比如二次反射波，这是因为在真实的海洋环境中，波浪遇到结构物产生的反射波会传播到远处，海洋无限大可看做无边界，那么反射波会有足够的空间传播并逐渐衰减至消失。但在海洋试验水池中，由于场地限制，水池边界作用明显，波浪传播至池壁处会产生反射波，部分反射波继而传至造波板处产生二次反射波，若不采取措施吸收二次反射波，那么水池内会形成多次反射，这种反射严重影响了实验中波浪特性，降低实验可靠度。因此，为了避免二次反射波影响试验所需波浪场，则必须进行反射波的吸收。

吸收技术分为两种——被动吸收和主动吸收。被动吸收是指在水池端部或侧向设置具有吸收波浪作用的消能器，通过拥有一定的透空率的装置消除反射波浪 [1]；而主动吸收是通过造波机自身来吸收反射波的方法，根据反馈信号的不同，主动吸收分为位移反馈式主动吸收和力反馈式主动吸收。力反馈式主动吸收一般应用于干背式造波机，其研究造波板受力仅考虑板前部分池水的作用，受力状态较为简单。因此接下来的研究基于力反馈式的干背式造波机主动吸收技术展开 [2]。

设计干背式吸收式造波系统，需要考虑几个关键的方面：输入信号的选择；控制器的结构；吸收滤波器的实现和优化。在造波过程中，造波板前的波浪场可以由几个物理量确定，其中，最重要的两个就是浪高和力，前者是通过浪高仪测量波浪高度，后者是通过力传感器测量施加的力，对于吸收式干背式造波机的实现主要就是依靠研究施加的力。

采集到力信号和造波板的速度信号后，有两种处理方式：一是，将力信号传入力信号吸收滤波器进行过滤，再与采集到的速度信号结合考虑，系统根据速度要求造波；二是，将被测速度信号进行过滤，与力信号结合，系统根据力要求造波。但是，最优吸收滤波器是反因果的，因此上述两种方式都是不切实际的，需要考虑因果近似，其细节取决于吸收滤波器的实现，关于吸收控制器的实现，目前，最常用的就是基于理想传递函数的无限脉冲响应滤波器，在IIR滤波器中，可以向较高频率分量的信号施加相位偏移或时间延迟，可以很好地达到吸收效果 [3]。

基于IIR无限脉冲响应滤波器的消波是基于时域的脉冲响应函数的近似，由于脉冲响应函数有因果关系，所以该法是可行的，然而，任何近似都会产生误差，直接计算脉冲响应进而研究频率吸收特性又是极其困难的，因此，须利用积分、迭代或经验数值等方式来设计最佳控制器。

IIR数字滤波器的设计就是使得所求的幅频特性尽可能地逼近理想幅频特性，设计方法主要分为两类：一是通过变换方法将模拟滤波器转换为数字滤波器，但是该方法设计效果不理想；二是直接利用优化算法设计数字滤波器，在优化准则下求解最优的滤波器系数，使设计的滤波器性能最优。目前，利用优化算法优化设计IIR数字滤波器是国内外研究的重点。

大多数字滤波器优化采用非线性最小二乘法，当独立系数的个数较少时，使用Matlab中的Isqnonlin函数可以有良好的收敛性(一阶低通滤波器一般有两个独立系数)，但是，当独立系数的个数较多时，就需要一个复杂的程序，这个程序可以对任何给定的输出参数进行优化，例如优化最优相位和幅值条件：滤波器的实部与流体动态响应的实部(水动力阻尼)相等，滤波器的虚部与流体动态响应的虚部(附加质量)互为相反数，使得滤波器和水动力系统的实部和虚部分别匹配；或者优化吸收功率系数，使得波浪吸收效果达到最优。优化过程必须限制特定的频率范围或带宽，经实验证明，上述两种方法结果差别不大，都可以很好地优化滤波器参数，接下来的研究主要在优化吸收功率的基础上利用改进的粒子群优化算法设计数字滤波器 [4]。

理想的吸收控制器使造波板冲程和力之间保持最优关系，但是，需要知道的是，一个理想的控制器是无法在实际中实现的，无论哪种技术都不可能达到理想的滤波器的幅频特性，只能是无限逼近。而且目前大部分算法存在早熟收敛问题，即在优化过程中会陷入局部最优，而无法得到全局最优解，收敛效果较差。粒子群优化算法是一种新的全局优化算法，算法模拟类似于鸟类捕食的行为，它在解决非线性优化问题方面有突出表现，同时算法简单，目前，它已成为一种重要的优化工具，基于粒子群优化算法，改进部分算法进程，将其应用在IIR数字滤波器的优化设计中，能够提高IIR数字滤波器系数的求解速度，减少误差，增强优化效果，为后期设计优良的吸收控制器提供算法依据，进一步改善主动吸收效果。

2. 吸收传递函数与数字滤波器

传递函数是由造波板的需要或要求输出与水池中实际产生的波浪之间的关系定义的。在需求信号被发送到造波板之前，传递函数被应用到需求信号，并对电学和机械的运行以及造波板的惯性进行校正。在理想情况下，对于任何给定的波浪分量，造波板的输出应该总是等于需求信号。

数字滤波器是一种离散时间系统，它可以根据期望响应特性对输入的离散数字信号进行滤波处理，它的作用主要分为两种，一是信号处于带通区间，此时频率响应幅值为1，信号幅值无失真输出；二是信号处于带阻区间，此时频率响应幅值为0，信号直接滤除。数字滤波器根据输入信号与截止频率的关系主要分为以下五种形式：低通滤波器、高通滤波器、带通滤波器、带阻滤波器、全通滤波器，数字滤波器系统的特性主要体现在幅频特性和相频特性，当数字滤波器系统的幅频特性为一恒定的常数(1或0)时，系统具有无限宽的通频带，且频率响应均匀，而相频特性曲线是经过原点的一条直线，这是为保证时域中各频率分量在时间轴上相对位置不变，只是整体滞后一个时间t₀。

根据脉冲响应的不同，数字滤波器又分为递归滤波器和非递归滤波器，递归滤波器即线性的、时不变系统，也称为无限脉冲响应滤波器(infinite impluse response filters，IIR滤波器)，非递归滤波器也称为有限脉冲响应滤波器(finite impluse response filters，FIR滤波器)。对于IIR滤波器，采用更少的滤波器系数模拟频率响应可以得到更快的数据进程，更适合应用于主动吸收系统中 [5]。

吸收传递函数代表了系统输入值与输出值之间的频域关系。在主动吸收系统中，吸收传递函数的幅频特性是一条变化的曲线，所以，为了达到主动吸收的效果，需要通过计算滤波器参数，使其在频率内与吸收传递函数具有相同的频率特性，使滤波器的频率特性曲线趋近于吸收传递函数的特性曲线，同时滤波器的相频特性曲线逼近吸收传递函数的相频特性曲线。

对于位移反馈式主动吸收系统，吸收传递函数表示为前进波的波幅与造波板冲程之比，即

$c_0=\frac{A}{X}$. (1)

对于力反馈式主动吸收系统，吸收传递函数表示为力与造波板冲程之比，即

$c_f=\frac{F}{X}$. (2)

二者之间的关系如下

$c_f=c_0\frac{1}{i\omega }\frac{1}{Z_d+Z_f}$. (3)

定义Z_d为系统动态反馈，主动系统的动态响应分析的阻抗分为两部分，一是驱动元件动态输出特性的分析；二是结构本身动态特性的分析，定义Z_f为吸收滤波器的特性。

在力反馈式主动吸收系统中，吸收传递函数代表造波板冲程与水动力之间的关系，而滤波器就起到了吸收传递的作用，对离散信号进行滤波处理，因此需要采用一种有效的优化设计方法，使得滤波器的频率响应和吸收传递函数的频率响应近似相等。

3. 粒子群优化算法

IIR数字滤波器在信号处理领域有非常重要的作用，目前滤波器设计方法主要分为两种：一是先设计模拟滤波器，然后利用脉冲响应不变法、阶跃响应不变法、双线性变换法等将其转换为数字滤波器，这种方法比较简单，但是设计误差大，结果不理想；二是直接使用优化算法，优化算法可以使得设计的滤波器性能最优，比如最小二乘法、粒子群算法、量子粒子群算法、人工蜂群算法、遗传算法、量子遗传算法等，其中，粒子群优化算法PSO (Partical Swarm Optimization)研究了鸟类等粒子所表现出的群体行为，如图1，将待优化的目标变量比作鸟类等粒子的位置信息，将目标函数的解比作鸟群要找的食物，最终通过鸟群内部一系列操作得到问题的最优解，该算法能够将更多的粒子聚集在最优粒子附近，从而增强优化能力。粒子群优化算法在解决非线性问题方面具有很好的表现，而且算法设计简单可行，成为解决许多非线性优化问题的有力工具。

与FIR滤波器相比，IIR滤波器的难点主要就是非线性和稳定性问题，但其可以通过更少的计算量来得到更优的结果。粒子群优化算法是一种更高级的优化方法，该算法简单，运行速度快，能够更好地解决非线性优化问题。但是该算法在运行时容易局限于局部最优而无法得到全局最优解，为了避免这个问题，做出如下改进：随着迭代的进行，让惯性权重递减，并基于动态调整权重下改进学习因子，从而更新粒子速度。通过改进，在算法初期，惯性权重较大，算法搜索范围也就大；到了算法晚期，惯性权重变小，算法的搜索能力变的更强更精细，更易找到全局最优解 [6]。

Figure 1. Principle of particle swarm optimization

图1. 粒子群优化算法原理图

粒子优化算法中粒子的每一个位置都代表问题的一个解，在改变粒子的位置的过程中，问题的一个个解也被搜索出来，粒子的适应度不同，解的优劣就不同。在每一次迭代中，算法中的随机粒子就根据下面两个公式来更新自己的速度和位置：

$V=w\ast V+c1\ast \text{rand}\left(\text{pBest}-\text{Present}\right)+c2\ast \text{rand}\left(\text{gBest}-\text{Present}\right)$,(4)

$\text{Present}=\text{Present}+V$.(5)

其中，pBest代表个体极值，即单个粒子在移动时所找到的最优解，gBest是整个粒子群在移动时所找的的最优解，代表全局极值，V是粒子的速度，Present是粒子的当前位置，rand( )是取值范围在(0, 1)之间的随机数，c1和c2是作学习因子，通常都取值为2。w是惯性权重，取值范围一般在(0.1, 0.9)内。

IIR滤波器的传递函数可以表示为 [7]：

$H\left(z\right)=A_0{\displaystyle \underset{k=1}{\overset{N}{\prod }}\frac{1+a_k{z}^{-1}+b_k{z}^{-2}}{1+c_k{z}^{-1}+d_k{z}^{-2}}}$. (6)

滤波器的频率响应为：

$H\left({\text{e}}^{j\omega }\right)=A_0{\displaystyle \underset{k=1}{\overset{N}{\prod }}\frac{1+a_k{z}^{-j\omega }+b_k{z}^{-2j\omega }}{1+c_k{z}^{-j\omega }+d_k{z}^{-2j\omega }}}=A_{0}G\left({\text{e}}^{j\omega }\right)$.(7)

采用频域最小均方误差优化设计IIR数字滤波器，假设滤波器的理想幅频特性为： $\left|H_d\left({\text{e}}^{j\omega }\right)\right|$，则在离散频率点上，使所设计的滤波器的幅频响应 $\left|H\left({\text{e}}^{j\omega }\right)\right|$ 与理想的幅频响应的均方误差：

$E={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\left[\left|H\left({\text{e}}^{j\omega }\right)\right|-\left|H_d\left({\text{e}}^{j\omega }\right)\right|\right]}^2}$ 为最小。从该式可以看出，目标函数E是一种非线性函数，它与增益A₀和

滤波器系数a_k，b_k，c_k，d_k ( $k=\text{1},\text{2},\text{3},\cdots$ )有关，其中N为二阶节个数。因此函数E有(4N + 1)个未知参数。

E可以表示为A₀和a_k，b_k，c_k，d_k的函数，即：

$E={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\left[\left|A_0\right|\cdot \left|G\left({\text{e}}^{j\omega }\right)\right|-\left|H_d\left({\text{e}}^{j\omega }\right)\right|\right]}^2}$. (8)

对于上式，先求出最佳增益A₀，再通过PSO算法得到滤波器系数a_k，b_k，c_k，d_k ( $k=\text{1},\text{2},\text{3},\cdots$ )。

增益A₀的解析表达式：

$A_0={\displaystyle \underset{i=1}{\overset{M}{\sum }}\left|G\left({\text{e}}^{j\omega }\right)\right|}\cdot \left|H_d\left({\text{e}}^{j\omega }\right)\right|/{\displaystyle \underset{i=1}{\overset{M}{\sum }}{\left|G\left({\text{e}}^{j\omega }\right)\right|}^2}$. (9)

在PSO算法中，适应度的大小决定了粒子位置的优劣，所以适应函数F必须符合实际问题要求，根据以上推导得到IIR滤波器的适应函数为：

$F=E={\displaystyle \underset{i=1}{\overset{M}{\sum }}{\left[\left|H\left({\text{e}}^{j\omega }\right)\right|-\left|H_d\left({\text{e}}^{j\omega }\right)\right|\right]}^2}$. (10)

从上式可以看出，当适应函数的值较小时，设计的滤波器的均方误差就小，滤波器系数就更优，优化设计更优良。算法运行过程就是寻找最小适应度的过程，当找到适应函数最小值时，最优解——即最优的滤波器参数也就得到了。

4. 改进惯性权重和学习因子

惯性权重代表着粒子保持前一刻运动状态的能力，学习因子代表粒子学习变化能力，如果粒子缺少变化，则容易导致局部最优，若粒子变化太大，则容易产生误差，引入惯性权重和学习因子这两种参数可以对PSO算法的运行进行全局调节，使得PSO算法的性能变得更加优良，不同的参数变化方式对算法运行的影响也不同 [8]。以下是三种改进方法。

1) 线性递减权重法

$w=w_{\max }-\left(w_{\max }-w_{\min }\right)\cdot \left(t/T\right)$. (11)

2) 非线性递减权重法

一般惯性权重采用线性递减的方法，但是运算过程的变化是非线性，二者无法很好地匹配，所以，将惯性权重非线性递减，非线性递减采用常用的指数递减法

$w=w_{\min }+\left(w_{\max }-w_{\min }\right)\cdot \exp \left(-20{\left(t/T\right)}^6\right)$. (12)

3) 基于惯性权重改进学习因子

在粒子搜索过程中，初期粒子速度大些，便于粒子进行全局搜索，避免陷入局部最优，后期粒子速度小些，便于局部搜索找到精确位置。根据公式(6)，让学习因子随着惯性权重变化而变化，在搜索前期，让学习因子c1取较大值，c2取较小值，粒子群向自身最好位置pBest学习，可以增强全局搜索能力，搜索后期，让学习因子c1取较小值，c2取较大值，粒子向全局最好位置gBest学习，可以增强局部搜索能力。现改进学习因子如下，并使其与惯性权重呈非线性函数关系，惯性权重采用指数递减法 [9]

$c1=0.5w^2+w+0.5$, (13)

$c2=2.5-c1$. (14)

粒子群优化算法的具体运算流程如图2。

5. 实验仿真与分析

为了验证改进后的PSO算法的优化性能，利用MATLAB设计一个低通IIR数字滤波器和带通IIR数字滤波器。对于仿真实验设计参数如下：群体大小为50，参数维数为12，最大加权因子 $w_{\max }=0.9$，最小加权因子 $w_{\min }=0.4$，滤波器系数取值范围限制在[−2, 2]，最大速率 $V_{\max }=4$。

Figure 2. Flow chart of particle swarm optimization

图2. 粒子群优化算法流程图

低通IIR数字滤波器的技术指标为：

$H\left({\text{e}}^{j\omega }\right)=\{\begin{array}{cc}\begin{array}{l}1\\ 0.5\\ 0\end{array}& \begin{array}{l}0\le \omega \le 0.4\pi \\ \omega =0.45\pi \\ 0.5\pi \le \omega \le \pi \end{array}\end{array}$. (15)

带通IIR数字滤波器的技术指标为：

$H_d\left({\text{e}}^{j\omega }\right)=\{\begin{array}{cc}\begin{array}{l}0\\ 1\\ 0.5\end{array}& \begin{array}{l}0\le \omega \le 0.28\pi ,0.72\pi \le \omega \le \pi \\ 0.32\pi \le \omega \le 0.68\pi \\ \omega =0.30\pi ,\omega =0.70\pi \end{array}\end{array}$. (16)

设计滤波器的二阶节个数为3个，频域采样点数为21，算法迭代次数为500次。循环运算10次，得到误差平均值。

图3画出了两种惯性权重递减策略函数随迭代次数的变化曲线，在方法三中，指数递减惯性权重在初始阶段的下降速度较小，这样使得学习因子变化幅度较慢，有利于粒子遍历整个搜索空间；随着迭代次数的增加，中期阶段的下降速度增大，使得粒子有能力向全局最优位置的周围区域逼近；到了迭代后期，惯性权重的下降速度变缓，使得粒子有足够时间在局部范围内搜索全局最优解。

Figure 3. Flow chart of particle swarm optimization

图3. 粒子群优化算法流程图

两种IIR数字滤波器的仿真试验结果如表1和表2。比较每个表中三种方法的误差和误差平均值可以看出，方法三的误差和误差平均值在三种方法中最小，优化结果最优，且方法三的优化结果优于方法二，方法二优于方法一；对比两种数字滤波器发现，相较于低通IIR数字滤波器，带通IIR数字滤波器的误差和误差平均值均偏大，带通的误差和误差平均值是低通的5倍~8倍，尤其是惯性权重和学习因子相关的方法三，带通的误差和误差平均值是低通的7倍~8倍，这说明三种方法更适合低通IIR数字滤波器的优化设计。

Table 1. Lowpass IIR filter simulation design results

表1. 低通IIR滤波器仿真设计结果

Table 2. Bandpass IIR digital filter simulation design results

表2. 带通IIR数字滤波器仿真设计结果

图4和图5是算法求解得到的两种滤波器的迭代过程曲线。从下图中可看出，两种滤波器的迭代过程曲线变化趋势相似，方法二在算法早期迅速趋近最优解，收敛速度快，方法三在早期较为缓慢，但在算法后期，经过较少次迭代后，方法三很快达到最优解，且优化结果优于其它两种方法。

Figure 4. The iterative process of the three improvement methods for lowpass IIR digital filter

图4. 低通IIR数字滤波器的三种改进方法的迭代过程

Figure 5. The iterative process of the three improvement methods for bandpass IIR digital filter

图5. 带通IIR数字滤波器的三种改进方法的迭代过程

方法一和方法二在算法早期，惯性权重都递减，算法全局最优性较强，但是两者比较之下，方法二的惯性权重的递减速度更快，算法的早期收敛速度也更快，所以应尽快让算法进入局部搜索，才能获得更好的求解效率。虽然许多调整惯性权重和学习因子的方法能够有效改善PSO算法性能，但是两者的调整方法相互脱离，一定程度上削弱了算法进化过程中的统一性，不利于算法的优化搜索。方法三中，学习因子与惯性权重相互关联，学习因子随着惯性权重的变化而变化，使其在迭代的不同时期取值不同，二者融合性能明显增强，算法的全局搜索能力有所改善，IIR数字滤波器的设计质量有所提高，更加趋近于理想幅频特性。

6. 结论

针对IIR数字滤波器逼近吸收传递函数频率响应问题，将一种基于惯性权重与学习因子相统一的改进后的粒子群优化算法应用于IIR数字滤波器设计。粒子群优化算法的运行速度也比较快，运行时间短，通过最小均方误差法设计适应度函数，在每一次迭代过程中，将惯性权重指数递减，且学习因子随之变化，两者相互统一，融合性能明显增强，改善算法的全局搜索能力，避免迭代陷入局部最优解，经过多次迭代后可以实现更小的误差逼近。利用改进后的粒子群优化算法，设计的IIR数字滤波器更加精确。

改进后的粒子群优化算法是一种通用的滤波器设计算法，它可以应用于设计低通、高通、带通、带阻滤波器，还可以用于设计多带滤波器，数字滤波器的设计决定着主动吸收系统中吸收控制器的准确度，脉冲响应函数的近似程度决定了二次反射波的吸收效果，改进后的粒子群优化算法在数值模拟中表现出优良的滤波效果。

在接下来的研究中，将考虑把更多变化的学习因子和惯性权重的方法结合到一起，找到一种更有效的新型全局粒子群优化算法，通过控制惯性权重和学习因子之间的相互关系来改善算法，增强两者之间的相互作用，分析惯性权重和学习因子的融合性能，使得它们在粒子的搜索过程中让粒子具有更强的局部和全局搜索能力，平衡算法的全局探索和局部开发能力，通过对IIR数字滤波器的优化分析，找到更优的粒子群优化算法，并将其应用于海洋模型实验中，探究其在实际试验中的滤波效果。

基金项目

中央高校基本科研业务费项目(DUT20LAB117)。