基于声学CT的温度场重建技术[1]是根据被测区域内多条路径上的声波飞行时间数据推算被测区域的温度分布，属于“由效果反求原因”的逆问题研究.该技术具有非接触不干扰被测温度场、测温范围广、环境适应能力强及可在线测量等优点.
重建算法在声学CT温度场重建中起到至关重要的作用.为实现温度场重建，通常要将被测区域划分成若干个网格，每个网格重建一个温度值，相当于一个原始温度采样点[2].典型的声学CT温度场重建算法按照正问题建模方式的不同主要有两大类：一类利用网格内声线路径长度建模，以最小二乘法(LSM)[3-5]为代表；另一类利用基函数对温度(声速)分布逼近建模，以Markov径向基函数Tikhonov正则化法(MTR)[6]为代表.LSM的优点是运算简单，重建速度快，但该算法要求被测区域网格剖分总数小于可获得的声波数据数；MTR算法的优点是被测区域网格剖分数不受声线数的限制，能够更有效地重建复杂温度场，但其重建质量仍然有较大的提升空间.卡尔曼滤波是一种实时递推算法，其实质是以最小均方误差为估计准则的最优估计，与最小二乘估计相比，卡尔曼滤波能够有效提升估计精度，抗噪能力也更为突出[7-8].
本文提出一种声学CT温度场重建新算法，利用Markov径向基函数逼近复杂的温度(声速)分布，建立正问题模型，利用卡尔曼滤波进行逆问题求解获得声速分布，进而利用温度与声速的关系得到温度场分布.
1 声学CT温度场重建原理
声波在气体介质中的传播速度c、气体介质的绝对温度T以及由气体介质组成成分所决定的声音常数B之间的关系[2]可以表示为
(1)
采用声学CT技术测量温度场需要在被测区域的边界处设置一定数量的声波收发器.各声波收发器依次发声，所有的收发器都可以接收该信号，从而在被测区域内部形成一定数量的声波传播路径.测量出声波在各声波传播路径上的飞行时间，利用某种重建算法以及声速与温度的对应关系，即可计算重建出被测区域内的温度分布.
2 重建BFAKF算法
声学CT重建属逆问题研究，通常要借助正问题模型来求解.本文提出的基于基函数逼近和卡尔曼滤波的重建算法(BFAKF)可分为正问题建模和逆问题求解两个环节.
2.1 基于基函数逼近的正问题建模
声学CT的正问题是在已知声波收发器的位置并确定有效声线后，根据温度(声速)分布，求解各路径上的声波飞行时间.
假设被测区域内声速倒数的分布可以表示为f(x，y，z)，则声波在第k条传播路径lk上的传播时间tk可以表示为
(2)
式中，N为有效声波路径数.将被测区域均匀地划分为M个网格，且M>N，利用M个Markov径向基函数的线性组合来表示f(x，y，z)，则有
(3)
式中：βi为待定系数；xi、yi、zi为第i个网格的中心点坐标；a为径向基函数的形状参数.将式(3)代入式(2)得
(4)
令A=(aki)k=1，2，…，N；i=1，2，…，M，t=[t1，t2，…，tN]T，β=[β1，β2，…，βM]T，可以得到用基函数逼近方法建立的声学CT温度场重建正问题模型，即
t=Aβ
(5)
2.2 基于卡尔曼滤波的逆问题求解
卡尔曼滤波[7-8]是一种递归的估计算法，在已知系统模型、噪声统计信息和初始状态的前提下，利用上一时刻的状态估计值和当前状态的观测值实现对当前时刻状态值的估计.假设考虑噪声的系统离散模型为
(6)
式中：Xk为k时刻的被估计状态；Φk/k-1为k-1时刻到k时刻状态转移矩阵；Ψk为系统噪声驱动矩阵；Wk为系统噪声；Hk为状态观测值；Zk为观测矩阵；Ok为观测噪声.
卡尔曼滤波由预测和校正两个过程组成.在预测阶段，用上时刻的状态估计
和估计误差协方差矩阵Pk-1求出当前时刻的状态预测
和预测误差协方差矩阵Pk/k-1.在校正阶段，利用对当前状态的观测值Zk修正在预测阶段获得的预测值，从而获得更接近真实值的新估计
和Pk.
针对式(6)的卡尔曼滤波方程可表示为
(7)
式中：Kk为滤波增益矩阵；Rk为观测噪声协方差矩阵.卡尔曼滤波正是通过不断调整对估计和观测的依赖程度，快速而有效地获得被估计量的最优估计值.
将卡尔曼滤波用于声学CT温度场重建逆问题求解时先要建立系统模型.假设测量声波飞行期间被测区域的温度(声速)不发生变化，即认为状态转移矩阵Φk/k-1为单位阵，系统噪声为0，则由k-1时刻到k时刻系统状态预测方程为
βk=βk-1
(8)
由于声波飞行时间测量过程中观测噪声的方差通常是比较稳定的，故假设观测噪声协方差矩阵Rk=R不随时间而变，则式(5)所对应的系统状态观测方程为
tk=Aβk+Ok
(9)
相应卡尔曼滤波方程可写为
(10)
在R已知的情况下，只要给定初值
就可根据该时刻的测量值tk，迭代计算出该时刻的状态估计值
没有先验知识时，估计误差协方差矩阵初值P0可取为bI，b取一个较大的值即可，本文取b=5.
3 重建误差评价指标
本文采用重建温度场的平均相对误差Rave、均方根误差Rrms以及热点温度相对误差Rh来评价重建质量，其定义分别为
(11)
(12)
(13)
式中：T(i)为模型温度场中第i个像素中心点的温度；Th为模型温度场中的热点温度；
为重建温度场中第i个像素的温度；
为重建温度场中的热点温度；M′为细化像素总数；Tave为像素中心点的平均温度.
4 重建算法仿真验证
本文将32个声波收发器均匀地布置在12 m×12 m×12 m的立方体区域周围，形成172条有效声波路径.采用最小二乘法(LSM)、Markov径向基函数Tikhonov正则化法(MTR)以及本文提出的基函数逼近卡尔曼滤波法(BFAKF)对四种模型温度场进行了仿真重建，各温度场表达式为
1) 热点位于(0，0，0)的单峰对称场
(14)
2) 热点位于(-2，-2，0)的单峰偏置场
(15)
3) 热点位于(-2，0，-2)和(2，0，2)的双峰对称场
(16)
4) 热点位于(-2，0，2)、(2，0，2)、(-2，0，-2)和(2，0，-2)的四峰对称场
(17)
采用LSM法重建时，被测区域均匀地划分为4×4×4=64个网格，即原始像素总数为64；采用MTR和BFAKF法重建时，被测区域均匀地划分为10×10×10=1 000个网格，即原始像素总数为1 000.MTR法的正则化参数μ=1×10-6，形状参数a=1×10-4；BFAKF法的形状参数a=1×10-4，估计误差协方差矩阵初值P0=5I，迭代次数为50，飞行时间数据无噪声时取观测噪声协方差矩阵R0=1×10-5I，飞行时间数据有噪声(噪声标准差δ=2×10-5)时取R0=7×10-5I.
图1以单峰对称和四峰对称模型温度场为例，给出了有噪声情况下BFAKF法的均方根误差、平均相对误差随迭代次数的变化曲线.表1、2分别给出了采用无噪声飞行时间和有噪声飞行时间时LSM法、MTR法以及BFAKF法对应的重建误差.由于LSM法无法正确地重建出四峰对称温度场特征，故本文不给出四峰对称温度场的LSM法热点温度相对误差.算法的重建时间分别是：LSM算法0.023 61 s、MTR算法0.116 46 s、BFAKF算法2.157 2 s.LSM法只有64个原始像素，正问题模型中的系数矩阵A只有172×64个元素，MTR法有1 000个原始像素，矩阵A有172×1 000个元素，所以逆问题求解需要更长的时间.
图1 误差随迭代次数变化曲线
Fig.1 Curves of errors varying with iteration number
表1 采用无噪声飞行时间和原始像素时的重建误差
Tab.1 Reconstruction errors with noise-free time-of-flight and original pixels %
算法重建误差单峰对称单峰偏置双峰对称四峰对称MTR算法BFAKF算法LSM算法Rave1.6962.1302.2552.653Rrms2.05012.87010.819.420Rh1.22712.87022.55,22.5520.40,20.40,20.40,20.40Rave1.2861.6821.9052.589Rrms1.6136.7904.2034.940Rh0.9269.8467.889,7.88910.65,10.65,10.65,10.65Rave3.0508.9797.4235.914Rrms3.13853.52012.395.622Rh2.63150.0200.083,0.083-,-,-,-
表2 采用有噪声飞行时间和原始像素时的重建误差
Tab.2 Reconstruction errors with noisy time-of-flight and original pixels %
算法重建误差单峰对称单峰偏置双峰对称四峰对称MTR算法BFAKF算法LSM算法Rave1.8582.1902.3532.764Rrms2.29412.90010.789.544Rh0.74819.75022.23,22.5320.58,20.44,20.37,20.03Rave1.6511.9152.1352.570Rrms2.1339.5647.0926.694Rh0.20814.36014.63,14.8015.27,15.13,14.83,14.64Rave3.0758.9487.4805.858Rrms3.11153.23012.485.556Rh2.92750.3408.326,8.139-,-,-,-
表3给出了飞行时间数据无噪声时LSM、MTR以及BFAKF法对应的重建温度场与模型温度场.为用尽可能少的篇幅显示温度场特征，单峰温度场只给出了x=0 m、y=0 m、z=0 m对应的切片图；双峰和四峰温度场只给出了y=0 m处切片图.由于仅用64或者1 000个像素难以细致地描述一个复杂的温度场，故而表4、5给出了细化像素对应的重建误差.具体做法是先将重建出的温度值赋予原始像素的中心点，然后将12 m×12 m×12 m的区域细化成31×31×31=29 791个像素，再用三次样条插值的方法求出细化像素中心点的温度值.LSM法的原始像素数M=4×4×4=64，细化像素数M′=23×23×23=12 167，细化像素的描述区域为9 m×9 m×9 m；MTR和BFAKF法的原始像素数M=10×10×10=1 000，细化像素数M′=27×27×27=19 683，细化像素的描述区域为10.8 m×10.8 m×10.8 m.
表3 模型温度场与重建温度场对比
Tab.3 Comparison between model and reconstruction temperature fields
模型温度场MTR法重建温度场BFAKF法重建温度场LSM法重建温度场单峰对称单峰偏置双峰对称四峰对称
5 结 论
本文提出一种利用基函数逼近和卡尔曼滤波的声学CT温度场重建算法.采用BFAKF算法以及有代表性的LSM、MTR算法分别对四种典型的温度场模型进行了仿真数据重建，重建结果表明：BFAKF法能使重建误差以较快的速度下降并趋于稳定；与LSM、MTR法相比，BFAKF法对应的重建误差更小，重建温度场与模型温度场更接近，BFAKF算法具有更好的复杂温度场重建能力.但BFAKF算法的重建时间较长.一方面是因为该算法采用了较多的原始像素导致正问题模型中的系数矩阵维数较大；另一方面由于该算法是迭代算法，系数矩阵维数大必然会导致逆问题求解时间长.
表4 采用无噪声飞行时间和细化像素时的重建误差
Tab.4 Reconstruction errors with noise-free time-of-flight and refined pixels %
算法重建误差单峰对称单峰偏置双峰对称四峰对称MTR算法BFAKF算法LSM算法Rave1.8232.3932.6193.042Rrms2.06513.71011.279.634Rh1.14128.85030.83,30.8325.20,25.20,25.20,25.20Rave1.4421.9672.2642.835Rrms1.6508.2845.3245.450Rh1.06019.26016.47,16.4717.84,17.84,17.84,17.84Rave5.67862.58022.7517.55Rrms6.67893.26030.9023.36Rh10.52021.82051.99,51.99-,-,-,-
表5 采用有噪声飞行时间和细化像素时的重建误差
Tab.5 Reconstruction errors with noisy time-of-flight and refined pixels %
算法重建误差单峰对称单峰偏置双峰对称四峰对称MTR算法BFAKF算法LSM算法Rave1.9382.4462.6753.124Rrms2.30113.76011.239.750Rh1.11228.75030.80,30.6925.27,25.27,25.24,25.07Rave1.7232.1922.5102.898Rrms2.15310.7107.8197.044Rh1.00323.64023.21,23.1921.67,21.63,21.13,21.06Rave5.59762.43022.8417.47Rrms6.61892.96030.9423.31Rh10.71021.96051.97,51.93-,-,-,-
[1] 顾梦楠.声学CT温度场测量算法重建研究 [D].沈阳：沈阳工业大学，2019.
(GU Meng-nan.Research on reconstruction algorithms of acoustic CT for temperature field measurement [D].Shenyang：Shenyang University of Technology，2019.)
[2] Yan H，Chen G N，Zhou Y，et al.Primary study of temperature distribution measurement in stored grain based on acoustic tomography [J].Experimental Thermal and Fluid Science，2012，42(5)：55-63.
[3] 董晨龙，周新志，白兴都，等.基于Reflected-Sigmoid径向基函数插值的温度场重建算法 [J].四川大学学报(自然科学版)，2019，56(5)：851-856.
(DONG Chen-long，ZHOU Xin-zhi，BAI Xing-du，et al.Temperature field reconstruction algorithm based on Reflected-Signmoid radial basis function interpolation [J].Journal of Sichuan University (Natural Science Edition)，2019，56(5)：851-856.)
[4] Zhang S，Shen G，An L，et al.Online monitoring of the two-dimensional temperature field in a boiler furnace based on acoustic computed tomography [J].Applied Thermal Engineering，2015，75：958-966.
[5] Ru X J，Qing Y X.Two-dimensional temperature field distribution reconstruction based on least square method and radial basis function approximation [J].Mathe-matical Problems in Engineering，2017，17：1-7.
[6] Yan H，Ma Z，Zhou Y G.An acoustic tomography system for online monitoring of temperature fields [J].IET Science Measurement & Technology，2017，11(5)：623-630.
[7] 李华，李莉，郭育艳.哈希编码和Kalman滤波的目标跟踪算法 [J].沈阳工业大学学报，2019，41(4)：406-411.
(LI Hua，LI Li，GUO Yu-yan.Object tracking algorithm based on Hash coding and Kalman filtering [J].Journal of Shenyang University of Technology，2019，41(4)：406-411.)
[8] 张立峰，李佳，田沛.Kalman滤波在电容层析成像图像重建中的应用 [J].计量学报，2017，38(3)：315-318.
(ZHANG Li-feng，LI Jia，TIAN Pei.Application of Kalman filter in image reconstruction of capacitance tomography [J].Metrology Journal，2017，38(3)：315-318.)