CN105844579A - 一种基于近红外光谱的原烟叶组配方维护方法 - Google Patents

Abstract

本发明公开了一种基于近红外光谱的原烟叶组配方维护方法，包括：步骤1，依据去年和今年烟叶品种的近红外光谱和化学值的一致性程度，将去年的烟叶品种划分为共性等级和个性等级；步骤2，针对共性等级，采用与去年相同的比例；针对个性等级，利用今年采集的烟叶品种组合成混合样，使混合样的近红外光谱以及化学值分别与对应的个性等级相一致，混合样的比例采用去年原烟叶组配方中对应的个性等级的比例，得到今年的初步原烟叶组配方；步骤3，若今年和去年制样的化学值和近红外光谱差异符合预期，得到今年的原烟叶组配方。本发明以近红外光谱为依据，以化学成分作为判断基础，当收购原烟的指标发生变化时，能够快速确定质量稳定的原烟叶组配方。

Description

一种基于近红外光谱的原烟叶组配方维护方法

技术领域

本发明涉及原烟技术领域，具体涉及一种基于近红外光谱的原烟叶组配方维护方法。

背景技术

烟叶近红外光谱中包含大量的物质信息，近红外光谱满足两个先决条件：

(1)相似的光谱必有相似的物理性质与化学性质；

(2)近红外光谱在一定条件下满足加和规律。

原烟作为一种农产品，历年的化学成分、以及烟叶外观质量波动较大，以人为主的原烟叶组配方设计，在实际的具体操作过程中通常会遇到很多问题，比如：(1)如何维护年度间原烟叶组配方质量的一致性？(2)烟叶收购的产地与烟叶等级发生变化时，如何快速的确定新的原烟叶组配方单？

随着打叶复烤均质化的逐渐推进，对原烟质量的要求也越来越高，不仅要求不同批次比较均质，也要求年度间的原烟叶组配方比较稳定，传统依靠人工类维护叶组配方的方式已经不能够满足实际的生产过程需要。

据公开资料显示，现在将近红外光谱运用在原烟叶组配方维护的工作主要是运用近红外投影后的MSPC(多变量过程控制)对烟叶配方的稳定性进行分析，还没有将近红外光谱运用于原烟叶组配方的维护。

发明内容

本发明提供了一种基于近红外光谱的原烟叶组配方维护方法，以近红外光谱为依据，以化学成分作为判断基础，当收购原烟的产地、等级、内在化学成分与感官发生变化时，能够快速确定质量保持稳定的新的原烟叶组配方。

一种基于近红外光谱的原烟叶组配方维护方法，包括：

步骤1，针对去年原烟叶组配方中的烟叶品种以及今年采集的烟叶品种，分别进行近红外光谱采集以及化学值测定，依据近红外光谱和化学值的一致性程度，将去年的烟叶品种划分为共性等级和个性等级；共性等级和个性等级的划分标准为：近红外光谱PCA投影第一得分值及化学值幅值的差异均在10％以内的作为共性等级，其余为个性等级。

步骤2，针对共性等级，在今年的原烟叶组配方中采用与去年相同的比例；

针对个性等级，利用今年采集的烟叶品种组合成混合样，使混合样的近红外光谱以及化学值分别与对应的个性等级相一致，在今年的原烟叶组配方中，混合样的比例采用去年原烟叶组配方中对应的个性等级的比例，得到今年的初步原烟叶组配方；

步骤3，分别依据去年的原烟叶组配方和今年的初步原烟叶组配方进行制样，并对制样分别进行化学值测定和近红外光谱采集，若化学值和近红外光谱差异符合预期，得到今年的原烟叶组配方。

本发明中以近红外光谱为基础，尊重传统的叶组配方维护工作的精华，首先，在叶组配方的维护工作中，没有直接针对已有叶组配方进行近红外光谱的组合，而是首先研究去年原烟叶组配方表中的具体烟叶品种的差异，找出跨年度化学值与近红外光谱均差异不大的原料作为跨年度原烟叶组配方重新设计的共性等级，对于这部分共性等级部分的烟叶比例不做改动，有效保证了跨年度原烟叶组配方设计的整体感官风格。

对于跨年度不能通用的样本等级，也称为跨年度原烟叶组配方设计的个性等级，对这部分等级，利用双权重的思想去替代。具体而言，首先找到需要替代的目标个性化等级在近红外光谱上与之最相似的一系列的光谱组，分别赋予一定的相似度权重，在此基础上，按照光谱加和为目标等级光谱的原则去计算具体的比例。

从具体的计算方式上来看，这种运用近红外光谱对原烟叶组配方维护的方式包含了传统人工对原烟叶组配方的维护方式，但又不仅局限于此，进一步采用了近红外光谱的加和规律，保证原烟叶组配方的质量与效率。

近红外光谱信息中包含烟草样本的物理信息、化学信息、产地信息以及烟叶部位信息等多种信息，因此，近红外光谱的差异在一定程度上能够反映出烟草样本的品质特性。

为了保证原烟叶组配方的可靠性，基于近红外光谱拟合得到的原烟叶组配方需要进一步通过验证，方能作为今年的原烟叶组配方，优选地，还包括步骤4，对制样进行感官评析验证，感官评析符合预期后，确定今年的原烟叶组配方。

作为优选，步骤1中，若去年原烟叶组配方中的烟叶品种和今年采集的烟叶品种的化学值和近红外光谱差异均不大，则判定为共有等级；去年原烟叶组配方中的其余烟叶品种判定为个性等级。在进行近红外光谱以及化学值的对比时，可以采用现有技术。

作为优选，步骤1中划分共性等级的步骤包括：

选取去年和今年共有的烟叶品种，对共有的烟叶品种进行近红外光谱和化学值的比对，若近红外光谱和化学值的差异均不大，则判定为共有等级。

为了保证化学值比对的结果能够充分反映原烟叶组配方的相似性，优选地，步骤1采集的化学值至少包括烟碱含量、总糖含量、还原糖含量、总氮含量、钾含量、氯含量。

若若化学值和近红外光谱差异不符合预期，通过人工进行微调，得到今年的原烟叶组配方。

本发明中，结合上年度的原烟叶组配方，各等级原烟的化学值与近红外光谱，与本年度的原烟各等级的基本化学值与近红外光谱的情况，以化学值作为判断原烟叶组配方维护一致性的前提，以近红外光谱作为具体的等级合并与烟叶等级比例的计算依据，再以最后的化学值，近红外光谱，感官评吸与外观数据作为验证本年度叶组配方与往年度叶组配方的一致性的可靠性；在原料化学值没有产生很大变化时，运用近红外光谱快速的维护叶组配方，在原料的化学值发生重大变化时，在保证烟叶可用性的前提下，尽量缩小本年度与往年度叶组配方质量的差异性。

附图说明

图1为2014年取样的烟叶样本的近红外光谱；

图2为2014年近红外光谱的投影向量计算结果；横坐标为光谱变量序号，纵坐标为投影向量系数值。

图3为2014年近红外光谱投影结果；横坐标为共性样本序号，纵坐标为光谱投影值。

图4为2014年的个性烟叶样本(近红外光谱PCA投影第一得分值及化学值幅值的差异均超过10％)的近红外光谱；横坐标为个性样本序号，纵坐标为光谱投影值。

图5为利用2015年烟叶样本拟合图4的近红外光谱得到的拟合光谱图；横坐标为个性样本序号，纵坐标为投影拟合差值。

图6为2015年的制样结果与2014年制样结果的对比图。横坐标为光谱变量序号，纵坐标为吸光值。

具体实施方式

下面结合附图，对本发明基于近红外光谱的原烟叶组配方维护方法做详细描述。

如图1所示，基于近红外光谱的原烟叶组配方维护方法，包括：

(1)获取2014年的原烟叶组配方表记为W1，W1如表1所示，表1中第一列为烟叶等级的编号，第二列为等级烟叶的数量，第三列为不同等级烟叶的比重，各种烟叶的等级比例记作P。

表1

按照2014年的原烟叶组配方进行各个等级的取样，并对取样得到的烟叶样本进行去梗，磨粉过40目筛后，进行实验室近红外光谱采集，获取实验室近红外光谱S1，如图1所示，；对该部分烟草样本同时进行流动检测，记录常规化学成分含量C1(即化学值)，如表2所示。

表2

(2)2015年的烟叶收购上来后，对烟叶按照分选批次进行分选，对每个分选等级都进行抽样，并对取样得到的烟叶样本进行去梗，磨粉过40目筛后，进行实验室近红外光谱采集，获取实验室近红外光谱S2；对该部分烟草样本同时进行流动检测，记录常规化学成分含量C2，如表3所示。

表3

寻找2014年原烟叶组配方表W1中具体的取样等级包含在2015年的烟叶等级中的烟叶等级记作D1；没有包含在2015年的烟叶收购中的等级记作D2。

2014年等级为D1的烟叶样本的化学成分含量为D1C，近红外光谱为D1S；将2015年的烟叶样本中与D1等级相对应的烟叶样本的等级的记作D1N，采集烟叶样本D1N的近红外光谱D1NS和化学成分含量D1NC，对比化学值D1C与化学值D1NC，近红外光谱D1S与D1NS的差异。

如果化学值D1C与化学值D1NC，近红外光谱D1S与D1NS的差异均不大，则对应的2014年的这部分烟叶样本记为D1G，D1G即为2014年烟叶中的共性等级，共性样本在2015年的原烟叶组配方比例中，按照2014年的原烟叶组配方比例进行设计，此比例记作PN1。

如果化学值D1C与化学值D1NC较大，或者近红外光谱D1S与D1NS的差异较大，则将差异较大的样本记为D2G，D2G为2014年烟叶中的个性样本。

以2014年和2015年为例，对2014年与2015年的近红外光谱进行投影，寻找化学值与近红外光谱差异均不大的等级，具体的计算方式如下：

a、将2014年与2015年的近红外光谱合并在一起形成光谱矩阵x，按照如下的公式进行投影：

$x_{ij}^{,}=\frac{x_{ij-}\min \left(x_j\right)}{\max \left(x_j\right)-\min \left(x_j\right)}$

其中：i为合并后光谱矩阵x的第i个样本，j为光谱矩阵的第j列；x_ij为光谱矩阵x的元素。

投影实质上是从不同的角度去观察数据，寻找最能充分挖掘数据特征的最优投影方向。在单位超球面中随机抽取若干个初始投影方向a＝(a₁,a₂,…,a_m)，计算各初始投影方向的投影指标大小，根据指标选大的原则，确定最大指标对应的解为最优投影方向。

若(a₁,a₂,…,a_m)为m维单位向量，则第i个烟叶样本在一维线性空间的第i个烟叶样本的烟叶投影特征值z_i的表达如下：

$z_i=\underset{j=1}{\overset{m}{\Σ}}{a}_j{x}_{ij}^{\′}$

其中，z_i为第i个烟叶样本的烟叶光谱投影特征值，a_j为第j个初始投影向量，x′_ij为归一化后的红外光谱数据。

b、确定目标函数。综合投影指标值时，要求投影特征值zi的散布特征应为：局部投影点尽可能密集，最好凝聚成若干个点团，而在整体上投影点团之间尽可能散开。在本实施例中，目标函数定位为初始投影向量的类内距离与投影向量的类内密度的乘积，即

Q(a)＝s(a)·d(a)

其中，a为初始投影向量，Q(a)为目标函数，s(a)为初始投影向量的类内距离，d(a)为初始投影向量的类内密度。

在本实施例中，初始投影向量的类内距离s(a)采用第i个烟叶样本的烟叶投影特征值zi的方差计算，计算公式如下：

$s\left(a\right)={\[\underset{i=1}{\overset{n}{\Σ}}{(z_i-{\stackrel{\‾}{z}}_a)}^2/n\]}^{1/2}$

其中，z_i为第i个烟叶样本的烟叶光谱投影特征值，为序列{z(i)|i＝1,2,…,n}的均值，n为烟叶样本的个数。

初始投影向量的类内密度d(a)利用任意两个烟叶样本的烟叶投影特征值间的距离r_ik计算，任意两个烟叶样本的烟叶光谱投影特征值间的距离r_ik的计算公式如下：

r_ik＝|z_i-z_k|(i,k＝1,2,…,n)

其中，z_i为第i个烟叶样本的烟叶光谱投影特征值，z_k为第k个烟叶样本的烟叶光谱投影特征值。

初始投影向量的类内密度d(a)的计算公式如下：

$d\left(a\right)=\underset{i=1}{\overset{n}{\Σ}}\underset{k=1}{\overset{n}{\Σ}}(R-r_{ik})f(R-r_{ik})$

其中，f(t)为一阶单位阶跃函数，t≥0时，其值为1；t＜0时，其值为0，R为估计局部散点密度的窗宽参数，i，k为烟叶样本的序号，n为烟叶样本的个数。

$f(R-r_{ik})=\left\{\begin{array}{cc}1& R\&GreaterEqual;r_{ik}\\ 0& R<r_{ik}\end{array}\right.$

其中，R为估计局部散点密度的窗宽参数，按宽度内至少包括一个散点的原则选定，其取值与样本数据结构有关，基本确定它的合理取值范围为r_max＜R≤2m，其中，r_max＝max(r_ik)(i，k＝1，2，…，n)。类内密度d(a)愈大，分类愈显著。

c、优化线性投影的方向。在本实施例中，当目标函数Q(a)取最大值时所对应的投影方向即为所需最优投影方向。因此，寻找最优投影方向的问题可转化为下列优化问题：

$\left\{\begin{array}{c}\max Q\left(a\right)=s\left(a\right)\&CenterDot;d\left(a\right)\\ ||a||=\underset{j=1}{\overset{m}{\&Sigma;}}{a}_j^2=1\end{array}\right.$

在本实施例中，这是以a_j为优化变量的复杂非线性优化问题，可采用遗传算法等优化方法求解。

共性样本可替代的原则需要投影光谱的差异比较小以及化学值的差异比较小。2014年光谱的投影向量计算结果如图2所示，图2中横坐标为光谱变量序号，纵坐标为投影向量系数值。光谱的投影结果如图3所示，图3中横坐标为共性样本序号，纵坐标为光谱投影值。

投影光谱差异的计算公式如下：

r_ik＝|z_i-z_k|(i，k＝1，2，…，n)。

化学值差异的计算公式如下：

设不同等级的化学值为CH，用CHP_ik＝|CH_i-CH_k|(i，k＝1，2，…，n)。2014年和2015年不同等级的烟叶样本的化学值差异如表4所示。

表4

序号	年份	编号2	尼古丁	共性样本序号	年份	尼古丁	绝对误差
								1	2015	2	2.50	2	2014	2.86	0.36
2	2015	2	2.57	6	2014	2.67	0.1
								3	2015	2	2.37	8	2014	2.4	0.03
4	2015	5	2.52	9	2014	2.34	0.18
								5	2015	8	2.51	10	2014	2.37	0.14
7	2015	10	2.70	12	2014	2.63	0.07
								8	2015	12	3.85	14	2014	3.45	0.4
9	2015	18	2.60	20	2014	2.63	0.03
								10	2015	19	2.83	21	2014	3.09	0.26
						平均绝对误差	0.174444

由表4可以得出2014年中共性样本的序号为2，6，8，9，10，12，14，20，21。

把D2G与D2合并，共同构成2014年原烟叶组配方表中的个性化样本等级。

以2014年和2015年为例，对2014年与2015年的近红外光谱进行投影，寻找化学值或近红外光谱差异较大的等级，具体的计算方式如下：

以2014年的原烟叶组配方的个性化样本等级的光谱D1CS为依据，以2015年所收购的全体烟叶等级的光谱D2S为基础，从D2S中筛选m个等级，比例分别为k1,k2,……km，合成D1CS的光谱，比例k1,k2,……km，记作PN2。

对光谱矩阵D2S以及光谱矩阵D1CS进行必要的预处理，预处理后的近红外光谱分别记为PD2S与PD1CS。

以PD1CS的每条光谱为目标光谱，PD2S为基准光谱，以PD2S的光谱合成PD1CS的每条光谱，最后再把运用PD2S组合时的各个比例组合，构成PD1CS的个性化烟叶品种替代配方模块(即混合样)。

求解个性化烟叶品种的替代比例有两个关键的点，第一个点是如何替代，第二个点是替代的烟叶的具体比例该如何确定。针对第一个点，选择基础烟叶的光谱与目标烟叶的光谱的相似度作为替代的权重，而相似度的具体表现又分为光谱的趋势性与距离。针对第二个点，替代的比例是在上述权重下根据光谱的二次规划方程式来约束求解。

求解个性化烟叶品种的具体比例的方式如下：

设PD1CS-k为PD1CS的第k条光谱；计算PD2S每条光谱与PD1CS-k的相似度。

对[PD1CS-k，PD2S]进行主成分投影后，计算每条光谱到PD1CS-k的马氏(Mahalanobis)距离D_i ²；

${D_i}^2={(s_i-{\stackrel{\&OverBar;}{s}}_i)}^T{V}^{-1}(s_i-{\stackrel{\&OverBar;}{s}}_i)$

其中，D_i为第i个烟叶样本的马氏距离，si为第i个样烟叶本的主成分得分矩阵，为第i个样烟叶本的主成分得分矩阵的平均值，V为主成分得分矩阵的协方差。

个性烟叶的原始光谱矩阵为[PD1CS-k，PD2S]；对于PD2S的第i条光谱，记为PD2S-i；PD2S-i从光谱的起始波长点开始，选择近红外波长点窗口宽度W,在W内的PD2S-i_j-w+1:j光谱段与PD1S-k_j-w+1:j的局部相关系数分析R_j；也可以根据从光谱矩阵中随机抽取个若干个波长点组成X_i-w+1:i光谱，按照加权计算每个波长点下的典型相关系数CCAR_j，其中典型相关系数的原理与计算方法如下：

设两组随机变量X_i-w+1:i＝(x₁,x₂,...x_w),Y＝y分别为w和1维的向量，根据典型相关分析的思想，用X和Y的线性组合α'X和β'Y之间的相关性来研究两组随机变量X和Y之间的相关性，找到α和β，使得ρ(α‘X,β'Y)最大。由相关系数的定义：

$\mathrm{\&rho;}({\mathrm{\&alpha;}}^{\&prime;}X,{\mathrm{\&beta;}}^{\&prime;}Y)=\frac{Cov({\mathrm{\&alpha;}}^{\&prime;}X,{\mathrm{\&beta;}}^{\&prime;}Y)}{\sqrt{Var\left({\mathrm{\&alpha;}}^{\&prime;}X\right)}\sqrt{Var\left({\mathrm{\&beta;}}^{\&prime;}X\right)}}$

易得出对任意常数e,f,c,d，均有

ρ[e(α'X)+f,c(β'Y)+d]＝ρ(α'X,β'Y)

将局部相关系数或者局部典型相关系数合成整体的相关系数

或者

构造PD2S每条光谱到PD1CS-k的相似性度量(比例解析权重)W_k,i：

$W_{k,i}=\frac{{\mathrm{MR}}_i}{1+D_i}$

在计算PD1CS-k的具体比例时，对PD2S的光谱先作用

WPD2S_k＝PD2S×W_k

然后按照二次规划的的方式求解方程式

$\left\{\begin{array}{c}\left(\right|{\min }^{pi}*W^{I}pi+\frac{1}{2}{\mathrm{pi}}^T*H*pi\left|\right)\&RightArrow;\min \\ Aeq*pi=beq;\\ Ua\&le;pi\&le;Ub\end{array}\right.$

对于PD1CS-k，Ua为PD2S的具体每个等级使用的配方下限，Ub为PD2S的具体每个等级使用的配方上限。

C_k＝W_k,i*pi

C_k即为PD2S光谱合成PD1CS-k的具体使用比例；

求解PD1CS-k的所有C_k＝W_k,i*pi，k＝1,2,……m；其中m为PD1CS个性烟叶样本的个数。

以PD2S(2015年收购的所有烟叶的等级)的等级为原则，计算个性烟叶PD1CS的替代配方模块PD2S每个等级的具体比例C_pk。

把2015年的共性烟叶的使用比例PN1以及2014年PD1CS的替代配方模块PD2S每个等级的具体比例C_pk按照PD2S等级为原则，计算出2015年新的叶组配方的具体等级比例PC_pk。

2014年的其中一个个性烟叶样本的红外光谱如图4所示，利用2015年的烟叶样本的近红外光谱拟合图4的近红外光谱，拟合结果如图5所示，各烟叶样本的拟合比例C_pk如表5所示。

表5

(3)对于2015年的模拟叶组配方，按照2015年原烟叶组配方的最终比例PC_pk制成叶组小样；并比较与2014年实际叶组配方的化学值差异。2014年与2015年制样中的尼古丁组分的差异值如表6所示。

表6

2014年与2015年制样的近红外光谱差异如图6所示，感官评吸差异如表7所示。

表7

由表7可以看出，采用本发明方法得到的制样与2014年的片烟相比，香气略饱满厚实、甜香略足、刺激略大、略涩口，整体品质满足要求。

(4)如果差异不在预期范围之内，利用人工对原烟叶组配方进行微调。

Claims (5)

1.一种基于近红外光谱的原烟叶组配方维护方法，其特征在于，包括：

步骤1，针对去年原烟叶组配方中的烟叶品种以及今年采集的烟叶品种，分别进行近红外光谱采集以及化学值测定，依据近红外光谱和化学值的一致性程度，将去年的烟叶品种划分为共性等级和个性等级；

步骤2，针对共性等级，在今年的原烟叶组配方中采用与去年相同的比例；

针对个性等级，利用今年采集的烟叶品种组合成混合样，使混合样的近红外光谱以及化学值分别与对应的个性等级相一致，在今年的原烟叶组配方中，混合样的比例采用去年原烟叶组配方中对应的个性等级的比例，得到今年的初步原烟叶组配方；

步骤3，分别依据去年的原烟叶组配方和今年的初步原烟叶组配方进行制样，并对制样分别进行化学值测定和近红外光谱采集，若化学值和近红外光谱差异符合预期，得到今年的原烟叶组配方。

2.如权利要求1所述的基于近红外光谱的原烟叶组配方维护方法，其特征在于，还包括步骤4，对制样进行感官评析验证，感官评析符合预期后，确定今年的原烟叶组配方。

3.如权利要求2所述的基于近红外光谱的原烟叶组配方维护方法，其特征在于，步骤1中，若去年原烟叶组配方中的烟叶品种和今年采集的烟叶品种的化学值和近红外光谱差异均不大，则判定为共有等级；去年原烟叶组配方中的其余烟叶品种判定为个性等级。

4.如权利要求3所述的基于近红外光谱的原烟叶组配方维护方法，其特征在于，步骤1采集的化学值至少包括烟碱含量、总糖含量、还原糖含量、总氮含量、钾含量、氯含量。

5.如权利要求4所述的基于近红外光谱的原烟叶组配方维护方法，其特征在于，步骤1中划分共性等级的步骤包括：

选取去年和今年共有的烟叶品种，对共有的烟叶品种进行近红外光谱和化学值的比对，若近红外光谱和化学值的差异均不大，则判定为共有等级。