CN110514608A - 一种基于光谱的反应动力学速率常数的无偏估计方法 - Google Patents

Abstract

本发明提供了一种基于光谱的反应动力学速率常数的无偏估计方法，属于化学反应动力学领域。该方法首先构建光谱数据结构，根据反应物初始浓度和反应动力学速率常数的初始估计值，计算浓度分布曲线，获得将光谱矩阵，与测量得到的光谱矩阵比较，得到误差矩阵；然后根据误差矩阵调整反应动力学速率常数的估计值，重新计算误差矩阵，重复迭代，直到误差达到最小值；最后用最终误差矩阵估计噪声方差，在测量得到的光谱中增加噪声，重新计算反应动力学速率常数估计值，使用无偏估计方法，得到的结果更加接近真实值。由于使用无偏估计方法，不依赖反应物质的纯光谱等其他信息，有效的减少了噪声产生的影响，使得到的反应动力学速率常数更加精确。

Description

一种基于光谱的反应动力学速率常数的无偏估计方法

技术领域

本发明属于化学反应动力学领域，具体地涉及一种基于光谱的反应动力学速率常数的无偏估计方法。

背景技术

建立精确的化学反应过程的动力学模型，以便对化学反应过程进行合理的预测和有效的控制，既能保证整个生产过程的安全性，也能提高化工生产效率，降低成本，最大化经济效益，对于现代化学工业具有重大意义。精确的反应动力学模型包括了模型结构和模型参数。确立模型结构不仅需要对化学反应过程的反应机理十分了解，还需要根据反应器的类型，所研究问题的特点进行合理的假设。同样十分重要的是如何优化模型参数，使得模型模拟得出的数据与实际值偏差最小，是当前化学反应动力学研究的一个重点。

经过多年的发展，大量学者从优化设定及最优控制等方面对模型参数估计展开了研究。基于光谱拟合动力学模型参数的方法也不断发展，并针对影响结果精度的因素提出了相应的策略。经典曲线分析方法是一种非线性最小二乘拟合方法，该方法是一种梯度方法，从初始值出发依据梯度信息迭代到最优值。使用该方法可以较为快速地从光谱中估计出反应速率常数，但是没有对实验误差和测量噪声进行进一步的研究。为了补充该缺陷，利用光谱估算反应速率常数的约束最小二乘法也被提出，但是该方法在参与反应物质的纯光谱难以获得的时候，效果大打折扣。因此现有的反应动力学速率常数估计方法仍存在缺陷，需要改进。

发明内容

为了解决现有的反应动力学参数估计方法存在的不足，本发明提出一种基于光谱的反应动力学速率参数的无偏估计方法。

本发明解决其技术问题所采用的技术方案是：一种基于光谱的反应动力学速率常数的无偏估计方法，具体包括以下步骤：

步骤一：基于Beer-Lambert定律构建光谱数据的结构，使用常微分方程组建立反应动力学模型，根据反应物的初始浓度和反应动力学速率常数初始估计值，计算化学反应中涉及的所有物质随时间变化的浓度分布曲线，利用最小二乘计算各物质的吸光系数，计算光谱数据矩阵，将计算得到的光谱数据矩阵与实验测量的光谱比较，得到误差矩阵E(k₀)。

步骤二：根据步骤一所述误差矩阵E(k₀)调整反应动力学速率常数初始估计值，计算得到新的浓度分布曲线，利用最小二乘计算各物质的吸光系数，计算新的光谱数据矩阵，将计算得到的新光谱数据矩阵与实验测量的光谱比较，重新计算得到新的误差矩阵E(k₁)。经上述方法重复迭代，获得最终误差矩阵E(k_n-1)，直到最终误差矩阵E(k_n-1)中每个元素的平方和与上一次误差矩阵E(k_n-2)中每个元素的平方和之差小于10^-6。n为迭代次数，n≥1。

步骤三：根据步骤二得到的最终误差矩阵E(k_n-1)中每个元素的平方和估计噪声的方差，在所述测量得到的光谱中增加具有所述方差的噪声，获得具有叠加噪声的光谱数据，重新计算反应动力学速率常数估计值，使用无偏估计方法，完成对反应动力学速率常数的估计。

进一步地，步骤一具体为：

所述光谱数据矩阵D为：D＝CS^T+E(k_n-1)，其中，C表示具有光吸收特性的物质随时间变化的浓度矩阵，S是每种物质在不同测量波长下的吸光系数矩阵，D表示反应溶液对于每种测量波长的吸光程度随时间的变化，E(k_n-1)与D维度相同。

以时间为自变量的常微分方程组建立反应动力学模型，给定光谱数据矩阵D和所有反应物的初始浓度，根据反应动力学速率常数的初始估计值k₀，计算出C、S、E(k_n-1)，其中，C⁺为C的伪逆：

E(k_n-1)＝D-CS^T

S^T＝C⁺D

C⁺＝(C^TC)^-1C^T。

进一步地，步骤二包括以下子步骤：

(2.1)E(k₀)中每个元素E_i,j的平方和ssq的求解过程为：

nt是采样点的数量，nλ为测量波长的数量。

(2.2)以ssq向其最小值移动的方式计算反应动力学速率常数k的偏移Δk，反应动力学速率常数初始值表示为k₀，经过Δk＝(Δk₁,Δk₂,…,Δk_nk)的微小偏移后，用泰勒展开来估计E(k₀)，其中E(k₀)的泰勒展开式中含有偏导数矩阵。

(2.3)将E(k₀)和所述偏导数矩阵矢量化，其中，矢量化后的偏导数矩阵称为雅可比矩阵J，通过下列公式求解Δk，e(k₀)为E(k₀)矢量化后的向量：

e(k₀)＝-JΔk

J⁺＝(J^TJ)^-1J^T

Δk＝-J⁺e(k₀)

(2.4)将计算所得的Δk与k₀相加，得到新的反应动力学速率常数估计值，重复步骤(2.1)-(2.3)，直到E(k_n-1)的ssq和与E(k_n-2)的ssq之差小于10^-6，循环终止，获得此时的反应动力学速率常数k(D)。

进一步地，所述步骤三包括以下子步骤：

(3.1)计算步骤(2.4)中得到的速率常数k(D)所对应的E(k_n-1)含有的所有元素的方差σ²，视为实验测量噪声的方差。

(3.2)在测量光谱数据矩阵D上加上均值为0、方差为σ²的高斯白噪声，记为Z，重新计算反应动力学速率常数的估计值k(Z)，得到无偏估计结果为k＝2k(D)-k(Z)。

(3.3)重复步骤(3.2)，将结果取平均，得到最终反应动力学速率常数无偏估计结果。

本发明的有益效果是：本发明从光谱出发，使用经典曲线分析方法实现对反应动力学速率常数的估计。由于实验中存在测量噪声，得到的结果含有一定的误差，因此使用无偏估计方法，不依赖反应物质的纯光谱等其他信息，有效的减少了噪声产生的影响，使得到的反应动力学速率常数更加精确。

附图说明

图1为本发明方法的流程示意图；

图2为步骤(2.2)中有限差分自适应步长选择方法的流程示意图。

具体实施方式

为了使本发明的目的和技术方案更加清楚明白，以下结合附图和实施例,对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

参照图1，一种基于光谱的反应动力学速率常数的无偏估计方法，包括以下步骤：

步骤一：基于Beer-Lambert定律构建光谱数据的结构，使用常微分方程组建立反应动力学模型，根据反应物的初始浓度和反应动力学速率常数初始估计值，计算化学反应中涉及的所有物质随时间变化的浓度分布曲线，计算化学反应中涉及的所有物质随时间变化的浓度分布曲线，利用最小二乘计算各物质的吸光系数，得到光谱数据矩阵，将计算得到的光谱矩阵与实验测量的光谱比较，得到误差矩阵E(k₀)。

步骤二：根据步骤一所述误差矩阵调整反应动力学速率常数初始估计值，计算得到新的浓度分布曲线，利用最小二乘计算各物质的吸光系数，得到新光谱数据矩阵，将计算得到的新光谱数据矩阵与所述测量得到的光谱比较，重新计算得到新的误差矩阵E(k₁)。经上述方法重复迭代，获得最终的误差矩阵E(k_n-1)，直到最终的误差矩阵E(k_n-1)中每个元素的平方和与上一次误差矩阵E(k_n-2)中每个元素的平方和之差小于10^-6。n为迭代次数，n≥1。

步骤三：根据步骤二得到的E(k_n-1)中每个元素的平方和估计噪声的方差，在所述实验测量的光谱中增加具有所述方差的噪声，获得具有叠加噪声的光谱数据，重新计算反应动力学速率常数估计值，使用无偏估计方法，完成对反应动力学速率常数的估计。

其中，步骤一的方法具体为：根据Beer-Lambert定律，时间为t时，在特定的波长λ处，溶液的吸光度是来自所有吸收该波长光线的物质的吸光总和。每种物质对吸光度的贡献与浓度和光程长度(入射光束在溶液中行进的距离)成线性比例。如果路径长度和比例常数由单个常数s_i(λ)表示，则Beer-Lambert定律由以下公式表示：

d(t,λ)＝c₁(t)s₁(λ)+c₂(t)s₂(λ)+…+c_nc(t)s_nc(λ)

nc为在波长λ处吸光物质的数量，c_i(t)是时间为t时i物质的浓度，d(t,λ)表示时间为t时，溶液对波长为λ的光的吸收总和。跟踪化学反应的变化，获得的光谱即为反应时间的函数，测量结果为二维光谱数据矩阵：D＝CS^T+E(k_n-1)，其中，C表示每种物质随时间变化的浓度，为nt×nc的矩阵；S是每种物质在不同测量波长下的吸光系数，为nλ×nc的矩阵；D为nt×nλ光谱数据矩阵，表示反应溶液对于每种测量波长的吸光程度随时间的变化；由于测量存在误差，引入E(k_n-1)，与D维度相同。

本发明研究的化学反应，假定反应器内各种物料均匀混合，该假设下研究对象的模型用以时间为自变量的常微分方程组(ordinary differential equations，ODES)建立反应动力学模型，给定光谱数据矩阵D和所有参与反应物质的初始浓度，根据反应动力学速率常数的初始估计值k₀，结合ODES，计算出C、E(k_n-1)，其中，C⁺为C的伪逆：

E(k_n-1)＝D-CS^T

S^T＝C⁺D

C⁺＝(C^TC)^-1C^T

步骤二还包括以下子步骤：

(2.1)使用经典曲线分析方法，计算E(k₀)中每个元素E_i,j的平方和ssq(sum ofsquares，ssq)的求解过程为：

nt是采样点的数量，nλ为测量波长的数量。

(2.2)以ssq向其最小值移动的方式计算反应动力学速率常数k的偏移Δk，为此，需要强调的是E(k₀)和随后的ssq仅是反应动力学速率常数k的函数，反应动力学速率常数初始值表示为k₀，经过Δk＝(Δk₁,Δk₂,…,Δk_nk)的微小偏移后，用泰勒展开来估计E(k₀)，其中E(k₀)的泰勒展开式中含有偏导数矩阵，仅使用一阶导数，得到如下线性表达式：

虽然上述公式为近似表示方法，由于其为线性表达式，这使得它易于处理。目标是使E(k₀+Δk)的值趋于0，因此用0代替E(k₀+Δk)，上式重新排列结果如下：

误差矩阵E(k₀)可以由上式表达，并且偏导数可以根据下式方法计算：

这里使用的是前向差分，本文根据Gill提出的算法进行选择Δk_i，算法流程参照图2，步骤详解如下：

ε_A＝|f(x)|ε_R

其中，h为有限差分自适应步长，为函数f(x)的前向差分，为函数f(x)的后向差分，为函数f(x)的中心差分，Φ(f,h)为函数f(x)的二阶导数，ε_R为机器精度，ε_A为最小步长，为第一判别变量，为第二判别变量，为第三判别变量。

步骤1：确认函数f(x)与自变量x，选择参数M，令m＝0，h_s＝-1， 如下式：

计算式(9)中的各个参数；如果则h_s←h₀。

步骤2：m＝m+1，h_m＝10h_m-1，重新计算式(9)中的各个参数。如果h_s<0并且则h_s←h_m。如果满足则h_Φ←h_m，步骤2停止，进入下一步；否则重复步骤2直到m＝M，进入下一步；

步骤3：m＝m+1，h_m＝h_m-1/10，重新计算式(9)中的各个参数。如果则h_s←h_m。如果则h_Φ←h_m，步骤3停止，进入下一步；否则重复步骤3直到m＝M，步骤3停止，进入下一步；

步骤4：h_F＝h_Φ，输出h_F，算法终止；

步骤5：如果h_s<0，则取如果则取h_F←h_s；如果则取h_F←h_M。输出h_F，算法终止。

(2.3)将E(k₀)和所述偏导数矩阵矢量化，即展开成长列向量，将所述长列向量折叠成矩阵与向量的乘积；矢量化后的偏导数矩阵称为雅可比矩阵J，通过下列公式求解Δk，e(k₀)为E(k₀)矢量化后的向量：

e(k₀)＝-JΔk

J⁺＝(J^TJ)^-1J^T

Δk＝-J⁺e(k₀)

(2.4)将计算所得的Δk与k₀相加，得到新的反应动力学速率常数估计值，重复步骤(2.1)-(2.3)，直到E(k_n-1)的ssq和与E(k_n-2)的ssq之差小于10^-6，循环终止，获得此时的反应动力学速率常数k(D)。

步骤三还包括以下子步骤：

(3.1)计算步骤(2.4)中得到的速率常数k(D)所对应的E(k_n-1)含有的所有元素的方差σ²，视为实验测量噪声的方差。

(3.2)在测量光谱数据矩阵D上加上均值为0、方差为σ²的高斯白噪声，记为Z，重新计算反应动力学速率常数的估计值k(Z)，得到无偏估计结果为k＝2k(D)-k(Z)。

(3.3)需要强调的是，无偏性的实际意义是指没有***性的偏差。无论测量几次同一化学反应过程，经过参数估计得到的反应速率常数相对于真实值而言，总有一些偏高，一些偏低。而无偏性表示，把这些误差平均起来，其值为零，即无偏估计量只存在随机误差。因此重复步骤(3.2)150次，将结果取平均，得到最终反应动力学速率常数无偏估计结果。

实施例

(1)根据Beer-Lambert定律，时间为t时，在特定的波长λ处，溶液的吸光度是来自所有吸收该波长光线的物质的吸光总和。每种物质对吸光度的贡献与浓度和光程长度(入射光束在溶液中行进的距离)成线性比例。跟踪化学反应的变化，获得的光谱即为反应时间的函数，测量结果是一个二维矩阵，即为光谱数据矩阵：D＝CS^T+E(k_n-1)，其中，C表示每种物质随时间变化的浓度，S是每种物质在不同测量波长下的吸光系数，D为光谱数据矩阵，表示反应溶液对于每种测量波长的吸光程度随时间的变化，E(k_n-1)表示误差矩阵，与D维度相同。

以时间为自变量的ODES建立反应动力学模型，这里选用MCR-ALS工具包中的演示算例，其ODES如下：

反应速率常数真实值k_r＝[2,0.2]，物质初始浓度[A]₀＝0.001mol·L^-1，其余物质为0。采样时间t＝10s，t₀＝0，采样点个数nt＝300。测量光谱矩阵D已知，反应动力学参数的初始估计值设为k₀＝[1,0.5]，结合ODES可以计算出三种物质浓度随时间分布的矩阵C，误差矩阵E的计算方法如下，其中C⁺为C的伪逆：

E＝D-CS^T

S^T＝C⁺D

C⁺＝(C^TC)^-1C^T

(2.1)使用经典曲线分析方法，计算所述误差矩阵E(k₀)中每个元素E_i,j的平方和ssq的求解过程为：

nt是采样点的数量，nλ为测量波长的数量。

(2.2)一旦确定了ssq，就以ssq向其最小值移动的方式计算反应速率常数的偏移。为此，需要强调的是E(k₀)和随后的ssq仅是速率常数k的函数，如果速率常数初始值表示为k₀，则参数经过Δk＝(Δk₁,Δk₂,…Δk_nk)的微小偏移之后，可以用泰勒展开来估计E(k₀)，如果仅使用一阶导数，得到如下线性表达式：

其中为E(k₀)对于k_nk的偏导数矩阵。目标是使E(k₀+Δk)的值趋于0，因此用0代替E(k₀+Δk)，上式重新排列结果如下：

误差矩阵E(k₀)可以由上式表达，并且偏导数可以根据下式方法计算：

(2.3)将E(k₀)和所述偏导数矩阵矢量化，即展开成长列向量，将所述长列向量折叠成矩阵与向量的乘积；矢量化后的偏导数矩阵称为雅可比矩阵J，通过下列公式求解Δk，e(k₀)为E(k₀)矢量化后的向量：

e(k₀)＝-JΔk

J⁺＝(J^TJ)^-1J^T

Δk＝-J⁺e(k₀)

(2.4)将计算所得的Δk与k₀相加，得到新的反应动力学速率常数估计值，重复步骤(2.1)-(2.3)，直到E(k_n-1)的ssq和与E(k_n-2)的ssq之差小于10^-6，循环终止，获得此时的反应动力学速率常数。

通常，如上所述方法在最小值附近快速收敛。然而，如果初始估计较差，则泰勒级数展开的函数逼近和问题的线性化可能导致ssq的发散和算法的失败。Marquardt基于Levenberg的思想提出的修改是在参数计算过程中将一定数量的Marquardt参数mp加到Hessian矩阵的对角元素上。Hessian矩阵及Δk新的表达式如下：

H＝J^TJ

Δk＝-(H+mp*I)^-1J^Te(k₀)

当mp的值显著大于Hessian矩阵的对角元素时，Δk表达式如下，具有线性下降的性质：

当mp很小时，Δk表达式不变。管理参数mp策略如下：最初设置为1，如果参数ssq开始发散，则mp不断增加(每次迭代乘以10)，ssq收敛则不断减小(每次迭代除以√10)，在迭代结束后重置为1。

(3.1)在步骤(2.4)中得到的速率常数记为k(D)＝[1.9766,0.19476]，计算它所对应的E(k_n-1)含有的所有元素的方差σ²＝4.3，可以视为实验测量噪声的方差；

(3.2)在测量光谱矩阵D上加上均值为0，方差为σ²的高斯白噪声，记为Z，重新计算反应动力学速率常数的估计值k(Z)，得到无偏估计结果为k＝2k(D)-k(Z)。

(3.3)重复步骤(3.2)150次，结果取平均，算作一次无偏估计计算。以真实值k_r＝[2,0.2]为标准，平均误差率可以作为评价一个结果是否优秀的指标。

未考虑测量噪声所得的结果为k(D)＝[1.9766,0.19476]，平均误差率为1.90％。使用了无偏估计方法之后，在100次的重复实验中，有69次k(Z)比k(D)更接近真实值，部分结果如表1。说明在此算例中，无偏估计能够较为有效地减少噪声带来的误差。

表1 MCR-ALS算例仿真无偏估计部分结果

以上显示和描述了本发明的基本原理和主要特征和本发明的优点，但本发明并不仅仅限于说明书和实施方式所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限制的一般概念下，本发明并不限于特定的细节和这里示出的实施例。

Claims (4)

1.一种基于光谱的反应动力学速率常数的无偏估计方法，其特征在于：具体包括以下步骤：

步骤一：基于Beer-Lambert定律构建光谱数据的结构，使用常微分方程组建立反应动力学模型，根据反应物的初始浓度和反应动力学速率常数初始估计值，计算化学反应中涉及的所有物质随时间变化的浓度分布曲线，利用最小二乘计算各物质的吸光系数，计算光谱数据矩阵，将计算得到的光谱数据矩阵与实验测量的光谱比较，得到误差矩阵E(k₀)。

步骤二：根据步骤一所述误差矩阵E(k₀)调整反应动力学速率常数初始估计值，计算得到新的浓度分布曲线，利用最小二乘计算各物质的吸光系数，计算新的光谱数据矩阵，将计算得到的新光谱数据矩阵与实验测量的光谱比较，重新计算得到新的误差矩阵E(k₁)。经上述方法重复迭代，获得最终误差矩阵E(k_n-1)，直到最终误差矩阵E(k_n-1)中每个元素的平方和与上一次误差矩阵E(k_n-2)中每个元素的平方和之差小于10^-6。n为迭代次数，n≥1。

步骤三：根据步骤二得到的最终误差矩阵E(k_n-1)中每个元素的平方和估计噪声的方差，在所述测量得到的光谱中增加具有所述方差的噪声，获得具有叠加噪声的光谱数据，重新计算反应动力学速率常数估计值，使用无偏估计方法，完成对反应动力学速率常数的估计。

2.根据权利要求1所述反应动力学速率常数的无偏估计方法，其特征在于，步骤一具体为：

所述光谱数据矩阵D为：D＝CS^T+E(k_n-1)，其中，C表示具有光吸收特性的物质随时间变化的浓度矩阵，S是每种物质在不同测量波长下的吸光系数矩阵，D表示反应溶液对于每种测量波长的吸光程度随时间的变化，E(k_n-1)与D维度相同。

以时间为自变量的常微分方程组建立反应动力学模型，给定光谱数据矩阵D和所有反应物的初始浓度，根据反应动力学速率常数的初始估计值k₀，计算出C、S、E(k_n-1)，其中，C⁺为C的伪逆：

E(k_n-1)＝D-CS^T

S^T＝C⁺D

C⁺＝(C^TC)^-1C^T。

3.根据权利要求1所述反应动力学速率常数的无偏估计方法，其特征在于，步骤二包括以下子步骤：

(2.1)E(k₀)中每个元素E_i,j的平方和ssq的求解过程为：

nt是采样点的数量，nλ为测量波长的数量。

(2.2)以ssq向其最小值移动的方式计算反应动力学速率常数k的偏移Δk，反应动力学速率常数初始值表示为k₀，经过Δk＝(Δk₁,Δk₂,…,Δk_nk)的微小偏移后，用泰勒展开来估计E(k₀)，其中E(k₀)的泰勒展开式中含有偏导数矩阵。

(2.3)将所述E(k₀)和所述偏导数矩阵矢量化，其中，矢量化后的偏导数矩阵称为雅可比矩阵J，通过下列公式求解Δk，e(k₀)为E(k₀)矢量化后的向量：

e(k₀)＝-JΔk

J⁺＝(J^TJ)^-1J^T

Δk＝-J⁺e(k₀)

(2.4)将计算所得的Δk与k₀相加，得到新的反应动力学速率常数估计值，重复步骤(2.1)-(2.3)，直到E(k_n-1)的ssq和与E(k_n-2)的ssq之差小于10^-6，循环终止，获得此时的反应动力学速率常数k(D)。

4.根据权利要求1所述反应动力学速率常数的无偏估计方法，其特征在于，所述步骤三包括以下子步骤：

(3.1)计算步骤(2.4)中得到的速率常数k(D)所对应的E(k_n-1)含有的所有元素的方差σ²，视为实验测量噪声的方差。

(3.2)在测量光谱数据矩阵D上加上均值为0、方差为σ²的高斯白噪声，记为Z，重新计算反应动力学速率常数的估计值k(Z)，得到无偏估计结果为k＝2k(D)-k(Z)。

(3.3)重复步骤(3.2)，将结果取平均，得到最终反应动力学速率常数无偏估计结果。