通过特征选择实现数据降维，构建结构精简的辨识模型，能够有效避免输入特征太多造成“维数灾难”以及给学习模型带来“过拟合”等问题^[1-3]。目前，特征选择的方法主要有偏最小二乘回归(partial least squares regression, PLSR)^[4]、灰色关联分析(grey relational analysis, GRA)^[5]和互信息(mutual information, MI)^[6-7]等。MI对样本的分布类型无特别要求，可有效捕捉特征间的非线性关系，特别适合多元序列特征选择问题。

互信息是一种衡量两个随机变量之间相互依赖强弱程度的准则，其来源于信息论中熵的概念^[8]。采用互信息进行特征选择主要有以下两个难点:1)特征评价策略；2)互信息的准确估计。针对评价策略，Battiti等^[9]提出了MIFS(mutual information feature selection)方法，通过平衡相关特征和冗余特征的选择，取得了较好的效果，此后，Peng等^[10]、Fleuret^[11]、Yang等^[12]、Brown^[13]和韩等^[14]都分别提出了各自的评价策略。上述方法主要是采用一维互信息来衡量特征好坏，当存在冗余特征时，并不能保证得到最优特征子集。互信息的估计一般通过非参数方式^[15-16]求解概率密度函数来实现，主要有直方图法和核函数法。近几年发展的k-近邻法(k-nearest neighbors, kNN) ^[17-18]构建的数据驱动型特征选择框架，无需假设概率密度函数形式，避免了对概率密度函数的估计，其基础理论完备，非常适合具有非线性不规则分布特点的实际数据，但是高维k-近邻互信息的估计存在一定困难。

因此，本文在k-近邻互信息基础上将其扩展用于高维互信息的估计，然后采用前向累加后向交叉(forward accumulation and backward cross, FABC)的最优特征子集选择策略，构建kNN-FABC的特征选择方法，最后以Friedman数据、Housing数据和实际污水处理出水总磷预测数据为例进行仿真实验，并结合多层感知器(multilayer perceptron, MLP)预测模型来验证所提特征选择方法的有效性。

1 互信息理论

互信息是一种衡量两个随机变量之间相互依赖强弱程度的准则，其来源于信息论中熵的概念。两个随机变量之间的互信息越大，则两者之间的相关性就越强^[17]。

对随机变量X和Y来说，设其联合概率密度函数为ρ_{X, Y}(x, y)，则其边缘概率密度函数为

$ {\rho _X}\left( x \right) = \int {{\rho _{X,Y}}\left( {x,y} \right){\rm{d}}y} $

(1)

$ {\rho _Y}\left( y \right) = \int {{\rho _{X,Y}}\left( {x,y} \right){\rm{d}}x} $

(2)

根据信息论的有关理论^[17]，随机变量X和Y之间的互信息定义为

$ I\left( {X;Y} \right) = \iint\limits_S {{\rho _{X,Y}}\left( {x,y} \right)\log \frac{{{\rho _{X,Y}}\left( {x,y} \right)}}{{{\rho _X}\left( x \right){\rho _Y}\left( y \right)}}{\rm{d}}x{\rm{d}}y} $

(3)

式中S为X和Y的定义域。

根据Shannon熵的定义^[17]，X、Y的熵和联合熵分别为

$ H\left( X \right) = - \int {{\rho _X}\left( x \right){\rm{log}}{\rho _X}\left( x \right){\rm{d}}x} $

(4)

$ H\left( Y \right) = - \int {{\rho _Y}\left( y \right){\rm{log}}{\rho _Y}\left( y \right){\rm{d}}y} $

(5)

$ H\left( {X,Y} \right) = - \iint\limits_S {{\rho _{X,Y}}\left( {x,y} \right)\log {\rho _{X,Y}}\left( {x,y} \right){\rm{d}}x{\rm{d}}y} $

(6)

根据式(4)、(5)和(6)，式(3)还可以表示为

$ I\left( {X;Y} \right) = H\left( X \right) + H\left( Y \right) - H\left( {X,Y} \right) $

(7)

2 高维k-近邻互信息理论 2.1 高维k-近邻互信息的估计

Kraskov等^[17-18]提出的k-近邻法通过数据直接近似估计特征的互信息，避免了对概率密度分布的近似估计。其基本思想是，在由随机变量X和Y构成的空间中对于给定的N个样本首先找出它的k个近邻，再找出其他样本分别在X和Y方向到当前样本的距离小于当前样本到k个近邻距离的最大值的数目，通过统计数目从而进行互信息的估计。

现将其扩展到高维互信息的估计，高维互信息计算的具体思路如下^[18]：

给定N个样本z⁽ⁱ⁾=[(x⁽ⁱ⁾)^T(y⁽ⁱ⁾)^T]^T (i=1, 2, …, N)，其中x⁽ⁱ⁾∈$\mathbb{R}$^d且y⁽ⁱ⁾∈$\mathbb{R}$^p。如果z和z′为数据集中两个不同的向量，则

$ \left\| {\mathit{\boldsymbol{z}} - \mathit{\boldsymbol{z'}}} \right\| = \max \left( {\left\| {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{x'}}} \right\|,\left\| {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{y'}}} \right\|} \right) $

(8)

式中‖·‖表示取最大范数。

若z⁽ⁱ⁾的k-近邻为z^(k(i))=[(x^(k(i)))^T(y^(k(i)))^T]^T，则z⁽ⁱ⁾与z^(k(i))之间的Euclidean距离为

$ \begin{array}{*{20}{c}} {{\varepsilon ^{\left( i \right)}} = \left\| {{\mathit{\boldsymbol{z}}^{\left( i \right)}} - {\mathit{\boldsymbol{z}}^{\left( {k\left( i \right)} \right)}}} \right\| = }\\ {\max \left( {\left\| {{\mathit{\boldsymbol{x}}^{\left( i \right)}} - {\mathit{\boldsymbol{x}}^{\left( {k\left( i \right)} \right)}}} \right\|,\left\| {{\mathit{\boldsymbol{y}}^{\left( i \right)}} - {\mathit{\boldsymbol{y}}^{\left( {k\left( i \right)} \right)}}} \right\|} \right)} \end{array} $

(9)

对于z⁽ⁱ⁾中的分量x⁽ⁱ⁾和y⁽ⁱ⁾有

$ \varepsilon _x^{\left( i \right)} = \left\| {{\mathit{\boldsymbol{x}}^{\left( i \right)}} - {\mathit{\boldsymbol{x}}^{\left( {k\left( i \right)} \right)}}} \right\| = \mathop {\max }\limits_{s = 1,2, \cdots ,d} \left| {\mathit{\boldsymbol{x}}_s^{\left( i \right)} - \mathit{\boldsymbol{x}}_s^{\left( {k\left( i \right)} \right)}} \right| $

(10)

$ \varepsilon _y^{\left( i \right)} = \left\| {{\mathit{\boldsymbol{y}}^{\left( i \right)}} - {\mathit{\boldsymbol{y}}^{\left( {k\left( i \right)} \right)}}} \right\| = \mathop {\max }\limits_{t = 1,2, \cdots ,p} \left| {\mathit{\boldsymbol{y}}_t^{\left( i \right)} - \mathit{\boldsymbol{y}}_t^{\left( {k\left( i \right)} \right)}} \right| $

(11)

将与x_s⁽ⁱ⁾和y_t⁽ⁱ⁾之间距离小于ε⁽ⁱ⁾的样本点个数分别记为N_xs⁽ⁱ⁾和N_yt⁽ⁱ⁾。因此，变量X和Y之间互信息的估计量为

$ \begin{array}{*{20}{c}} {\hat I\left( {X;Y} \right) = \mathit{\Psi }\left( k \right) - \frac{{d + p - 1}}{k} + \left( {d + p - 1} \right)\mathit{\Psi }\left( N \right) - }\\ {\frac{1}{N}\sum\limits_{p = 1}^N {\left( {\sum\limits_{s = 1}^d {\mathit{\Psi }\left( {N_{{x_s}}^{\left( i \right)}} \right)} + \sum\limits_{t = 1}^p {\mathit{\Psi }\left( {N_{{y_t}}^{\left( i \right)}} \right)} } \right)} } \end{array} $

(12)

式中：Ψ(·)为digamma函数，Ψ(t+1)= Ψ(t)+ $\frac{1}{t}$ ，Ψ(1)≈-0.577 215 6；k一般取2~6。通过式(12)可以发现， $\hat I\left( {X,Y} \right)$ 不仅与k和N有关，而且与变量X和Y的维数有关，能够实现高维互信息的估计。

2.2 高维k-近邻互信息特征选择流程

由于输入特征之间并非互相独立，存在非线性耦合，同时根据信息论知识，增加相关特征，能更好表征输出特征。基于高维k-近邻互信息特征选择的基本思想为，首先利用前向累加搜索策略找出由强相关特征和弱相关特征组成的相关特征候选子集，然后采用后向交叉搜索策略剔除候选子集中的冗余特征，最终形成最优强相关特征子集。具体的流程如下。

1) 参数初始化，设置k值和无关特征个数。

2) 根据式(12)计算输入特征X的每一维分量与输出特征Y之间的一维互信息，提取出互信息最大值对应的特征，作为第1个相关特征。

3) 在2)得到的特征基础上，利用式(12)计算其他输入特征与输出特征Y之间的高维互信息，找出第2个相关特征。

4) 重复2)~3)，直至所有特征处理完毕，得到全部特征变量的排序，根据预设无关特征个数，剔除无关特征，得到相关特征子集。

5) 根据相关特征子集，计算两两之间的互信息，找出互信息最大值对应的特征组。

6) 结合2)的一维互信息值，剔除冗余特征，最终得到最优强相关特征子集。

3 仿真实验

为了验证本文特征选择方法的有效性，本文分别采用标准Friedman函数数据和实际污水处理出水总磷预测数据进行特征选择的实验分析，并与PLSR、MRMR和JMI等方法进行比较，最后利用MLP预测模型^[19]进行误差分析。实验时，k设置为5，MLP预测模型的隐含层节点取20个，迭代次数为20 000次，学习率η取为0.001。由于神经网络建模受初值影响具有随机性，本文实验结果为运行30次取平均值和标准差。

3.1 Friedman数据

Friedman由式(13)描述：

$ \begin{array}{*{20}{c}} {Y = 10\sin \left( {\pi {\mathit{\boldsymbol{x}}_1}{\mathit{\boldsymbol{x}}_2}} \right) + 20{{\left( {{\mathit{\boldsymbol{x}}_3} - \frac{1}{2}} \right)}^2} + 10{\mathit{\boldsymbol{x}}_4} + }\\ {5{\mathit{\boldsymbol{x}}_5} + 0 \times \sum\limits_{i = 6}^{10} {{\mathit{\boldsymbol{x}}_i} + {\mathit{\boldsymbol{x}}_{11}} + {\mathit{\boldsymbol{x}}_{12}} + \varepsilon } } \end{array} $

(13)

式中：x₁~x₁₀服从[0, 1]的均匀分布，x₁₁=0.5x₁+ ε，x₁₂=0.5x₂+ε。随机产生500个含高斯噪声数据作样本，且高斯噪声ε满足N(0, 0.1)。由式(13)可知，Y只与x₁~x₅有关，x₆~x₁₀是无关特征，x₁₁和x₁₂属于冗余特征。

利用本文方法进行特征选择，第1步得到的互信息如图 1所示，可见仅依靠一维互信息会误选x₁₁，误剔除x₃和x₅；第2步得到的高维互信息如图 2所示，设定无关变量个数为5，则剔除x₆~x₁₀得到最优相关特征子集为x₁~x₅、x₁₁和x₁₂；第3步得到的互信息如图 3所示，确定x₁、x₁₁和x₂、x₁₂中存在冗余特征，再结合一维互信息值，可以判定x₁₁和x₁₂为冗余特征，剔除后得到最优强相关特征子集为x₁~x₅，完全与式(13)特征属性相吻合。

利用本文特征选择方法，并利用MLP网络对Friedman数据进行测试，随机产生240个带噪声数据作训练样本，1 000个不含噪声数据作测试样本，变量选择及网络测试结果如表 1所示。

由表 1可知，kNN-FABC方法和文献[17]的STEP方法都成功筛选出全部仅有的5个相关特征，MRMR、CMIM和MIFS误剔相关特征x₃、误选无关特征x₆，JMI和MMI误剔相关特征x₃、误选冗余特征x₁₂，PLSR误选冗余特征x₁₁和x₁₂。同时，kNN-FABC方法的RMSE均值和标准差都要远小于其他方法，表明其极大地提高了网络泛化能力。

3.2 Housing数据

Housing数据由13个房屋属性的输入特征和1个房屋价格的输出特征组成，共有506组样本。随机选择338组作训练样本，剩余168组作测试样本。利用本文方法进行特征选择，并与其他文献特征选择结果进行了比较，特征选择及MLP网络测试结果如表 2所示。

由表 2可知，在归一化方法、MLP模型参数选取一致的情况下，kNN-FABC的RMSE均值和标准差均表明了采用本文方法所选特征建立的预测模型精度更高。

3.3 污水处理出水总磷预测数据

近年来，总磷(total phosphorus, TP)污染问题开始凸显。研究表明，水体富营养化程度与总磷等指标密切相关。因此加强对TP的检测有助于解决水体富营养化问题^[20]。目前，污水处理厂普遍采用的生化方法成本高、耗时长，在线检测仪表对检测环境要求高并且维护成本昂贵，而基于数据驱动的软测量技术能够在线、快速、准确预测，受到研究者青睐。但是污水处理生化反应过程中可测变量较多，在建立TP软测量模型时，若选取全部可用特征，则会产生冗余，降低网络泛化能力，同时也大幅增加检测成本，故对TP预测模型进行输入特征选择具有重要意义。污水处理生化反应过程中可测变量及其意义如表 3所示。本文采用kNN-FABC方法实现对TP预测特征变量的选择。

从北京市某污水处理厂获取了2015年6~8月共492组数据，首先剔除明显异常值，然后采用小波包(sym8小波、2层分解、软阈值方式)进行降噪处理，处理效果见图 4。

从原始数据集中每隔3个样本选取1个样本共123个组成测试集，其余369个作为训练集。TP数据特征选择及MLP网络测试结果如表 4所示。相比PLSR，kNN-FABC获得了较少的特征、较低的RMSE及其标准差，寻找到了最有效的特征子集，构建了更稳定的软测量模型，验证了kNN-FABC方法在特征选择上的有效性。

图 5和图 6分别给出了MLP模型对TP数据预测结果和预测误差。可见，通过kNN-FABC方法选出的特征能够很好地表征原始数据，建模误差主要集中在-0.3~0.3 mg/L之间，能够较好地满足污水处理厂对总磷检测精度的要求。

4 结束语

通过将数据驱动型k-近邻方法扩展用于多维特征之间相关性计算，结合前向累加策略，能够较为准确地获得全部特征的最优排序，进而剔除无关特征，再结合后向交叉策略，能够定位并删除冗余特征，最终得到最优强相关特征子集。Friedman、Housing和实际污水处理出水TP预测实验证实了本文方法的有效性。下一步的研究工作是实现无关特征的自动判定，避免人为设置带来的无关特征误剔的可能。