面向工业自组织网络混响衰落信道模拟的双目标优化方法*

杨璇，祁伟，秦飞

（中国科学院大学电子电气与通信工程学院，北京 101408）

0 引言

随着通信技术的发展和社会生产需要，海量的网络设备需要更宽广的无线网络覆盖面积，更高效的网络连接和网络可靠性，随着这一发展趋势，无线赋能的工业应用为工业环境和工业系统提供高速，可靠，高覆盖率的网络连接，支持工业4.0 的战略部署要求[1]。工业无线网络通常构建传感器和控制器间的闭环反馈系统，来提高工业控制系统的感知边界和控制效能[2]，因此工业无线网络通常保障闭环反馈系统下传感器和控制器之间的实时、稳定、准确的信息交互。然而，工业环境的特殊性使得工业无线网所面临的挑战与传统的无线网络有所不同[3]。工业现场环境复杂，设备密集，信号传播路径复杂，非平稳多径混响效应明显，物理层可靠传输覆盖范围有限。因此工业场景下通常采用自组织网络以提高可靠通信性能。然而工业信道的复杂性也对自组网的设备和算法提出了新的要求[4]。且不同于蜂窝网络测试，工业场景下的现场测试通常代价较大，为了研究和评估工业移动自组网[5]在真实信道下的性能和行为，同时节省测试成本和时间，可以使用信道模拟器模拟各种信道特性，以提供可复现且有真值的测试环境模拟，从而有效的提升工业移动自组网的设备研制和系统测试效能[6][7]，高效验证其在复杂环境下的可靠性和效能[8]。

目前，商用信道模拟器的通常面向室外移动场景研制，由图1(a)可以看到此场景下的多径信道呈现稀疏离散特征。因此，现有商业信道模拟器通常采用正弦叠加法研制[9][10]，其结构本质即为FIR 滤波器，如图1(b) 所示。然而与室外场景相比，工业现场存在海量金属体表面，发射机辐射信号将有较大概率与被这些金属表面多次交互，形成混响效应，从而在信道冲激响应上呈现出如图1(c) 所示连续分布，显著异于图1(a) 所示的经典模型[11][12]。此时，虽然可以依然采用FIR 滤波器架构对该系统进行模拟，但是显然将需要极高的滤波器阶数以满足模拟误差要求，给硬件设计部署带来了极大的挑战。因为IIR（Infinite Impulse Response，无限冲激响应）滤波器有递归的特性，对于同一系统，IIR 结构的滤波器的复杂度远远低于FIR结构滤波器。因此在理论上可以使用图1(d) 的较低阶IIR结构来实现工业混响信道模拟，以减少模型的复杂度和计算开销，同时尽可能地保持原始信道特性。

图1 不同场景的信道冲激响应和对应的滤波器结构

为达到该目标，一个直观的方法是首先根据采样定理将混响信道模型抽样成一个高阶FIR 滤波器架构，再利用经典的平衡截断法将高阶FIR 滤波器转换为较低阶的IIR 滤波器[13]。但是该方案首先在截断时必然面临信息损失，且无法补偿。其次该方法也缺乏在阶数和模拟误差间调整的机制，无法满足工业混响信道的模拟要求。与此同时，多项式计算法可以在给定目标冲激响应函数与阶数后，自适应迭代寻找到最优的滤波器系数组。但是该方法需要人为给定阶数作为超参，才能寻找到最优参数组，因此亦无法满足工业混响信道的模拟要求。针对这一问题，本文通过构建包含多目标优化函数，同时考虑滤波器的精度和滤波器的复杂度，将该问题建模为凸优化求最优值的问题。具体的，本文先使用平衡截断法得到一个次优的初始阶数，再通过遗传算法迭代滤波器的阶数，通过多项式计算法得出该阶数下的滤波器系数，不断更新目标函数值，搜索到全局最优值。数值实验验证了该方案的有效性。

1 现有的方案

1.1 平衡降阶法

平衡截断法[13]是1991 年由Moore 和Grimble 提出的一种降阶技术，其主要原理是通过保留原系统动态行为中的主要特征，将高维系统降至较低的维度，从而减少系统计算量和复杂度，同时保证系统的控制性能。在平衡截断法中，系统的动态特性和控制性能是通过对系统传递函数矩阵的奇异值分解（Singular Value Decomposition，SVD）来分析的。通过分析系统的左奇异向量和右奇异向量，可以得到系统的平衡矩阵，并基于此矩阵进行模型简化。但是平衡截断法也存在一些不足，为实现合理的降阶，平衡截断法忽略了系统高阶的状态信息，可能会导致降阶后的模型丧失一定的表现能力，无法通过模型完全描述系统的行为特征。并且平衡截断法在对原始系统进行截断时，并不能选择合适的截断长度来保证系统的准确性。

1.2 多项式计算法

多项式系数计算法[14]基于最小二乘法，通过最小化滤波器的频率响应与目标响应之间的误差平方和，来求解多项式系数。多项式系数计算法根据滤波器的阶数使用多项式比的形式来拟合原始滤波器的频率响应。多项式系数计算算法是一种高效而准确的方法，它可以方便地根据已知的频率响应计算出滤波器的多项式系数。该算法克服了直接从频率响应反推滤波器多项式系数的困难，提供了一个便捷和灵活的工具。但是数字滤波器的频率响应与它的多项式系数之间存在映射关系。频率响应是一个连续的函数，它的不同频率处的值于频率是连续的，并且随着频率的变化而缓慢变化。而多项式系数则是一组离散的参数，用来描述数字滤波器的性能，多项式系数的微小变化可以引起频率响应的剧烈变化，这种映射关系是一个高度非线性的过程。多项式计算法通过最小二乘法可以得到初步的近似解，但迭代计算的过程中目标函数可能存在非线性特性，因此采用高斯牛顿法来进一步优化参数[15]。高斯牛顿法[20]通过计算目标函数的梯度来更新参数，迭代进行优化，直到满足收敛条件，得到当前阶数下的与原始滤波器误差最小的滤波器的系数。但是这种方法同时也存在不足：多项式计算法不能选择准确的迭代次数，这种方法在求解多项式系数时需要通过迭代来逼近最优解，对于某些复杂的滤波器设计问题，可能需要进行大量的迭代才能达到满意的结果。因此，如何选择合适的迭代次数是该方法中需要注意的核心问题。

2 问题的构建以及求解

2.1 问题构建

通用信道模拟器主要采用正弦叠加法进行模拟，在对工业场景进行模拟时，需要采用高阶FIR 系统，该系统可以描述为：

其中y(n)是输出变量，u(n)是输入变量，bk是冲激响应的一个抽头，N表示冲激响应的阶数。根据经验知识和实际的工业混响衰落信道的参数，工业无线信道对应的离散化后的冲激响应中的阶数将会非常大，通常是上百阶的，也使得简单应用FIR 方案进行的工业连续混响信道的模拟会导致极大的硬件代价和算力要求，仅具备理论上的可行性。

现有的平衡截断法难以选择一个合适的长度来保证准确度，而多项式计算法为了找到准确率最高的解需要遍历阶数。为了解决上述两种方法的不足之处，如前文所述，在模型复杂度和模型对原始FIR 滤波器的描述能力之间寻求最佳平衡，需要选择一个较优的模型，以在保证仿真的滤波器准确性的同时考虑滤波器的阶数。理论上，随着IIR 滤波器阶数的增加，模型的精度也会提高，与原始FIR 滤波器的误差越来越小，这个关系应该是一个单调函数，如图2(a) 所示。然而，为了避免精度过高所带来的模型复杂度过大，将在这个单调函数中考虑模型的复杂度，即IIR 滤波器的阶数，并将其进行凸化处理。通过权衡模型复杂度和信道误差，构建一个目标函数，保证找到具有全局最优解的滤波器。本文设计目标函数如下：

图2 误差E和目标函数仿真示意图

其中，E是仿真后的IIR 滤波器与原FIR 滤波器的误差，a和na分别是IIR 滤波器前向系数的系数和阶数，b和nb是IIR 滤波器后向系数的系数和阶数。该表达式综合考虑了降阶和误差对结果的影响。然而，由于误差值本身较小，需要增加误差在目标函数中的权重。如图2(b)所示，可以看到目标函数此时呈现凸函数的形式，即存在一个唯一的全局最小值。下一节中将给出该全局最优解的求解方法。

2.2 问题求解

（1）初始定阶

根据上述滤波器结构，其传递函数H(z)可以表示为：

其中，h[n] 为离散的信道冲激响应，h[n]=bk,0≤n≤M。参考[16]所提方法，该滤波器也可以利用状态空间方程表示如下：

其中x(n) 是状态变量，y(n) 是输出变量，u(n) 是输入变量。根据[13]，一个线性时不变系统{A,B,C,D}，当且仅当其可控Gram 矩阵P正定时，系统是完全可控的，当且仅当Gram 矩阵为Q正定时，系统是完全可观的，当系统完全可控和完全可观时，该系统是最小实现的。其中可控矩阵P和可观矩阵Q分别是李雅普诺夫方程（Lyapunov）的解。

为了求得可控矩阵P和可观矩阵Q，对{A,B,C,D}进行平衡变换，即使用一个非奇异的变换矩阵T，将原始的线性系统{A,B,C,D}变为等价系统，其中相应的可控和可观矩阵均为对角阵并且,这里的是公式的解。这时的变换T 便称平衡变化，相应的称为平衡系统。并且平衡变换后的系统仍是最小实现的，经平衡变换后的可控矩阵和可观矩阵是正定的对角阵。

系统的状态空间达到平衡：

构建Hankel 奇异值矩阵[19]，其对角元素为系统可控矩阵和可观矩阵乘积的特征值的平方根，即有。直观意义上，Hankel 奇异值向量内每一元素代表了该目标系统模拟能力的贡献。因此对平衡系统进行截断，也就是截掉Hankel 矩阵中对角元素较小的奇异值所对应的状态，即可得到保持系统大部分特性的降阶系统[16]。显然，该方法的截断过程是一个有损过程，且无法选择合适的阶数来保证精度。所以本文只依靠平衡截断法完成初始定阶，降低后续迭代求解的开销。

（2）误差E 的求解

在得到初始次优阶数后，在算法第二步中采用多项式系数计算法得到的滤波器的误差。由于z域中z=ejω，为了方便的对ω进行求和，将其写频域的形式[14]。多项式系数计算法需要根据目标频率响应和初步定阶的阶数，写出用多项式比形式表示的滤波器：

其中原始频率响应:

两个频率响应G(jω) 和F(jω) 之间的数值差表示拟合中的误差：

将方程两边乘D(ω) 得到：

公式(12)实数和虚数的函数，可以将他们分离得到：

该函数的绝对值：

在不同频率值下表示：

定义|D(ωk)ε(ωk)|2为E，然后在不同频率上求和：

未知多项式系数ai和bi可以在最小化函数E上求得。使用最小二乘法将多项式比形式的频率响应和原始频率响应的误差平方和最小化，得到滤波器的差分方程，从而求解滤波器系数的初值。注意，最小二乘法取得最优解的前提条件是误差满足高斯分布[18]，在此场景下并不成立。因此需要根据高斯牛顿法迭代更新滤波器的系数，直至满足收敛条件，找到当前阶数下与原滤波器频率响应误差最小的滤波器系数,进而得到目标函数中的E、na和nb，代入目标函数，可以得到目标函数的解。

（3）基于遗传算法的最优解迭代

注意到目标函数中拟合误差项E 本身是一个关于滤波器阶数na和nb的表达式，但由于尚无法直接推导出E关于na和nb的解析式，因此难以构建解析法来求解目标方程式的最小值。但是通过初始定阶定位至最优解附近后，可以通过暴力搜索法来搜索最优解。显然，可以通过应用遗传算法等成熟方法来加速暴力搜索过程遗传算法[17]是一种模拟自然生物进化过程的优化算法，它通过模拟自然选择、遗传、变等过程，对变量空间中的解进行优化搜索，从而加速找到最优解的过程。

该算法第一步通过随机选择10 个不同滤波器阶数作为一个种群，第二步通过计算种群中每个个体对应的目标函数的值来评价个体适应度并选择适应度高的参数组保留，第三步对目标函数值小的个体赋予更大的保留概率，对目标函数值大的个体赋予小的保留概率.第四步是两个个体一组基于概率选择交叉点，将两个个体交叉点后的阶数进行互换，例如两个个体(a,b) 和(c,d)，如果交叉点是第一位，则两个个体变成(a,d) 和(b,c)。第五步基于设定变异率对保留下的种群进行变异，将个体(a,b) 随机变异为(a±range,b±range) 范围内的某一个体(a1,b1)，最后生成新的种群，再重复第二到第四步,若迭代的过程平均适应度变化率小于设定门限即认为满足收敛条件[21]，即找到目标函数值最小的个体，也就是问题的最优解。

3 实验验证

本节中将对所提工业混响信道模拟的多目标优化方法进行仿真验证。目标信道模型采用如图3 所示的在中国科学院大学微电子工艺间实测所得混响信道模型，并采用如文献[16] 所述方法对混响信道模型进行抽样后得到一个300 阶的FIR 滤波器。

图3 无线信道冲激响应

经过平衡截断法初始定阶后，系统初始阶数为na=33,nb=33。此时目标函数值为613.219 5。随机构建种群为[(24,31),(34,26),(36,42),(33,30),(42,38),(35,38),(41,41),(26,35),(26,38),(39,30)]。如图4 所示，经过一次迭代后的新种群为[(33,38),(26,26),(41,30),(39,41),(41,41),(41,41),(34,30),(39,37),(41,41),(42,41)]，此时目标函数值为91.3。经过4 次迭代后，目标函数值变化为67.3，此时种群为[(35,30),(26,41),(41,41),(41,28),(42,30),(34,41),(26,30),(26,30),(26,30),(25,30)]。由于变化率在第二次到第三次迭代和第三次到第四次迭代均小于0.1 %，算法即认为寻找到系统最优解，中止迭代，此时系统阶数为(26,30)，目标函数值为67.3。与暴力搜索相比显著减少了计算开销，提高算法效率。

图4 遗传算法演化图

图5 给出了本文所提算法的收敛结果。与平衡截断法得到的(30,30) 阶的滤波器在幅值响应和频率响应相比，本文所提算法得到的冲激响应与原来的冲激响应相比，明显具有较高的拟合度。平衡截断法得到的(30,30)阶RMSE 为10.87%，而本文得到的(26,30) 阶的RMSE仅为1.13%，实现了数量级的性能提升。而与多项式计算法得到的(30,30) 阶的滤波器在幅值响应和频率响应进行对比，多项式计算法的RNMSE 则为2.18%，两者拟合度接近。但是本文所提方法明显通过更低的滤波器阶数实现了更高的拟合度。

图5 两种方法与原始冲激响应的幅值响应和相位响应对比图

4 结束语

本文提出的一种双目标优化方法具备用低阶IIR 系统模拟工业混响信道连续冲激响应的能力，在保证系统精度的同时，大大降低了信道模拟系统的复杂度。实验验证表明，与现有的平衡截断法相比，本文所提方法不但成功收敛到系统最优阶数，且模拟精度的RNMSE 仅为1.13%，远远低于平衡截断方法RNMSE 的10.87%，达到了数量级的性能提升。在下一步工作，我们将进一步构建统一的优化迭代模型，合并分步迭的过程，以进一步提高工业混响信道的模拟效能。