一种基于GBDT和改进MFO的水质污染物检测方法与流程

一种基于gbdt和改进mfo的水质污染物检测方法
技术领域
1.本发明具体涉及一种基于gbdt和改进mfo的水质污染物检测方法，属于污染物检测技术领域技术领域。


背景技术：

2.水资源是人类生存的基础，水质安全直接影响着人民的生活质量。随着分子光谱技术的发展，基于激光诱导击穿光谱（laser-induced breakdown spectroscopy，libs）的水质检测方法具有检测速度快、无需试剂、对有机物灵敏度高等优点，在水质污染物原位在线分析等方面已展现出明显优势。
3.libs在测量过程中，由于黑体辐射、负荷辐射等产生的连续背景辐射干扰，会使光谱基线发生偏移，影响数据分析的准确性。目前常用的基线校正方法有分段拟合法、迭代多项式拟合法、小波变换算法等。但是，现有的基线校正方法存在计算复杂、不稳定、无法完成特定元素的基底去除等问题。
4.将libs技术与机器学习方法结合可实现对物质的快速分析。梯度提升树（gradient boosting decision tree，gbdt）算法能对光谱峰值与污染元素含量百分比之间进行精准预测，但其准确度受学习率、树的最大深度和弱回归树的个数影响很大，若参数选择不当，会降低分析结果的准确度。故需要使用优化算法对gbdt进行参数寻优。mfo算法中的飞蛾种群位置是随机生成，易造成飞蛾位置分布不均匀，影响其收敛速度和求解精度。


技术实现要素：

5.本发明的目的在于针对现有技术不足，提出一种基于gbdt和改进mfo的水质污染物检测方法，采用混沌初始化的方法对原始飞蛾火焰优化（moth-flame optimization algorithm，mfo）进行改进，能有效避免算法陷入局部最优解。然后，将改进后的混沌初始化飞蛾火焰优化（moth-flame optimization algorithm based on chaotic initialization，cmfo）算法与gbdt结合，能有效提高libs光谱数据的分析精度，为水质污染物的在线检测提供一种新方法。
6.为解决上述技术问题，本发明采用的技术方案为：一种基于gbdt和改进mfo的水质污染物检测方法，包括如下步骤：步骤1：采集污染水源样品，对污染水源样品进行libs光谱检测，形成光谱图；步骤2：基于待检测污染物波长，获取待检测污染物在该波长处的光谱强度和实际浓度；步骤3：对待检测污染物的光谱强度和实际浓度进行归一化处理，构建回归梯度提升树模型gbdt；步骤4：采用混沌初始化映射对mfo算法进行改进，得到cmfo优化算法；步骤5：基于改进后的cmfo算法，优化gbdt模型参数，对污染物含量进行预测。
7.进一步的，步骤2用波长i处的光谱面积替代光谱强度，以更全面的运用光谱信息，
提高精度，波长i处的光谱面积的计算公式为：其中，s
raw
为波长i在波长范围内的面积，m为波长范围的第一个像素点， xm为m像素点处的光谱强度，n为波长范围的最后一个像素点， xn为n像素点处的光谱强度，xk为k像素处的光谱强度，k为[m+1，n]范围内的整数。
[0008]
进一步的，步骤2对待测污染物波长i处的光谱面积进行基底噪音处理，提高模型精度。
[0009]
优选的，基底噪音处理的具体操作步骤为：步骤a：基于污染水源的libs光谱图，将待检测污染物波长i对应像素点前后v个波长分成两个光谱段a、b，其中v是正整数；步骤b：分别计算两个光谱段a和b的面积；步骤c：通过两个光谱段a和b的平均面积与待检测污染物波长i像素长度的乘积，计算基底噪音面积；步骤d：波长i处的光谱面积去除基底噪音面积，即为修正后波长i处的光谱面积。
[0010]
进一步的，基于下述公式修正波长i处的光谱面积：上式中，x
i-v
为波长为i-v处的光谱强度，x
m-1
为波长m-1处的光谱强度，x
l
为波长为l处的光谱强度，x
n+1
为波长n+1处的光谱强度，x
i＋v
为波长i＋v处的光谱强度，s
background
为基底噪音面积，s
correction_i
为修正后波长i处的光谱面积。
[0011]
进一步的，步骤4采用logistics映射进行混沌初始化，解决初始化不均匀问题。
[0012]
优选的，在mfo中飞蛾的位置定义为待解决的问题，即污染物预测值，火焰被定义为飞蛾所搜索到的局部最优位置，通过迭代更新飞蛾的位置求解最优解，采用混沌初始化映射对mfo算法进行改进，混沌初始化公式为：其中，
µ
∈（0,4]为控制变量，y
p
表示第p只飞蛾在混沌空间中的位置，其取值范围为[0,1]；y
p+1
表示第p+1只飞蛾在混沌空间中的位置；在混沌空间中随机生成第一只飞蛾的位置y
p1
，对y
p1
进行n0-1次迭代，得到其余n0-1只飞蛾位置，再将所有飞蛾位置映射到搜索空间，作为迭代初值：
其中，ub、lb为搜索空间的搜索范围的上、下限，m
p
表示第p只飞蛾的位置，p=1,2,
⋯
,n0，n0是飞蛾最大迭代次数。
[0013]
进一步的，第p只飞蛾m
p
会向距离最近火焰f(n1
×
d)逼近，其中n1表示飞蛾种群数，d表示飞蛾种群维度，f是飞蛾矩阵，逼近的路线为螺旋曲线。
[0014]
优选的，逼近函数为：其中，fj为f(n1
×
d)中距离m
p
最近的火焰；d
p
=|f
j-m
p |为飞蛾距离火焰的距离；为对数螺旋曲线相关的常量；随机数t定义了飞蛾下一个位置离火焰的远近程度，介于[-1，1]之间；为对数螺旋曲线表达式，与d
p
相乘后是一个常量，与fj相加后是向量，表示更新后的飞蛾位置，即将飞蛾位置m
p
代入逼近函数后，自动更新为m
p
'，完成新位置的更新。
[0015]
进一步的，为了算法寻优能力更强，在迭代过程中火焰数量需逐渐减少，火焰减少过程如下：其中，flameno为每轮迭代减少后的火焰数量；round为四舍五入函数；为当前迭代次数，为迭代次数的最大值，n为火焰数量的最大值。
[0016]
本发明与现有技术相比，具有如下有益效果：本发明的目的在于针对现有技术不足，提出一种基于gbdt和改进mfo的水质污染物检测方法，采用混沌初始化的方法对原始飞蛾火焰优化（moth-flame optimization algorithm，mfo）进行改进，能有效避免算法陷入局部最优解。然后，将改进后的混沌初始化飞蛾火焰优化（moth-flame optimization algorithm based on chaotic initialization，cmfo）算法与gbdt结合，能有效提高libs光谱数据的分析精度，为水质污染物的在线检测提供一种新方法。
[0017]
1、本发明提出用波长i处的光谱面积替代光谱强度，以更全面的运用光谱信息，提高精度。
[0018]
2、本发明鉴于水质污染物分析是针对特定污染元素进行测，对待测污染物波长i处的光谱面积进行基底噪音处理，相比于传统导数法、迭代多项式拟合法、分段拟合算法、小波变换算法等整体基底去除方法，该方法计算简单、稳定，能很好的完成特定元素的基底去除，提高模型精度。
[0019]
3、本发明设计了修正波长i处的光谱面积公式，考虑前后波段的情况，可以更准确的计算光谱强度，提高模型精度。
[0020]
4、本发明采用混沌初始化的方法对原始mfo进行改进，能有效避免算法陷入局部最优解。将改进后的混沌初始化飞蛾火焰优化算法与gbdt结合，能有效提高libs光谱数据的分析精度。
附图说明
[0021]
图1为本发明的步骤框架图；图2为本发明的流程图；图3为本发明实施例中fe在238.204nm处的谱线图。
具体实施方式
[0022]
下面将结合本发明中的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例，都属于本发明保护的范围。
[0023]
参见图1-图3，本实施例提供一种基于gbdt和改进mfo的水质污染物检测方法， 包括以下步骤：步骤s1：采集污染水源样品，对污染水源样品进行libs光谱检测，形成光谱图；步骤s2：基于待检测污染物波长，获取待检测污染物在该波长处的光谱强度和实际浓度；采集污染水源（实验室样品，已知各污染物浓度）的libs光谱数据进行处理，提取预检测（技术人员可根据实际检测需求，在nist数据库找到待测量污染物的波长）污染物的光谱波段，构成满足分析需求的libs原始数据集｛(x
i，
yi)｝。其中，i为根据实际需求选取的待测污染物波长（通常每次选1个波长）；xi为待测污染物在波长i处的光谱强度，yi为待测污染物实际浓度。待测量污染物包括硫元素、铁元素、铝元素、镍元素、钡元素等等。选取特定波长进行分析，而不采用全光谱进行分析，即特征选择过程，可以减少后期模型训练的复杂度。
[0024]
具体包括以下步骤：s201：libs原始数据集｛(x
i，
yi)｝由波长峰值（即待测污染物在波长i处的光谱强度）和对应元素含量百分比（即待测污染物实际浓度）组成，其中不同波长对应不同元素，峰值代表该元素含量的高低。由于libs的单次检测中，激光的能量，样品表面不平整度，仪器其他方面的不稳定，都会随时间存在一定的波动和不确定性而导致光谱数据存在偏差，故取5组数据进行平均，去除异常数据影响。由此获得平均后的的数据集｛(x
i，
yi)｝；s202：根据nist数据库提供的元素-波长对应关系选取污染物波长：硫元素（545.383nm、563.997nm）、铁元素（238.204nm、238.863nm、239.563nm）、铝元素（308.215nm、309.271nm）、镍元素（221.648nm）、钡元素(495.709nm)，该选择方案可根据污水实际情况做出调整，元素和波长的对应关系是不唯一时，通常选取峰值最大的那个波长，或根据光谱图像，人工判断出受干扰最小的波长。libs仪器通常由美国制造，所以本发明通常采用美国国家标准与技术研究院(nist)数据库。
[0025]
s203：去除libs光谱数据的基底噪声：取所选波长i对应像素点范围前后v波长的两个光谱段a、b，通常v为[20,30]间的整数。其中i为需要去除基底噪声的元素波长，即步骤s102中选择的波长。如附图3所示，图3为fe在238.204nm处的谱线，其波长范围对应libs光谱仪5798-5808像素单元。求出面积sa、sb，其中sa是光谱段a的面积，sb是光谱段b的面积。将平均面积（sa+sb）/2与i波长像素长度的
乘积定义为基底噪音面积。具体计算过程为：其中，s
raw
为波长i在波长范围内的面积；m为波长范围的第一个像素点，图3中对应像素点5798；xm为m像素点处的光谱强度；n为波长范围的最后一个像素点，图3中对应像素点5808；xn为n像素点处的光谱强度；xk为k像素处的光谱强度，k为[m+1，n]范围内的整数；上式中，x
i-v
为波长为i-v处的光谱强度，x
m-1
为波长m-1处的光谱强度，x
l
为波长为l处的光谱强度，x
n+1
为波长n+1处的光谱强度，x
i＋v
为波长i＋v处的光谱强度，s
background
为基底噪音面积，s
correction_i
为修正后波长i处的光谱面积。
[0026]
s204：使用波长对应面积代替光谱强度，即将所选波长i的光谱强度替换为修正后的光谱面积s
correction_i
。对多组libs数据重复步骤s101-s104，得到满足训练模型所需求的数据集t'，其中，t'中包含多组光谱面积和污染物实际浓度数据组合{s
correction_i_r
,y
i_r
},r=1,2,3,
⋯
，其中{s
correction_i_r
,y
i_r
}为t'的第r组数据；s
correction_i_r
为第r组修正后波长i处的光谱面积；y
i_r
为第r组波长i所对应元素的实际含量。
[0027]
s205：对步骤s104的数据t'进行归一化处理，得到符合模型训练的数据集t
″
={xr,yr},r=1,2,3,
⋯
。其中xr为归一化后的光谱面积，yr为归一化后的元素含量实际值。
[0028]
步骤s3：对待检测污染物的光谱强度和实际浓度进行归一化处理，构建回归梯度提升树模型gbdt；步骤s3包括如下步骤：s301：定义回归梯度提升树模型gbdt。gbdt需经过多次迭代构建多个决策树来组成最终的模型。每次迭代过程中，每一颗决策树沿着梯度下降最快的方向降低残差。
[0029]
将步骤s1处理后的数据t
″
分成训练集和测试集，gbdt模型使用的训练样本集为t
″
中选取出的训练集，将多组xr作为模型输入，模型将多组yr划分为a个区域r1，r2，
…
，ra，即叶子节点个数，并在每个区域上确定输出值ca，ca为gbdt模型的输出。提升树可表示为：为gbdt模型的输出。提升树可表示为：其中，xr为该棵提升树的输入，θ表示输入的区域划分和各区域上的取值；i(xr∈
ra)是一个指示函数，当括号里的条件成立时取1，否则取0。
[0030]
s302：初始化。损失函数为：确定初始提升树，求得常数c使损失函数达到极小值：其中，损失函数中的是模型预测值，是样本真实值，即归一化后的元素含量实际值；为∑l(yr,c)取最小值时c对应的值；f0(xr)为第0轮的gbdt模型；因为第一次迭代会有模型f1(xr)，故须先定义f0(xr)。
[0031]
代入损失函数得：对∑(y
r-c)2中的c求导，令其导数为0，得到关系式：其中m为输入样本的个数。即:s303：gbdt模型迭代。每轮迭代都在上一轮迭代完成的模型基础上再增加一颗回归提升树，用于减小上一轮迭代的残差。对b=1,2,
…
,b（b表示迭代次数），第b步所得模型为：其中，fb(x)是第b轮迭代的模型；f
b-1
(x)是第b-1轮迭代得到得模型，为已知表达式；t(x;θb)为该轮迭代中需求解的模型，即第b颗树，其参数θb通过经验风险极小化来求解；其中，等号右侧为∑l[yi,f
b-1
(x)+t(x;θb)]取最小值时的取值。
[0032]
采用步骤s302定义的损失函数，其损失为：其中，r
bj
为当前模型拟合数据的残差。在第b轮迭代中需要求解的模型t(x;θb)= r
bj
时，
损失函数最小。对m个样本t
″
，迭代轮数b=1,2,
…
,b，使用r
bj
拟合t(x;θb)。拟合样本为{(xr,r
bi
)}，得到第b棵树的叶子节点区域r
ba
,a=1,2,
⋯
,a，a的取值由算法自动生成。每个区域上的取值计算如下：采用步骤s302定义的损失函数时，对c
ba
求导，令其导数等于0即可求得c
ba
。
[0033]
s304：更新模型。
[0034]
此处x=xr,若第b轮时最后一次迭代，则最终的回归提升树模型为：式中，r
ba
是第b轮迭代中的第a个叶子节点；c
ba
为r
ba
叶子节点中的取值；即任意输入libs光谱强度转换后的面积值xr,若xr被划分到叶子节点r
ba
中，则模型给出预测值，预测xr光谱对应的污染物含量为c
ba
%。
[0035]
步骤s4：采用混沌初始化映射对mfo算法进行改进，得到cmfo优化算法；步骤s4包括如下步骤：s401：为了解决初始化不均匀带来的问题，采用logistics映射进行混沌初始化：其中，
µ
∈（0,4]为控制变量，y
p
表示第p只飞蛾在混沌空间中的位置，其取值范围为[0,1]；y
p+1
表示第p+1只飞蛾在混沌空间中的位置；s402：在混沌空间中随机生成第一只飞蛾的位置y
p1
，对y
p1
进行n0-1次迭代，得到其余n0-1只飞蛾位置，再将所有飞蛾位置映射到搜索空间，作为迭代初值：其中，ub、lb为搜索空间的搜索范围的上、下限，m
p
,p=1,2,
⋯
,n0为第p只飞蛾的位置；s403：在cmfo中，飞蛾的位置定义为待解决的问题，火焰被定义为飞蛾所搜索到的局部最优位置，通过迭代更新飞蛾的位置求解最优解，飞蛾种群的空间位置用矩阵表示，火焰的空间位置用矩阵f(n1
×
d)表示，飞蛾和火焰的适应度值分别为和of(n
×
1)。其中，n为飞蛾种群数量，即为n只飞蛾总数，d为搜索空间的维度，即待寻优参数的个数。
[0036]
cmfo算法流程表述如下：函数i’代表初始飞蛾种群以及相应的适应度函数值。其系统模型如下所示：
式中，代表飞蛾空间位置矩阵；代表存储每只飞蛾适应度的矩阵；φ抽象了函数i’的运算过程，即函数i’实现了初始化飞蛾种群和计算每只初始位置飞蛾的适应度。
[0037]
函数p决定飞蛾在解空间中的运动轨迹，飞蛾空间位置矩阵输入p后更新为m’，通过不断迭代，进而更新自己的位置。
[0038]
函数t
″
判断算法是否达到预设的停止条件，达到则为真true，否则为假false。
[0039]
综上所述，函数i’初始化飞蛾位置，具体初始化方法如s401-s402所述，函数p不断迭代，直至函数t
″
返回为真，即达到预期效果。
[0040]
在飞蛾种群数量、位置和适应度函数确定后，第p只飞蛾m
p
会向距离最近火焰f(n1
×
d)逼近，其中n1表示飞蛾种群数，d表示飞蛾种群维度，f是飞蛾矩阵，逼近的路线为螺旋曲线。cmfo的对数螺旋函数表示为：其中，fj为f(n1
×
d)中距离mp最近的火焰；dp=|fj-mp |为飞蛾距离火焰的距离；为对数螺旋曲线相关的常量，本方案中取1；随机数t定义了飞蛾下一个位置离火焰的远近程度，介于[-1，1]之间；为对数螺旋曲线表达式，与d
p
相乘后是一个常量，与fj相加后是向量，表示更新后的飞蛾位置。即将飞蛾位置m
p
代入逼近函数后，自动更新为m
p
'，完成新位置的更新。
[0041]
为了算法寻优能力更强，在迭代过程中火焰数量需逐渐减少，火焰减少过程如下：其中，flameno为每轮迭代减少后的火焰数量；round为四舍五入函数；为当前迭代次数，为迭代次数的最大值，n为火焰数量的最大值。
[0042]
步骤s5：引入改进后的cmfo算法，优化gbdt模型参数，对污染物含量进行预测；步骤s5包括如下步骤：s501：根据需求，gbdt需要寻优的参数有学习率σ、树的最大深度g和弱回归树的个数γ，故步骤s403中的搜索空间的维度d=3；s502：步骤s501中飞蛾和火焰数目n=60，最大迭代次数l
max
=100，飞蛾搜索范围的上限ub=10，搜索范围下限lb=10；该处取值可根据实际需求进行调整；s503：按s401和s402的公式初始化飞蛾种群，在搜索空间内随机生成飞蛾位置:其中为飞蛾矩阵的第p行向量；每只飞蛾的位置都是存在于一个矩阵中，例
如矩阵为n
×
d列，n为飞蛾的只数（矩阵行数），d为维度，即待寻优参数的个数（矩阵列数）；第一只飞蛾的位置存在矩阵的第一行，故第p只飞蛾的位置存在第p行；s504：将m
p1
、m
p2
、m
p3
输入gbdt模型进行污染物含量预测，根据预测结果，可以计算出相应的适应度函数值，本发明采用平均绝对百分比误差（mean absolute percentage error，mape）为适应度函数，计算公式为：其中，m为污染物预测样本数；xi对应上文s
correction
；f(x)于s304中最终的回归提升树的公式定义；yi、f(xi)分别为污染物的实际含量与预测值；s505：将飞蛾位置按适应度值从小到大排序，赋给火焰；s506：飞蛾按照火焰减少过程公式更新其位置，其中=1，该处取值可根据实际需求进行调整；s507：将更新后的飞蛾位置带入gbdt模型，重复步骤s504、s505；s508：根据火焰减少过程公式减少火焰数量，迭代次数；s509：判断是否达到最大迭代次数（此处最大迭代次数设为100，当时，则达到最大迭代次数），若达到则输出最优火焰位置即全局最优解（输出值为三个常数，分别为学习率σ、树的最大深度g和弱回归树的个数γ取值），否则转到步骤s506。
[0043]
使用训练好的cmfo-gbdt模型对污染物的含量进行预测，根据污染物含量指标来判断水质是否满足要求。
[0044]
本发明提供一种基于gbdt和改进mfo的水质污染物检测方法，该法首先去除基底，相比于传统导数法、迭代多项式拟合法、分段拟合算法、小波变换算法等基底去除方法，该方法计算简单、稳定，能很好的完成特定元素的基底去除。本发明采用混沌初始化的方法对原始mfo进行改进，能有效避免算法陷入局部最优解。将改进后的混沌初始化飞蛾火焰优化算法与gbdt结合，能有效提高libs光谱数据的分析精度。
[0045]
在此处所提供的说明书中，说明了大量具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下被实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。