荧光寿命测温多指数数据处理方法

Multi-Exponential Data Processing Method in Fluorescence Lifetime Temperature Measurement

摘要: 基于荧光寿命的光纤测温技术关键点在于将荧光信号中的直流分量与噪声信号等背景干扰去除，尽可能准确地求得荧光寿命τ，但在高精度测量时，运用传统的单指数模型算法，如最小二乘拟合、积分面积比、FFT拟合等方法处理多指数的荧光信号时，求取荧光寿命 会带来很大误差。将LM算法和Prony算法引入荧光测温中，通过仿真分析可知，二者在处理多指数的荧光信号时较传统方法有明显优势，LM算法在低噪条件下，处理双指数模型的最大误差在±10⁻⁴ ms之间，响应速度3 × 10⁻¹ s，适合用于高精度测量之中；Prony算法受噪声影响较为明显，在低噪环境下的最大误差为±10⁻³ ms，响应速度可以达到2 × 10⁻³ s，适合应用在快速响应条件下。

Abstract: The key point of optical fiber temperature measurement technology based on fluorescence lifetime is to remove the background interference such as the dc component and the noise signal in the fluorescence signal so as to obtain the fluorescence life τ as accurately as possible. However, in the case of high-precision measurement, when the traditional single exponential model algorithms, such as the least squares fitting, the integral area ratio, the FFT fitting and other methods are used to process the multi-exponential fluorescence signal, the fluorescence life error obtained will be large. LM algorithm and Prony algorithm are introduced into fluorescence temperature measurement. Through simulation analysis, it can be seen that the two algorithms have obvious advantages over the traditional methods in processing multi-exponential fluorescence signals. Under low noise condition, the maximum error of LM algorithm is between ±10⁻⁴ ms, and the re-sponse speed is 3 × 10⁻¹ s. It is very suitable for high-precision measurement. Prony algorithm is significantly affected by noise. The maximum error is ±10⁻³ ms in the low noise environment, and the response speed can reach 2 × 10⁻³ s. It is suitable for application in the case of fast response.

1. 引言

荧光测温是一种可以应用在工业、农业、医疗等行业的温度测量技术，因其抗电磁干扰的特点 [1] [2]，较传统热电偶测温有更大优势。同时基于荧光寿命的光纤测温系统有着不受激励光强干扰的特点，具有很好的发展前景 [3]。

实际应用时荧光寿命的理论衰减形式为：

$I\left(t\right)=I_0\exp \left(-t/\tau \right)+I_d+I_n$ (1)

其中，I_d为直流分量，I _n为噪声信号，二者均属于背景干扰 [4] [5]。 $\tau$ 为荧光寿命，是温度的单值函数，求取 $\tau$ 值即可得到被测温度值 [6] [7]。为保证求 $\tau$ 的准确性，需要对接收的荧光信号进行去直流，消除背景噪声等预处理，通常采用最小二乘拟合、积分面积比、FFT拟合等方法 [8]。

但在高精度测量时发现，荧光信号为多指数形式，即：

$I\left(t\right)=I_0\exp \left(-t/{\tau }_1\right)+I_1\exp \left(-t/{\tau }_2\right)+\cdots +I_N\exp \left(-t/{\tau }_J\right)+I_d+I_n$ (2)

此时传统针对单指数的各算法存在弊端，需寻找更准确的算法来进行提高测量的精度。

LM算法和Prony算法是针对多指数模型的数据处理方法 [9] [10]。本文通过仿真分析，研究将LM算法和Prony算法引入荧光测温的可行性和测量精度的问题。

2. L-M算法

2.1. L-M算法在测温系统中的实现

定义待测荧光信号的模型为：

$I\left(t\right)={\displaystyle \underset{i=1}{\overset{k}{\sum }}{A}_i\exp \left(-t/{\tau }_i\right)}$ (3)

在这里，构造一个代价函数来估计衰减残差的平方和：

$\begin{array}{c}S\left(\beta ;t\right)=\frac{1}{2}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\left[I_j\left(\beta ;t\right)-y_j\right]}^2}\\ =\frac{1}{2}{\displaystyle \underset{j=1}{\overset{n}{\sum }}{f}_j{\left(x\right)}^2}\\ \equiv \frac{1}{2}F{\left(x\right)}^{\text{T}}F\left(x\right)\end{array}$ (4)

式(4)中的 $F\left(x\right)$ 是由 $f\left(x\right)$ 组成的向量值函数， $\beta$ 是由A_i和 $\tau$ 组成的参数。这样一来，求取荧光寿命 $\tau$ 就变成了求取参数 $\beta$ 来使得S最小化的问题。考虑到S不是呈线性关系，需要更加精确、快速的收敛方法，因此提出了使用LM算法来进行拟合。

LM算法的迭代格式为：

$x_{k+1}=x_k+d_k$ (5)

其中d_k为搜索方向，在高斯牛顿法的迭代条件基础上，Levenberg和Marquardt提出了通过对角截断来对J^TJ矩阵进行衰减，即引进一个阻尼因子 $\lambda$，使其能够在较大的初始范围内进行收敛：

$x_{k+1}=x_k-{\left(J^{\text{T}}J+\lambda \text{diag}\left(J^{\text{T}}J\right)\right)}^{-1}JF\left(x\right)$ (6)

式中：J为I对 $\beta$ 的梯度；J^TJ为Hessian矩阵；用Hessian的对角矩阵来代替单位矩阵，可以减小计算机字长短带来的累积舍入误差的影响，提高计算结果的精度。

阻尼因子 $\lambda$ 在每次迭代过程中会进行自我调整，当S减小时， $\lambda$ 相应减小，当 $\lambda \to 0$ 时，上式近似为高斯牛顿法；当S增大时， $\lambda$ 相应增大，当 $\lambda \to \infty$ 时，上式近似于最速下降法。当目标函数的精度达到预设值时，迭代结束。因此只要选择合适的迭代初值，就可以兼顾高斯牛顿法和最速下降法的优点，使得迭代收敛速度快，全局稳定性高。

在仿真过程中，使用LM算法分别对单指数模型和双指数模型进行了分析。

2.2. L-M算法单指数仿真分析

针对单指数模型，此时荧光信号为：

$I\left(t\right)=A_1\exp \left(-t/\tau \right)+I_d+I_n$

利用MATLAB做一条单指数函数曲线，设其模型为：

$y=A\exp \left(-i\Delta t/\tau \right)+I_d+I_n$

采样时间间隔 $\Delta t=0.01\text{ms}$。假设荧光寿命 $\tau =0.\text{5}\text{ms}$，通过仿真反求 $\tau$ 值偏离0.5 ms的大小来分析LM算法对直流分量和噪声分量的敏感程度及算法处理速度。

2.2.1. 直流分量的影响

令 $A=1$，噪声信号 $I_n=0$，改变直流分量I_d的幅度与荧光信号强度的比，结果如表1所示。可以看出，利用LM算法求取荧光寿命 $\tau$ 时，会受到直流分量的影响，随着直流分量的增大，荧光寿命 $\tau$ 也逐渐偏离理论值。因此若想应用此法得到较高的测量精度需要想办法去掉直流分量，这是此法的一个弊端。

Table 1. The τ value table of the DC component is 0.5% - 2.5% of the fluorescence intensity

表1. 直流分量为荧光强度的0.5%~2.5%时的 $\tau$ 值表

2.2.2. 噪声信号的影响

令直流分量 $I_d=0$，改变噪声信号I_n的幅度与荧光信号强度的比，结果如表2所示。可以看出，LM算法对噪声不敏感，是一种精度较高的算法。采样点的选取不仅决定了数据处理速度，还与噪声有关，采样长度并不是越大越好，例如在0.5%的幅度比情况下，采样长度120的仿真结果最接近理论值，因此要根据实际情况确定采样长度。

Table 2. The τ value table of the noise amplitude is 0.5% - 2.5% of the fluorescence intensity

表2. 噪声幅度为荧光强度的0.5%~2.5%时的 $\tau$ 值表

2.3. LM算法双指数仿真分析

用MATLAB做一条双指数曲线，设其模型为：

$y=A_1\exp \left(-i\Delta t/{\tau }_1\right)+A_2\text{exp}\left(-i\Delta t/{\tau }_2\right)+I_d+I_n$

假设荧光寿命 ${\tau }_{\text{1}}=0.\text{5}\text{ms}$， ${\tau }_{\text{2}}=0.\text{3}\text{ms}$，通过仿真反求两个 $\tau$ 值偏离程度来进行对比分析。

噪声信号的影响

通过单指数模型仿真分析可知，算法受直流分量的影响，因此此处讨论无直流情况下对噪声的敏感程度。选择采样长度为120，改变噪声幅度与荧光强度的比例，得到表3所示的结果。

Table 3. The τ value table of the fitting length of l-m algorithm is 120 points

表3. L-M算法拟合长度为120点时的 $\tau$ 值表

这里对比传统的最小二乘拟合法和FFT拟合法，分别在采样长度相同的情况下求取荧光寿命 ${\tau }_{\text{1}}$ 进行对比，结果如表4所示。

Table 4. The τ value table of the fitting length of traditional algorithm is 120 points

表4. 传统算法拟合长度为120点时的 $\tau$ 值表

对比以上两表数据可知，传统最小二乘法和FFT拟合法等算法在针对双指数模型时 $\tau$ 值计算效果较差，精度很低，误差最小在0.12 ms，受噪声影响明显，并不适合处理多指数荧光信号；而LM算法的误差维持在 0⁻⁴ ms到10⁻³ ms，测量精度较高，但由于算法需要使用到大量的迭代技术，对初值的猜测也有一定的要求，导致运行时间较长。针对这种情况，提出另一种扩充Prony算法进行改进。

3. Prony算法

3.1. Prony算法在测温系统中的实现

Prony算法是使用多指数函数的一种线性组合来描述等间距采样数据的数学模型，可以根据采样值直接估算出信号频率、衰减、幅值和初相位的分析方法。

该算法针对等间距采样点，在荧光测温中，实际光谱的荧光寿命曲线可看做纯指数或一系列指数函数的叠加，形式为：

$\begin{array}{c}I\left(t\right)=I_1\exp \left(-t/{\tau }_1\right)+I_2\exp \left(-t/{\tau }_2\right)+\cdots +I_n\exp \left(-t/{\tau }_n\right)\\ ={\displaystyle \underset{j=1}{\overset{n}{\sum }}{I}_{0j}\exp \left(-t/{\tau }_j\right)}\end{array}$ (7)

这样，当把采样值带入上式时，即可得到联立方程组：

$\{\begin{array}{l}{I}_0=I_{01}+I_{02}+\cdots +I_{0n}\\ I_1=I_{01}\exp \left(-\Delta t/{\tau }_1\right)+I_{02}\exp \left(-\Delta t/{\tau }_2\right)+\cdots +I_{0n}\exp \left(-\Delta t/{\tau }_n\right)\\ I_2=I_{01}\exp \left(-2\Delta t/{\tau }_1\right)+I_{02}\exp \left(-2\Delta t/{\tau }_2\right)+\cdots +I_{0n}\exp \left(-2\Delta t/{\tau }_n\right)\\ ⋮\\ I_{2n-1}=I_{01}\exp \left(-\left(2n-1\right)\Delta t/{\tau }_1\right)+I_{02}\exp \left(-\left(2n-1\right)\Delta t/{\tau }_2\right)+\cdots +I_{0n}\exp \left(-\left(2n-1\right)\Delta t/{\tau }_n\right)\end{array}$ (8)

即：

$I_i=I\left(i\Delta t\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}{I}_{0j}\exp \left(-i\Delta t/{\tau }_j\right)}$ (9)

其中，I_i为采样值，Δt为采样间隔，求解该非线性联立方程组，即可求得荧光寿命 $\tau$ 和幅值I_0j，求解该方程组非常复杂，因此引入Prony算法，

设测量数据的估计模型为：

${\hat{I}}_n={\displaystyle \underset{i=1}{\overset{p}{\sum }}{b}_i{z}_i^n}\left(n=0,1,\cdots ,N-1\right)$ (10)

其中更一般的表示模型有：

$b_m=A_m$ (11)

$z_m=\exp \left(a_m\Delta t\right),a=-1/\tau$ (12)

这里，A_m为振幅，a_m为衰减因子，Δt为采样间隔，prony算法采用平方误差最小原则来使模拟信号向真实信号进行逼近，即：

$\min \left(\epsilon ={\displaystyle \underset{n=0}{\overset{N-1}{\sum }}{\left|I\left(n\right)-\hat{I}\left(n\right)\right|}^2}\right)$ (13)

同时考虑到估计模型函数是一个常系数线性差分方程的齐次解，因此定义了多项式：

$\phi \left(z\right)={\displaystyle \underset{t=0}{\overset{p}{\prod }}\left(z-z_k\right)}={\displaystyle \underset{i=0}{\overset{p}{\sum }}{a}_i{z}^{p-1}}$ (14)

改造式(13)有：

$\hat{I}\left(n-m\right)={\displaystyle \underset{i=1}{\overset{p}{\sum }}{b}_i{z}^{n-m}}\left(0\le n-m\le N-1\right)$ (15)

将该式乘a_m并对p + 1个乘积求和：

${\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m\hat{I}\left(n-m\right)}={\displaystyle \underset{l=1}{\overset{p}{\sum }}{b}_l{\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m{z}_l^{n-m}}},p\le n\le N-1$ (16)

代入 $z_l^{n-m}=z_l^{n-p}{z}_l^{p-m}$，有：

${\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m\hat{I}\left(n-m\right)}={\displaystyle \underset{l=1}{\overset{p}{\sum }}{b}_l{z}_l^{n-p}}{\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m{z}_l^{p-m}}=0$ (17)

该式等于零正是因为第二项求和恰好是 $\phi \left(z\right)={\displaystyle \underset{t=0}{\overset{p}{\prod }}\left(z-z_k\right)}={\displaystyle \underset{i=0}{\overset{p}{\sum }}{a}_i{z}^{p-1}}$ 位于根z_i处的多项式 $\varphi \left(z_i\right)$，而

$\varphi \left(z_i\right)=0$。因此可知 $\hat{I}\left(n\right)$ 满足递推的差分方程：

$\hat{I}\left(n\right)=-{\displaystyle \underset{m=1}{\overset{p}{\sum }}\hat{I}\left(n-m\right)}\left(p\le n\le N-1\right)$ (18)

假设测量值 $\stackrel{⌢}{I}\left(n\right)$ 与真实值I(n)之间关系为： $I\left(n\right)=\stackrel{⌢}{I}\left(n\right)+e\left(n\right)$，其中为e(n)一个加性高斯白噪声，这样就有了：

$\begin{array}{c}I\left(n\right)=-{\displaystyle \underset{m=1}{\overset{p}{\sum }}{a}_m\hat{I}\left(n-m\right)+e\left(n\right)}\\ ={\displaystyle \underset{m=1}{\overset{p}{\sum }}{a}_{m}I\left(n-m\right)}+{\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_{m}e\left(n-m\right)}\left(p\le n\le N-1\right)\end{array}$ (19)

其中 $a_0=\text{1}$，上式表明白噪声中的指数过程是一个ARMA(p,p)模型，它具有相同的AR和MA参数，且激

励噪声为原加性白噪声e(n)。参数的最小二乘估计使 ${\displaystyle {\sum }_{n=p}^{N-1}{\left|e\left(n\right)\right|}^2}$ 最小，这又将得到一组困难的非线性方程。

估计AR参数的线性方法定义：

$\epsilon \left(n\right)={\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_{m}e\left(n-m\right)}\left(n=p,\cdots ,N-1\right)$ (20)

这样式(19)变为：

$I\left(n\right)=-{\displaystyle \underset{m=1}{\overset{p}{\sum }}{a}_{m}I\left(n-m\right)}+\epsilon \left(n\right)$ (21)

所以求使 ${\displaystyle \underset{n=p}{\overset{N-1}{\sum }}{\left|e\left(n\right)\right|}^2}$ 最小就转换成了求使 ${\displaystyle \underset{n=p}{\overset{N-1}{\sum }}{\left|\epsilon \left(n\right)\right|}^2}$ 最小，这就是扩充的Prony算法，即扩充的Prony算法

就是求解下列矩阵方程：

$\left[\begin{array}{cccc}x\left(p\right)& I\left(p-1\right)& \cdots & I\left(0\right)\\ x\left(p+1\right)& I\left(p\right)& \cdots & I\left(1\right)\\ ⋮& ⋮& \ddots & ⋮\\ x\left(p\right)& I\left(p\right)& \cdots & I\left(p\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}\epsilon \left(0\right)\\ \epsilon \left(1\right)\\ ⋮\\ \epsilon \left(N-1\right)\end{array}\right]$

或

$Ia=\epsilon$ (22)

为使 ${\epsilon }_p={\displaystyle \underset{n=p}{\overset{N-1}{\sum }}{\left|\epsilon \left(n\right)\right|}^2=}{\displaystyle \underset{n=p}{\overset{N-1}{\sum }}{\left|{\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_{m}I\left(n-m\right)}\right|}^2}$ 最小，求 $\partial {\epsilon }_p/\partial a_m$ 令其等于零，有

${\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m\left[{\displaystyle \underset{n=p}{\overset{N-1}{\sum }}I\left(n-m\right)I^{\ast }\left(n-i\right)}\right]=0}$ (23)

对应的最小误差能量：

${\epsilon }_p={\displaystyle \underset{m=0}{\overset{p}{\sum }}{a}_m\left[{\displaystyle \underset{n=p}{\overset{N-1}{\sum }}I\left(n-m\right)I^{\ast }\left(n\right)}\right]}$ (24)

在这里定义：

$r\left(i,j\right)={\displaystyle \underset{n=p}{\overset{N-1}{\sum }}I\left(n-j\right)I^{\ast }\left(n-i\right)}$ (25)

则上述偏导式可写为

$\left[\begin{array}{cccc}r\left(0,0\right)& r\left(0,1\right)& \cdots & r\left(0,p\right)\\ r\left(1,0\right)& r\left(1,1\right)& \cdots & r\left(1,p\right)\\ ⋮& ⋮& \ddots & ⋮\\ r\left(p,0\right)& r\left(p,1\right)& \cdots & r\left(p,p\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}{\epsilon }_p\\ 0\\ ⋮\\ 0\end{array}\right]$ (26)

求解这一方程即可得到AR参数 $a_1,\cdots ,a_p$ 和最小误差能量 ${\epsilon }_p$ 的估计值。进一步可以求解如下特征方程的根Z：

$1+a_1{z}^{-1}+\cdots +a_p{z}^{-p}=0$ (27)

这样，上述联立方程组模型就可以简化为：

$\left[\begin{array}{cccc}1& 1& \cdots & 1\\ z_1& z_2& \cdots & z_p\\ ⋮& ⋮& \ddots & ⋮\\ z_1^{N-1}& z_2^{N-2}& \cdots & z_p^{N-1}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_p\end{array}\right]=\left[\begin{array}{c}\hat{x}\left(1\right)\\ \hat{x}\left(2\right)\\ ⋮\\ \hat{x}\left(p\right)\end{array}\right]$ (28)

其中，幅值A_i为和衰减因子a_i (即荧光寿命 $\tau$ )可表示为：

$\{\begin{array}{l}{A}_i=\left|b_i\right|\\ a_i=\ln \left(z_i\right)/\Delta t\end{array}$ (29)

由上推导，实现了Prony算法在荧光测温中的实现，现使用MATLAB软件分别对单指数模型和双指数模型进行分析。

3.2. Prony算法的单指数仿真

用MATLAB软件做一条单指数荧光曲线，设其模型为：

$y=A\exp \left(-i\Delta t/\tau \right)+I_d+I_n$，

采样时间间隔 $\Delta t=0.0\text{1}\text{ms}$。假设荧光寿命 $\tau =0.\text{5}\text{ms}$，通过仿真反求 $\tau$ 值偏离 .5 ms的大小进行分析。

针对单指数模型，令 $A=\text{1}$，噪声信号 $I_n=0$，采样长度选择120，改变直流分量I_d的幅度与荧光信号强度，结果如表5所示。可以看出，Prony算法同样会受到直流分量的背景干扰，在直流分量较小时，仿真结果较好。

Table 5. The τ value table of the DC component is 0.5% - 2.5% of the fluorescence intensity (sampling point 120)

表5. 直流分量为荧光强度的0.5%~2.5%时的 $\tau$ 值表(采样点120)

再令直流分量 $I_d=0$，采样长度选择120，改变噪声信号幅度与荧光信号幅度的比，结果如表6所示。可以看出，Prony算法对噪声很敏感，随着噪声的增加，仿真结果也越来越偏离理论值，在低噪环境下精度较好。

Table 6. The τ value table of the noise amplitude is 0.5% - 2.5% of the fluorescence intensity (sampling point 120)

表6. 噪声幅度为荧光强度的0.5%~2.5%时的 $\tau$ 值表(采样点120)

3.3. Prony算法的双指数仿真

用MATLAB做一条双指数荧光曲线，设其模型：

$y=A_{}\exp \left(-i\Delta t/{\tau }_1\right)+A_2\exp \left(-i\Delta t/{\tau }_2\right)+I_d+I_n$

假设荧光寿命 ${\tau }_{\text{1}}=0.\text{5}\text{ms}$， ${\tau }_{\text{2}}=0.\text{3}\text{ms}$，探究噪声对其的干扰。设拟合数据120点，改变噪声幅度与荧光强度的比例，结果如表7所示。由此可知在低噪环境下，Prony算法表现较好，仿真时长在10⁻³ s数量级，比起LM算法更适合于快速响应环境下；同时其受噪声干扰较大，当噪声干扰较大时，精度不如LM算法。

Table 7. The τ value table of the noise amplitude is 0.5% - 2.5% of the fluorescence intensity (sampling point 120)

表7. 噪声幅度为荧光强度的0.5%~2.5%时的 $\tau$ 值表(采样点120)

4. 结论

论证了LM算法与Prony算法可以运用在荧光测温系统中，比起传统数据处理方法，在处理多指数函数模型时有着明显优势。LM算法在低噪条件下，处理双指数模型的最大误差在±10⁻⁴ ms之间，响应速度3 × 10⁻¹ s，非常适合用于高精度测量之中；Prony算法受噪声影响较为明显，在低噪环境下的最大误差为±10⁻³ ms，相应速度可以达到2 × 10⁻³ s，适合应用在快速响应条件下。两种算法各有优势所在，具有较好发展前景。