文章目录
一、理论基础1、BP神经网络概述2、RBF神经网络概述二、案例背景1、问题描述2、解题思路与步骤三、MATLAB程序实现1、产生训练集/测试集2、创建/训练BP神经网络及仿真测试3、创建RBF神经网络及仿真测试4、性能评价5、绘图四、网络参数的影响与选择五、参考文献补充一、理论基础
1、BP神经网络概述
（1）BP神经网络的结构
BP神经网络由Rumelhard和McClelland于1986年提出，从结构上讲，它是一种典型的多层前向型神经网络，具有一个输入层、数个隐含层（可以是一层，也可以是多层）和一个输出层。层与层之间采用全连接的方式，同一层的神经元之间不存在相互连接。理论上已经证明，具有一个隐含层的三层网络可以逼近任意非线性函数。
隐含层中的神经元多采用S型传递函数，输出层的神经元多采用线性传递函数。图1所示为一个典型的BP神经网络结构，该网络具有一个隐含层，输入层神经元数目为mmm，隐含层神经元数目为lll，输出层神经元数目为nnn，隐含层采用S型传递函数tansig，输出层传递函数为purelin。
图1 典型的BP神经网络
（2）BP神经网络的学习算法
BP神经网络的误差反向传播算法是典型的有导师指导的学习算法，其基本思想是对一定数量的样本对（输入和期望输出）进行学习，即将样本的输入送至网络输入层的各个神经元，经隐含层和输出层计算后，输出层各个神经元输出对应的预测值。若预测值与期望输出之间的误差不满足精度要求，则从输出层反向传播该误差，从而进行权值和阈值的调整，使得网络的输出和期望输出间的误差逐渐减小，直至满足精度要求。
BP网络的精髓是将网络的输出与期望输出间的误差归结为权值和阈值的“过错”，通过反向传播把误差“分摊”给各个神经元的权值和阈值。BP网络学习算法的指导思想是权值和阈值的调整要沿着误差函数下降最快的方向——负梯度方向。
下面详细推导利用BP网络学习算法对权值和阈值进行调整的公式。
设一样本对(X,Y)(\boldsymbol X,\boldsymbol Y)(X,Y)为X=[x1,x2,⋯,xm]′\boldsymbol X=[x_1,x_2,\dotsm,x_m]'X=[x1,x2,⋯,xm]′，Y=[y1,y2,⋯,yn]′\boldsymbol Y=[y_1,y_2,\dotsm,y_n]'Y=[y1,y2,⋯,yn]′，隐含层神经元为O=[O1,O2,⋯,Ol]\boldsymbol O=[O_1,O_2,\dotsm,O_l]O=[O1,O2,⋯,Ol]。输入层与隐含层神经元间的网络权值距离W1\boldsymbol W^1W1和隐含层与输出层神经元间的网络权值W2\boldsymbol W^2W2分别为W1=[w111w121⋯w1m1w211w221⋯w2m1⋮⋮⋮wl11wl21⋯wlm1],W2=[w112w122⋯w1l2w212w222⋯w2l2⋮⋮⋮wn12wn22⋯wnl2](1)\boldsymbol W^1=\begin{bmatrix} w_{11}^1 & w_{12}^1 & \cdots & w_{1m}^1 \\ w_{21}^1 & w_{22}^1 & \cdots & w_{2m}^1 \\ \vdots & \vdots & & \vdots \\ w_{l1}^1 & w_{l2}^1 & \cdots & w_{lm}^1\\ \end{bmatrix}, \boldsymbol W^2=\begin{bmatrix} w_{11}^2 & w_{12}^2 & \cdots & w_{1l}^2 \\ w_{21}^2 & w_{22}^2 & \cdots & w_{2l}^2 \\ \vdots & \vdots & & \vdots \\ w_{n1}^2 & w_{n2}^2 & \cdots & w_{nl}^2\\ \end{bmatrix} \tag{1}W1=⎣⎢⎢⎢⎡w111w211⋮wl11w121w221⋮wl21⋯⋯⋯w1m1w2m1⋮wlm1⎦⎥⎥⎥⎤,W2=⎣⎢⎢⎢⎡w112w212⋮wn12w122w222⋮wn22⋯⋯⋯w1l2w2l2⋮wnl2⎦⎥⎥⎥⎤(1)隐含层神经元的阈值θ1\boldsymbol \theta^1θ1和输出层神经元的阈值θ2\boldsymbol \theta^2θ2分别为θ1=[θ11,θ21,⋯,θl1]′,θ2=[θ12,θ22,⋯,θn2]′(2)\boldsymbol \theta^1=[\theta_1^1,\theta_2^1,\cdots,\theta_l^1]',\boldsymbol \theta^2=[\theta_1^2,\theta_2^2,\cdots,\theta_n^2]'\tag{2}θ1=[θ11,θ21,⋯,θl1]′,θ2=[θ12,θ22,⋯,θn2]′(2)则隐含层神经元的输出为Oj=f(∑i=1mwji1−θj1)=f(netj),j=1,2,⋯,l(3)O_j=f(\sum_{i=1}^m w_{ji}^1-\theta_j^1)=f(net_j),\quad j=1,2,\cdots,l \tag{3}Oj=f(i=1∑mwji1−θj1)=f(netj),j=1,2,⋯,l(3)其中，netj=∑i=1mwji1−θj1,j=1,2,⋯,l;f(⋅)net_j=\displaystyle\sum_{i=1}^m w_{ji}^1-\theta_j^1,j=1,2,\cdots,l;f(\boldsymbol \sdot)netj=i=1∑mwji1−θj1,j=1,2,⋯,l;f(⋅)为隐含层的传递函数。
输出层神经元的输出为zk=g(∑j=1lwkj2Oj−θk2)=g(netk),k=1,2,⋯,n(4)z_k=g(\sum_{j=1}^l w_{kj}^2O_j-\theta_k^2)=g(net_k),\quad k=1,2,\cdots,n\tag{4}zk=g(j=1∑lwkj2Oj−θk2)=g(netk),k=1,2,⋯,n(4)其中，netk=∑j=1lwkj2Oj−θk2,k=1,2,⋯,n;g(⋅)net_k=\displaystyle\sum_{j=1}^l w_{kj}^2O_j-\theta_k^2,k=1,2,\cdots,n;g(\boldsymbol\sdot)netk=j=1∑lwkj2Oj−θk2,k=1,2,⋯,n;g(⋅)为输出层的传递函数。
网络输出与期望输出的误差为E=12∑k=1n(yk−zk)2=12∑k=1n[yk−g(∑j=1lwkj2Oj−θk2)]2=E=\frac12\sum_{k=1}^n(y_k-z_k)^2=\frac12\sum_{k=1}^n[y_k-g(\sum_{j=1}^l w_{kj}^2O_j-\theta_k^2)]^2=E=21k=1∑n(yk−zk)2=21k=1∑n[yk−g(j=1∑lwkj2Oj−θk2)]2=12∑k=1n{yk−g[∑j=1lwkj2f(∑i=1mwij1xi−θj1)−θk2]}2(5)\frac12\sum_{k=1}^n\{y_k-g[\sum_{j=1}^lw_{kj}^2f(\sum_{i=1}^mw_{ij}^1x_i-\theta_j^1)-\theta_k^2]\}^2\tag{5}21k=1∑n{yk−g[j=1∑lwkj2f(i=1∑mwij1xi−θj1)−θk2]}2(5)误差EEE对隐含层与输出层神经元间的权值wkj2w_{kj}^2wkj2的偏导数为∂E∂wkj2=∂E∂zk∂zk∂wkj2=−(yk−zk)g′(netk)Oj=−δk2Oj(6)\frac{\partial E}{\partial w_{kj}^2}=\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial w_{kj}^2}=-(y_k-z_k)g'(net_k)O_j=-\delta_k^2O_j\tag{6}∂wkj2∂E=∂zk∂E∂wkj2∂zk=−(yk−zk)g′(netk)Oj=−δk2Oj(6)其中，δk2=(yk−zk)g′(netk)\delta_k^2=(y_k-z_k)g'(net_k)δk2=(yk−zk)g′(netk)。
误差EEE对输入层与隐含层神经元间的权值wji1w_{ji}^1wji1的偏导数为∂E∂wji1=∑k=1n∑j=1l∂E∂zk∂zk∂Oj∂Oj∂wji1=\frac{\partial E}{\partial w_{ji}^1}=\sum_{k=1}^n\sum_{j=1}^l\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial O_j}\frac{\partial O_j}{\partial w_{ji}^1}=∂wji1∂E=k=1∑nj=1∑l∂zk∂E∂Oj∂zk∂wji1∂Oj=−∑k=1n(yk−zk)g′(netk)wkj2f′(netj)xi=−δj1xi(7)-\sum_{k=1}^n(y_k-z_k)g'(net_k)w_{kj}^2f'(net_j)x_i=-\delta_j^1x_i\tag{7}−k=1∑n(yk−zk)g′(netk)wkj2f′(netj)xi=−δj1xi(7)其中，δj1=∑k=1n(yk−zk)g′(netk)wkj2f′(netj)=f′(netj)∑k=1nδk2wkj2\delta_j^1=\displaystyle\sum_{k=1}^n(y_k-z_k)g'(net_k)w_{kj}^2f'(net_j)=f'(net_j)\displaystyle\sum_{k=1}^n\delta_k^2w_{kj}^2δj1=k=1∑n(yk−zk)g′(netk)wkj2f′(netj)=f′(netj)k=1∑nδk2wkj2。
由式（6）和式（7）可得权值的调整公式为{wji1(t+1)=wji1(t)+Δwji1=wji1(t)−η1∂E∂wji1=wji1(t)+η1δj1xiwkj2(t+1)=wkj2(t)+Δwkj2=wki2(t)−η2∂E∂wkj2=wkj2(t)+η2δj2Oj(8)\begin{dcases} w_{ji}^1(t+1)=w_{ji}^1(t)+\Delta w_{ji}^1=w_{ji}^1(t)-\eta^1\frac{\partial E}{\partial w_{ji}^1}=w_{ji}^1(t)+\eta^1\delta_j^1x_i\\ w_{kj}^2(t+1)=w_{kj}^2(t)+\Delta w_{kj}^2=w_{ki}^2(t)-\eta^2\frac{\partial E}{\partial w_{kj}^2}=w_{kj}^2(t)+\eta^2\delta_j^2O_j \tag{8} \end{dcases}⎩⎪⎪⎨⎪⎪⎧wji1(t+1)=wji1(t)+Δwji1=wji1(t)−η1∂wji1∂E=wji1(t)+η1δj1xiwkj2(t+1)=wkj2(t)+Δwkj2=wki2(t)−η2∂wkj2∂E=wkj2(t)+η2δj2Oj(8)其中，η1\eta^1η1和η2\eta^2η2分别为隐含层和输出层的学习步长。
同理，误差EEE对输出层神经元的阈值θk2\theta_k^2θk2的偏导数为∂E∂θk2=∂E∂zk∂zk∂θk2=−(yk−zk)g′(netk)(−1)=(yk−zk)g′(netk)=δk2(9)\frac{\partial E}{\partial \theta_k^2}=\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial \theta_k^2}=-(y_k-z_k)g'(net_k)(-1)=(y_k-z_k)g'(net_k)=\delta_k^2\tag{9}∂θk2∂E=∂zk∂E∂θk2∂zk=−(yk−zk)g′(netk)(−1)=(yk−zk)g′(netk)=δk2(9)误差EEE对隐含层神经元的阈值θj1\theta_j^1θj1的偏导数为∂E∂θj1=∑k=1n∂E∂zk∂zk∂Oj∂Oj∂θj1=−∑k=1n(yk−zk)g′(netk)wkj2f′(netj)(−1)=\frac{\partial E}{\partial \theta_j^1}=\sum_{k=1}^n\frac{\partial E}{\partial z_k}\frac{\partial z_k}{\partial O_j}\frac{\partial O_j}{\partial \theta_j^1}=-\sum_{k=1}^n(y_k-z_k)g'(net_k)w_{kj}^2f'(net_j)(-1)=∂θj1∂E=k=1∑n∂zk∂E∂Oj∂zk∂θj1∂Oj=−k=1∑n(yk−zk)g′(netk)wkj2f′(netj)(−1)=∑k=1n(yk−zk)g′(netk)wkj2f′(netj)=δj1(10)\sum_{k=1}^n(y_k-z_k)g'(net_k)w_{kj}^2f'(net_j)=\delta_j^1\tag{10}k=1∑n(yk−zk)g′(netk)wkj2f′(netj)=δj1(10)由式（9）和式（10）可得阈值的调整公式为{θj1(t+1)=θj1(t)+Δθj1=θj1(t)+η1∂E∂θj1=θj1(t)+η1δj1θk2(t+1)=θk2(t)+Δθk2=θk2(t)−η2∂E∂θk2=θk2(t)+η2δk2(11)\begin{dcases} \theta_j^1(t+1)=\theta_j^1(t)+\Delta\theta_j^1=\theta_j^1(t)+\eta^1\frac{\partial E}{\partial \theta_j^1}=\theta_j^1(t)+\eta^1\delta_j^1\\ \theta_k^2(t+1)=\theta_k^2(t)+\Delta\theta_k^2=\theta_k^2(t)-\eta^2\frac{\partial E}{\partial \theta_k^2}=\theta_k^2(t)+\eta^2\delta_k^2\tag{11} \end{dcases}⎩⎪⎪⎨⎪⎪⎧θj1(t+1)=θj1(t)+Δθj1=θj1(t)+η1∂θj1∂E=θj1(t)+η1δj1θk2(t+1)=θk2(t)+Δθk2=θk2(t)−η2∂θk2∂E=θk2(t)+η2δk2(11)
（3）BP神经网络的MATLAB工具箱函数
MATLAB神经网络工具箱中包含了许多用于BP神经网络分析与设计的函数，本节将详细说明几个主要函数的功能、调用格式、参数意义及注意事项等。
BP神经网络创建函数
自Rb版本以后，MATLAB神经网络工具箱对BP神经网络的创建函数进行了更新，将函数feedforwardnet用于创建一个前向传播神经网络，其调用格式为net=feedforwardnet(hiddenSizes,trainFcn)\boldsymbol {net=feedforwardnet(hiddenSizes,trainFcn)}net=feedforwardnet(hiddenSizes,trainFcn)其中，hiddenSizes为一个行向量，表征一个或多个隐含层所包含的神经元个数（默认为10，即仅有一个包含10个神经元的隐含层）；trainFcn为网络训练函数（默认为trainlm）。
可以发现，新版本的BP神经网络创建函数无需给定输入和输出向量的信息，且输入参数明显减少并均有默认值，用户调用该函数时甚至无需给出任何参数。同时，由于新版本创建的网络中无需存储输入和输出向量的信息，使得在内存管理上更加高效。
BP神经网络训练函数
函数train用于训练已经创建好的BP神经网络，其调用格式为[net,tr]=train(net,P,T,Pi,Ai,EW)\boldsymbol{[net,tr]=train(net,P,T,Pi,Ai,EW)}[net,tr]=train(net,P,T,Pi,Ai,EW)其中，net为训练前及训练后的网络；P为网络输入向量；T为网络目标向量（默认为0）；Pi为初始的输入层延迟条件（默认为0）；Ai为初始的输出层延迟条件（默认为0）；tr为训练记录（包含步数及性能）。需要注意的是，自Rb版本以后，MATLAB神经网络工具箱对train函数新增了一个输入参数EW，这意味着用户可以通过设置该参数调整输出目标向量中各个元素的重要程度。
BP网络预测函数
函数sim用于利用已经训练好的BP神经网络进行仿真预测，其调用格式为[Y,Pf,Af,E,perf]=sim(net,P,Pi,Ai,T)\boldsymbol{[Y,Pf,Af,E,perf]=sim(net,P,Pi,Ai,T)}[Y,Pf,Af,E,perf]=sim(net,P,Pi,Ai,T)其中，net为训练好的网络；P为网络输入向量；Pi为初始的输入层延迟条件（默认为0）；Ai为初始的输出层延迟条件（默认为0）；T为网络目标向量（默认为0）；Y为网络输出向量；Pf为最终的输入层延迟条件；Af为最终的隐含层延迟条件；E为网络误差向量；perf为网络的性能。值得一提的是，在新版本的神经网络工具箱中，除了上述方式，还可以做如下调用y=net(P,Pi,Ai)\boldsymbol{y=net(P,Pi,Ai)}y=net(P,Pi,Ai)
2、RBF神经网络概述
（1）RBF神经网络的结构
1985年，Powell提出了多变量插值的径向基函数（radial basis function,RBF）方法。1988年，Moody和Darken提出了一种神经网络结构，即RBF神经网络，它能够以任意精度逼近任意连续函数。
RBF网络的结构与多层前向型网络类似，是一种三层前向型网络，其网络结构如图2所示。输入层由信号源结点组成；第二层为隐含层，隐含层神经元数目视所描述问题的需要而定，隐含层神经元的传递函数是对中心点径向对称且衰减的非负非线性函数；第三层为输出层，它对输入模式的作用做出响应。从输入空间到隐含层空间的变换是非线性的，而从隐含层空间到输出层空间的变换是线性的。
图2 RBF神经网络结构
RBF网络的基本思想是：用RBF作为隐含层神经元的“基”构成隐含层空间，这样就可以将输入矢量直接映射到隐含层空间，而不需要通过权连接。当RBF的中心点确定以后，这种映射关系也就确定了。而隐含层空间到输出层空间的映射是线性的，即网络的输出是隐含层神经元输出的线性加权和，此处的权即为网络可调参数。由此可见，从总体上看，网络由输入到输出的映射是非线性的，而网络输出对可调参数而言却又是线性的。这样网络的权就可以由线性方程直接解出，从而大大加快学习速度并避免局部极小值问题。
（2）RBF神经网络的学习算法
根据隐含层神经元数目的不同，RBF神经网络的学习算法总体上可以分为两种：
<1>隐含层神经元数目逐渐增加，经过不断的循环迭代，实现权值和阈值的调整与修正。
<2>隐含层神经元数目确定（与训练集样本数目相同），权值和阈值由线性方程组直接解出。
对比不难发现，第二种学习算法速度更快、精度更高，因此下面仅讨论第二种学习算法，即隐含层神经元数目确定且等于训练集样本数目这一类型。
具体的学习算法步骤如下;
<1>确定隐含层神经元径向基函数中心
为不失一般性，设训练集样本输入矩阵P\boldsymbol PP和输出矩阵T\boldsymbol TT分别为P=[p11p12⋯p1Qp21p22⋯p2Q⋮⋮⋮pM1pM2⋯pMQ],T=[t11t12⋯t1Qt21t22⋯t2Q⋮⋮⋮tN1tN2⋯tNQ](12)\boldsymbol P=\begin{bmatrix} p_{11} & p_{12} & \cdots & p_{1Q} \\ p_{21} & p_{22} & \cdots & p_{2Q} \\ \vdots & \vdots & & \vdots \\ p_{M1} & p_{M2} & \cdots & p_{MQ}\\ \end{bmatrix},\boldsymbol T=\begin{bmatrix} t_{11} & t_{12} & \cdots & t_{1Q} \\ t_{21} & t_{22} & \cdots & t_{2Q} \\ \vdots & \vdots & & \vdots \\ t_{N1} & t_{N2} & \cdots & t_{NQ}\\ \end{bmatrix}\tag{12}P=⎣⎢⎢⎢⎡p11p21⋮pM1p12p22⋮pM2⋯⋯⋯p1Qp2Q⋮pMQ⎦⎥⎥⎥⎤,T=⎣⎢⎢⎢⎡t11t21⋮tN1t12t22⋮tN2⋯⋯⋯t1Qt2Q⋮tNQ⎦⎥⎥⎥⎤(12)其中，pijp_{ij}pij表示第jjj个训练样本的第iii个输入变量；tijt_{ij}tij表示第jjj个训练样本的第iii个输出变量；MMM为输入变量的维数；NNN为输出变量的维数；QQQ为训练集样本数。
则如上文所述，QQQ个隐含层神经元对应的径向基函数中心为C=P′(13)\boldsymbol {C=P'}\tag{13}C=P′(13)
<2>确定隐含层神经元阈值
为了简便起见，QQQ个隐含层神经元对应的阈值为b1=[b11,b12,⋯,b1Q]′(14)\boldsymbol {b}_1=[b_{11},b_{12},\cdots,b_{1Q}]'\tag{14}b1=[b11,b12,⋯,b1Q]′(14)其中，b11=b12=⋯=b1Q=0.8326spread,spreadb_{11}=b_{12}=\cdots=b_{1Q}=\frac{0.8326}{spread},spreadb11=b12=⋯=b1Q=spread0.8326,spread为径向基函数的扩展速度。
<3>确定隐含层与输出层间权值和阈值
当隐含层神经元的径向基函数中心及阈值确定后，隐含层神经元的输出便可以由式（15）计算：ai=exp(−∣∣C−pi∣∣2bi),i=1,2,⋯,Q(15)a_i=exp(-||\boldsymbol {C-p_i}||^2b_i),\quad i=1,2,\cdots,Q\tag{15}ai=exp(−∣∣C−pi∣∣2bi),i=1,2,⋯,Q(15)其中，pi=[pi1,pi2,⋯,piM]′\boldsymbol p_i=[p_{i1},p_{i2},\cdots,p_{iM}]'pi=[pi1,pi2,⋯,piM]′为第iii个训练样本向量。并记A=[a1,a2,⋯,aQ]\boldsymbol A=[a_1,a_2,\cdots,a_Q]A=[a1,a2,⋯,aQ]。
设隐含层与输出层间的连接权值W\boldsymbol WW为W=[w11w12⋯w1Qw21w22⋯w2Q⋮⋮⋮wN1wN2⋯wNQ](16)\boldsymbol W=\begin{bmatrix} w_{11} & w_{12} & \cdots & w_{1Q} \\ w_{21} & w_{22} & \cdots & w_{2Q} \\ \vdots & \vdots & & \vdots \\ w_{N1} & w_{N2} & \cdots & w_{NQ}\\ \end{bmatrix}\tag{16}W=⎣⎢⎢⎢⎡w11w21⋮wN1w12w22⋮wN2⋯⋯⋯w1Qw2Q⋮wNQ⎦⎥⎥⎥⎤(16)其中，wijw_{ij}wij表示第jjj个隐含层神经元与第iii个输出层神经元间的连接权值。
设NNN个输出层神经元的阈值b2\boldsymbol b_2b2为b2=[b21,b22,⋯,b2N]′(17)\boldsymbol b_2=[b_{21},b_{22},\cdots,b_{2N}]'\tag{17}b2=[b21,b22,⋯,b2N]′(17)
则由图2可得[Wb2]⋅[A;I]=T(18)[\boldsymbol{W\quad}\boldsymbol b_2]\boldsymbol\sdot[\boldsymbol{A;I}]=\boldsymbol T\tag{18}[Wb2]⋅[A;I]=T(18)其中，I=[1,1,⋯,1]1×Q\boldsymbol I=[1,1,\cdots,1]_{1×Q}I=[1,1,⋯,1]1×Q。
解线性方程组（18），可得隐含层神经元与输出层权值W\boldsymbol WW和阈值b2\boldsymbol b_2b2，即{Wb=T/[A;I]W=Wb(:,1:Q)b2=Wb(:,Q+1)(19)\begin{dcases} \boldsymbol{Wb}=\boldsymbol T/[\boldsymbol{A;I}]\\ \boldsymbol W=\boldsymbol{Wb}(:,1:Q) \\ \boldsymbol b_2=\boldsymbol{Wb}(:,Q+1)\tag{19} \end{dcases}⎩⎪⎨⎪⎧Wb=T/[A;I]W=Wb(:,1:Q)b2=Wb(:,Q+1)(19)
（3）RBF神经网络的MATLAB工具箱函数
函数newrbe用于创建一个精确的RBF神经网络，其调用格式如下：
new=newrbe(P,T,spread)\boldsymbol {new=newrbe(P,T,spread)}new=newrbe(P,T,spread)其中，PPP为网络输入量；TTT为网络目标量；spreadspreadspread为径向基函数的扩展速度（默认为1.0）；netnetnet为创建好的RBF网络。
二、案例背景
1、问题描述
辛烷值是汽油最重要的品质指标，传统的实验室检测方法存在样品用量大、测试周期长和费用高等问题，不适用于生产控制，特别是在线测试。近年来发展起来的近红外光谱分析方法（NIR），作为一种快速分析方法，已广泛应用于农业、制药、生物化工、石油产品等领域。其优越性是无损检测、低成本、无污染、能在线分析，更适合于生产和控制的需要。
针对采集得到的60组汽油样品，利用傅里叶近红外变换光谱仪对其进行扫描，扫描范围为0~450nm，扫描间隔为2nm，每个样品的光谱曲线共含401个波长点。样品的近红外光谱曲线如图3所示。同时，利用传统实验室检测方法测定辛烷值含量，现要求利用BP神经网络及RBF神经网络分别建立汽油样品近红外光谱及其辛烷值间的数学模型，并对模型的性能进行评价。
图3 60个样品的近红外光谱图
2、解题思路与步骤
依据问题描述中的要求，实现BP神经网络及RBF神经网络的模型建立及性能评价，大体上可以分为以下几个步骤，如图4所示。
图4 模型建立步骤
（1）产生训练集/测试集
为了保证建立的模型具有良好的泛化能力，要求训练集样本数量足够多，且具有良好的代表性。一般认为，训练集样本数量占总体样本数量的2/3~3/4为宜， 剩余的1/4~1/3作为预测集样本。同时，尽量使得训练集与预测集样本的分布规律近似相同。
（2）创建/训练BP神经网络
创建BP神经网络前需要确定网络的结构，即需要确定以下几个参数：输入变量个数、隐含层数及各层神经元个数、输出变量个数。如前文所述，只含有一个隐含层的三层BP神经网络可以逼近任意非线性函数，因此我们仅讨论单隐含层BP神经网络。从问题描述可知，输入变量个数为401，输出变量个数为1，隐含层神经元个数对网络性能的影响较大，具体讨论见本博客第四节。
（3）创建RBF神经网络
创建RBF神经网络时需要考虑spreadspreadspread的值对网络性能的影响，具体讨论参见第四节。
（4）仿真测试
模型建立后，将测试集的输入变量送入模型，模型的输出便是对应的预测结果。
（5）性能评价
通过计算测试集预测值与真实值的误差，可以对模型的泛化能力进行评价。在此基础上，可以进行进一步的研究和改善。
三、MATLAB程序实现
利用MATLAB神经网络工具箱提供的函数，可以方便地在MATLAB环境下实现上述步骤。
1、产生训练集/测试集
60个样品的光谱及辛烷值数据保存在spectra_data.mat中，该文件包含两个变量矩阵：NIR为60行401列的样品光谱数据，octane为60行1列的辛烷值数据。为不失一般性，这里采用随机法产生训练集和测试集，即随机产生50个样品作为训练集，剩余的10个样品作为测试集。具体程序如下：
%% 产生训练集/测试集load spectra_data.mat% 随机产生训练集和测试集temp = randperm(size(NIR,1));% 训练集——50个样本P_train = NIR(temp(1:50),:)';T_train = octane(temp(1:50),:)';% 测试集——10个样本P_test = NIR(temp(51:end),:)';T_test = octane(temp(51:end),:)';N = size(P_test,2);
说明：
（1）由于训练集/测试集产生的随机性，故每次运行的结果有可能不同；
（2）函数randperm(n)用于随机产生一个长度为n的正整数随机序列。
2、创建/训练BP神经网络及仿真测试
利用MATLAB神经网络自带工具箱的函数，可以方便地进行BP神经网络创建、训练及仿真测试。在训练之前，可以对相关的训练参数进行设置，也可以采取默认设置。具体程序如下：
%% BP神经网络创建、训练及仿真测试% 创建网络net = feedforwardnet(9);% 设置训练参数net.trainParam.epochs = 1000;net.trainParam.goal = 1e-3;net.trainParam.lr = 0.01;% 训练网络net = train(net,P_train,T_train);% 仿真测试T_sim_bp = net(P_test);
3、创建RBF神经网络及仿真测试
利用MATLAB神经网络自带工具箱的函数，可以方便地进行RBF神经网络创建及仿真测试。具体程序如下：
%% RBF神经网络创建及仿真测试% 创建网络net = newrbe(P_train,T_train,0.3);% 仿真测试T_sim_rbf = sim(net,P_test); % T_sim_bp = sim(net, P_test);旧版本
4、性能评价
BP神经网络及RBF神经网络仿真测试结束后，通过计算预测值与真实值的偏差情况，可以对网络的泛化能力进行评价。具体程序如下：
%% 性能评价% 相对误差errorerror_bp = abs(T_sim_bp - T_test)./T_test;error_rbf = abs(T_sim_rbf - T_test)./T_test;% 决定系数R^2R2_bp = (N * sum(T_sim_bp .* T_test) - sum(T_sim_bp) * sum(T_test))^2 / ((N * sum((T_sim_bp).^2) - (sum(T_sim_bp))^2) * (N * sum((T_test).^2) - (sum(T_test))^2)); R2_rbf = (N * sum(T_sim_rbf .* T_test) - sum(T_sim_rbf) * sum(T_test))^2 / ((N * sum((T_sim_rbf).^2) - (sum(T_sim_rbf))^2) * (N * sum((T_test).^2) - (sum(T_test))^2));% 结果对比result = [T_test' T_sim_bp' T_sim_rbf' error_bp' error_rbf']
本节选用的两个评价指标为相对误差和决定系数，其计算公式分别如下：
Ei=∣y^i−yi∣yi,i=1,2,⋯,n(20)E_i=\frac{|\hat y_i-y_i|}{y_i},\quad i=1,2,\cdots,n\tag{20}Ei=yi∣y^i−yi∣,i=1,2,⋯,n(20)R2=(l∑i=1ly^iyi−∑i=1ly^i∑i=1lyi)2(l∑i=1ly^i2−(∑i=1ly^i)2)(l∑i=1lyi−(∑i=1lyi)2)(21)R^2=\frac{\left(l\displaystyle\sum_{i=1}^l\hat y_iy_i-\sum_{i=1}^l\hat y_i\sum_{i=1}^ly_i\right)^2}{\left(l\displaystyle\sum_{i=1}^l\hat y_i^2-\left(\sum_{i=1}^l\hat y_i\right)^2\right)\left(l\displaystyle\sum_{i=1}^ly_i-\left(\displaystyle\sum_{i=1}^ly_i\right)^2\right)}\tag{21}R2=⎝⎛li=1∑ly^i2−(i=1∑ly^i)2⎠⎞⎝⎛li=1∑lyi−(i=1∑lyi)2⎠⎞(li=1∑ly^iyi−i=1∑ly^ii=1∑lyi)2(21)其中，y^i(i=1,2,⋯,n)\hat y_i(i=1,2,\cdots,n)y^i(i=1,2,⋯,n)为第iii个样品的预测值；yi(i=1,2,⋯,n)y_i(i=1,2,\cdots,n)yi(i=1,2,⋯,n)为第iii个样品的真实值；nnn为样品的数目。
说明：
（1）相对误差越小，表明模型的性能越好；
（2）决定系数在[0,1]内，越接近于1，表明模型的性能越好；反之，越接近于0，表明模型的性能越差。
5、绘图
为了更为直观地对结果进行观察和分析，这里以图形的形式将结果呈现出来。具体程序如下：
%% 绘图figureplot(1:N,T_test,'b:*',1:N,T_sim_bp,'r-o',1:N,T_sim_rbf,'k-.^')legend('真实值','BP预测值','RBF预测值')xlabel('预测样本')ylabel('辛烷值')string = {'测试集辛烷值含量预测结果对比(BP vs RBF)';['R^2=' num2str(R2_bp) '(BP)' ' R^2=' num2str(R2_rbf) '(RBF)']};title(string)
由于训练集和测试集都是随机产生的，每次运行结果都会不同。某次运行结果如图5所示。从图中可以清晰地看到，BP神经网络和RBF神经网络均能较好地实现辛烷值含量的预测。
图5 测试集辛烷值含量预测结果对比(BP vs RBF)
四、网络参数的影响与选择
1.隐含层神经元个数的选择
如前所述，隐含层神经元个数对BP神经网络的性能影响较大。若隐含层神经元的个数较少，则网络不能充分描述和输入变量之间的关系；相反，若隐含层神经元的个数较多，则会导致网络的学习时间变长，甚至会出现过拟合的问题。一般地，确定隐含层神经元个数的方法是在经验公式的基础上，对比隐含层不同神经元个数对模型性能的影响，从而进行选择。
图6描述了隐含层神经元个数对BP神经网络性能的影响，这里选取的评价指标为程序运行10次对应的决定系数的平均值。从图中可以看出，当隐含层神经元个数为7时，测试集的决定系数平均值最大，为0.9504。10次运行对应的具体结果如表1所列。
图6 隐含层神经元个数对BP神经网络性能的影响表1 10次运行对应的具体结果
2.spreadspreadspread值的选择
一般而言，spreadspreadspread值越大，函数的拟合就越平滑。然而，过大的spreadspreadspread将需要非常多的神经元以适应函数的快速变化；反之，若spreadspreadspread太小，则意味着需要许多的神经元来适应函数的缓慢变化，从而导致网络性能不好。
图7描述了不同spreadspreadspread值对RBF神经网络性能的影响，从图中可以看出，当spreadspreadspread取0.3时，网络的性能最好，此时对应的测试集决定系数为0.8608。
图7 spread值对RBF神经网络性能的影响
五、参考文献
[1] John H. Kalivas. Two data sets of near infrared spectra[J]. Chemometrics and Intelligent Laboratory Systems, 1997, 37(2): 255-259.
[2] Jeong H I, Lee H S, Jeon J H. Determination of Research Octane Number using NIR Spectral Data and Ridge Regression[J]. Bulletin, 2001, 22(1).
[3] 郁磊, 史峰, 王辉, 等. MATLAB智能算法30个案例分析(第2版)[M]. 北京: 北京航空航天大学出版社, .
补充
数据归一化是将数据映射到[0,1]或[-1,1]区间或其他的区间。由于（1）输入数据的单位不一样，有些数据的范围可能特别大，导致的结果是神经网络收敛慢、训练时间长。
（2）数据范围大的输入在模式分类中的作用可能会偏大，而数据范围小的输入作用就可能会偏小。
（3）神经网络输出层的激活函数的值域是有限制的，因此需要将网络训练的目标数据映射到激活函数的值域。例如神经网络的输出层若采用S形激活函数，由于S形函数的值域限制在(0,1)，也就是说神经网络的输出只能限制在(0,1)，所以训练数据的输出就要归一化到[0,1]区间。
（4）S形激活函数在(0,1)区间以外区域很平缓，区分度太小。例如S形函数f(X)在参数a=1时，f(100)与f(5)只相差0.0067。
以上四个原因，所以数据需要归一化。常用的归一化算法：
y=x−xminxmax−xmin,y∈(0,1)(22)y=\frac{x-x_{min}}{x_{max}-x_{min}},\quad y∈(0,1)\tag{22}y=xmax−xminx−xmin,y∈(0,1)(22)y=2(x−xmin)xmax−xmin+1,y∈(−1,1)(23)y=\frac{2(x-x_{min})}{x_{max}-x_{min}}+1,\quad y∈(-1,1)\tag{23}y=xmax−xmin2(x−xmin)+1,y∈(−1,1)(23)
由于本案例的数据不多，且范围有限，故没有进行数据归一化操作。如果对数据进行归一化操作的话，需要在神经网络创建之前加上归一化的程序，即：
%% 数据归一化[p_train, ps_input] = mapminmax(P_train, 0, 1);p_test = mapminmax('apply', P_test, ps_input);[t_train, ps_output] = mapminmax(T_train, 0, 1);
之后就使用归一化的训练集和测试集进行神经网络的创建、训练和学习；最后，一定要进行数据的反归一化，即：
%% 数据反归一化T_sim_bp = mapminmax('reverse', t_sim_bp, ps_output);
