基于稀疏重构的声学CT温度场重建

魏元焜,颜 华,周英钢

(沈阳工业大学信息科学与工程学院,辽宁沈阳 110870)

0 引言

声学CT温度场重建技术[1]具有非接触、不干扰被测温度场、测量范围广、环境适应能力强、维护方便等优点,作为一种新型工业测温仪表的关键技术,具有广泛的应用前景。该技术可用于声学高温计,用来监测锅炉炉膛、烟道燃烧室的温度分布[1-2]。

以声学CT重建温度场时,需要在被测区域周围设置多个声波收发器,形成M条从不同方向穿越被测区域的声波路径(简称声线),并测得声波在各声线上的飞行时间(time of flight,TOF)。根据温度分布求TOF是声学CT温度场重建的正问题,由测得的TOF反推温度分布是其逆问题。为求解逆问题,通常将测温区域划分成N个网格,利用正问题建立数学模型,然后利用重建算法根据TOF反推N个网格内的声慢(声速倒数)值,最后以插值算法获得细密的温度分布。经典的重建算法有最小二乘法(least square method,LSA)[3-4]、联合迭代重建技术(simultaneous iterative reconstruction technique,SIRT)[5-6]等。实际应用中可获得的TOF数据非常有限,因此N大于M的重建算法更适合于复杂声慢度(温度)分布的重建。但N大于M时逆问题是欠定的(即线性方程组的方程数少于未知数个数),且通常N值越大,声学CT投影矩阵条件数越大,即重建结果对TOF的测量误差越敏感。为解决这一问题,取得稳定的逆问题解,经典算法往往要求N

为提升声学CT温度场重建精度,应对N>M导致的逆问题求解的不稳定,本文提出了不同的思路,提出了基于稀疏重构(sparse reconstruction)的声学CT温度场重建算法,简称SR算法。

1 SR算法

1.1 数学模型的建立

以二维声学CT温度场重建为例考虑声学CT正问题,设测温平面内声慢度分布函数为f(x,y),则被测区域内第j(1

(1)

式中:fi为网格i内的声慢度值;aji为声线pj在网格i内的长度。

同时考虑M条声线,定义t=(t1,…,tM)T,A=(aji),f=(f1,…,fN)T,j=1,…,M,i=1,…,N,则声学CT系统获得TOF的过程可表示为

t=Af

(2)

式中:f为N×1维向量,由N个网格内的声慢值构成,可看作系统的观测对象;A为M×N维矩阵,为描述系统的观测矩阵,表征了声学CT系统对声慢度的观测过程。

分析式(1)、式(2)可知,矩阵A的元素aji可视为声线pj对网格i内声慢度信号的采样;M×1维向量t为系统的观测结果,由M条声线上的声波TOF构成,可由声学CT系统测得。收发器位置固定后,矩阵A为已知量。因此在式(2)所示模型下,声学CT逆问题的实质是根据已知的t和A求f。

由压缩感知理论[15-17]可知,只要信号在某个变换域是稀疏的,那么就可以由一个与变换基不相干的观测矩阵所获得的少量观测数据,然后利用稀疏重构算法极高概率地重构出原信号。这一过程与声学CT温度场重建相似,可作为N>M的逆问题求解思路。但声速不可能为无穷大,所以声慢度不会是0,即声慢度向量f不是稀疏的,因此需要按式(3)将f变换到稀疏域:

f=Ψθ

(3)

式中:Ψ为N×N维矩阵,通常Ψ是一个标准正交基;θ为N×1维向量,且θ中大部分元素为0或非常接近0。

称Ψ为f的稀疏基,θ为f在Ψ域的稀疏表示,θ中非零元素的个数为其稀疏度。常用的稀疏基有离散傅里叶变换(DFT)基和离散余弦变换(DCT)基,本文研究发现声慢信号在DCT域中具有更好的稀疏性,本文以DFT作为SR算法的稀疏基Ψ。

基于f在DFT域的稀疏性,将式(3)带入式(2),即可得到基于稀疏重构的声学CT数学模型:

t=AΨθ=Uθ

(4)

式中:θ为逻辑上的被测稀疏信号;U为等效观测矩阵,U=AΨ。

(5)

式中z为由气体组分决定的常数,空气中z=20.05。

最后利用三次样条插值算法,即可获得更加细致像素下的温度场描述。

在稀疏重构中有2个重要环节——观测矩阵的优化和稀疏重构,将在下文展开论述。

1.2 观测矩阵的优化

由压缩感知理论可知,为实现稀疏信号θ的高精度重构,观测矩阵A和稀疏基Ψ的相干性要尽可能低。高斯随机矩阵因其随机性,在很大概率上与任意稀疏基不相干,因此压缩感知中常采用独立同分布的高斯随机矩阵作为观测矩阵[18-19]。但声学CT系统的结构化特质使其无法实现高斯随机矩阵,因此需要对矩阵A进行优化,使其尽量接近高斯随机矩阵的随机特性。本文通过提高声线穿过各个网格的“无差别化”和“无偏颇性”,使声线在各网格内的分布更充分、更均匀,从而等效提高观测矩阵A对声慢度信号采样的随机性。

为定量评价声线在各网格内布局的充分、均匀程度,定义N×1维向量Γ如式(6)所示:

(6)

式中:Γi为网格i内所有声线长度之和,i=1,2,…,N;当测温区域和网格划分确定后,aji为各收发器位置坐标构成的向量X的函数。

定义声线布局充分均匀度F(X)如式(7)所示,

F(X)=μF(Γ)-σF(Γ)

(7)

式中:X为各收发器位置坐标构成的向量;μF(Γ)为向量Γ各元素的均值,反映声线布局的充分性,其数值越大充分性越好;σF(Γ)为向量Γ各元素的标准差,反映声线布局的均匀性,其数值越小均匀性越好。

F(X)数值越大,声线布局充分均匀性越好,F(X)数值越小,声线充分均匀性越差。在给定被测区域及N值和M值的情况下,F(X)仅由收发器位置坐标向量X决定。

本文以求F(X)最大值为优化目标,用遗传算法求解X。种群规模设置为200,最大代数为300,迁移率为0.2,杂交率为0.8。设置被测区域为边长1.3 m的正方形,划分N=15×15=225个网格。4个边框上各布置3个声波收发器,形成M=54条有效穿越被测区域的声线。由于收发器固定在正方形的边框上,因此每个收发器必有一个坐标是常数,要优化的收发器位置坐标X可缩减为12×1维向量。

图1(a)给出了常用的收发器标准布局[3,4,10,13]以及对应的声线。该布局也被用作遗传算法的初值。图1(b)给出了用遗传算法获得的优化布局。优化前、后的声线布局充分均匀度F(X)分别为0.122 7和0.138 7。相对于优化前,优化后的充分均匀度提高11.5%。

图1 优化前后的收发器布局

1.3 稀疏信号的重构

根据式(4)由已知的t和U求θ的过程可归结为式(8)所示的l0范数最小化问题[19]。

(8)

式(8)可由OMP算法[15-16]求解,具体流程如下:

(1)初始化,令残差r0=t,索引向量Λ0为空,支撑矩阵V0为空,迭代计数m=1。

(2)列选择,计算矩阵U各列与残差rm-1的内积。将内积最大的列向量加入支撑矩阵Vm-1,构成新的支撑集矩阵Vm;将选中的列向量在矩阵U中的列索引号加入索引向量Λm-1,构成新的索引向量Λm;最后,以零向量替换U中被选中列。

(3)更新残差,以最小二乘法求解式(9),计算θm并以式(10)更新残差。

(9)

rm=t-Vmθm

(10)

(4)确定是否满足迭代终止条件,若满足,则停止迭代并跳转到步骤(5);否则令m=m+1,返回步骤(2)。

2 温度场重建实验

在自行设计的声学CT温度场重建装置上实施了重建实验,被测区域为边长1.3 m的正方形。实验装置的系统框图如图2所示,图2中关键模块实物如图3所示。计算机控制12个扬声器(TC9FD-19-08)轮流发声,12个传声器(CHZ 221)将接收到的声波信号转换成电信号,经信号调理器后送入数据采集器(PCI 6123),利用互相关时延估计算法可获取各声线上的TOF,重建程序调用重建算法重建测温层面上的温度分布。

图2 硬件系统框图

图3 硬件实物图

实验选取了4种声学CT重建算法:本文提出的SR算法,其观测矩阵对应图1(b)所示的优化布局,记为SR-b;采用了标准收发器布局的SR算法,其观测矩阵对应图1(a)所示的标准布局,记为SR-a算法;采用了标准收发器布局的LSA算法以及SIRT算法,其观测矩阵对应图1(a)所示的标准布局。SR算法采用N=15×15的网格划分,LSA和SIRT为避免逆问题求解的不稳定,采用了N=5×5的网格划分。SIRT和LSA算法按文献[5,7]实现,其中SIRT算法的松弛因子为0.3,当相邻2次迭代结果之差相对上次迭代结果之比小于1×10-4时,停止SIRT迭代。SR算法估计σ参数的方法如下:数据采集卡PCI 6123的最大采样率为500 kHz,由此可先确定采样间隔所对应的最大量化误差σs=2 μs;然后在均温场TF1下连续进行了40次TOF测量,计算出所有声线TOF标准差的平均值σm=1.1 μs;最后取σs和σm中的数值较大者作为

TOF所含噪声标准差σ的估计值,即σ=2×10-6。

以4种算法分别重建301.2 K的均温场(TF1)和单峰对称温度场(TF2)。为保证标准布局和优化布局下TF2的近似性,电热炉的位置保持不变。利用空调调节室内温度,在保证加热前测温区温度均为301.2 K的情况下,开启电热炉。图4给出了两种收发器布局下的TF2温度场。

(a)标准布局

(b)优化布局图4 2种收发器布局下的单峰温度场TF2

图5给出了4种算法对均温场TF1、单峰场TF2的重建图像。由图5可以看出:4种算法的重建结果均能反映实际温度场的主要特征;对于TF1而言,SR-b和SR-a算法重建的图像呈现出非常好的均匀性,而SIRT和LSA算法重建的图像有明显的非均匀性;对于TF2,SR-b算法重建图像中伪影的颜色与背景的颜色最接近,而其他3种算法在被测区域的四角或边界附近伪影的颜色与背景颜色差异明显。

(a)SR-b重建TF1

(b)SR-a重建TF1

(c)LSA重建TF1

(d)SIRT重建TF1

(e)SR-b重建TF2

(f)SR-a重建TF2

(g)LSA重建TF2

为对重建温度场进行误差评价,通常以热电偶进行接触式测量。由于在测试平面布置大量的热电偶会影响声学CT的精确重建,因此重建误差评价通常只针对关键点进行[10-13]。本文将80PK-1型热电偶安装在高精度手持式数字温度计(Fluke 54 II)上实现温度的实时显示。对于TF1温度场,用热电偶测量了被测层面的最高温度Tmax和最低温度Tmin。对于TF2温度场,用热电偶测量了被测层面上电热炉上方的最高温度。表1给出了热电偶测量和算法重建温度的对比。重建温度相对于热电偶测量值的误差e也在表1中给出。为便于整体评价,表1中的最后一行还给出了同一算法所对应的|e|的平均值。由表1可以看出:SR-b算法的重建温度与热电偶的测量温度最为接近,SR-a算法与LSA算法重建温度与热电偶测温偏差逐渐增大;从|e|的平均值来看,SR-b的重建误差平均值最小,SR-b、SR-a的重建误差明显低于SIRT、LSA算法,SR-b和SR-a的重建误差低于LSA误差最为明显。

表1 热电偶测量与算法重建温度对比

温度场TF2是在301.2 K的均温场基础上通过电热炉局部加热形成的,因此TF2温度场中所有像素(三次样条插值的细密像素37 pixel×37 pixel)的真实温度都不应出现低于301.2 K的异常结果。但受重建误差的影响,4种算法的TF2重建图像中都存在温度值低于301.2K的“异常”像素。显然,异常像素数目越少,异常像素的温度越接近301.2 K,重建温度场越接近实际温度场,则算法性能越好。表2给出了4种算法所对应的异常像素的数目Nabn、异常像素温度的平均值Tavgabn和最小值Tminabn。由表2可以看出,4种算法中SR-b算法的Nabn最小,Tavgabn和Tminabn与301.2 K也最接近。

表2 异常像素统计

综合图5、表1和表2可知,4种算法中SR-b算法的重建质量最优,LSA算法的重建质量最差,SR-a算法的重建质量优于SIRT、LSA算法。因此从重建误差来看,本文提出的稀疏重构算法优于经典的LSA和SIRT算法,且本文的观测矩阵优化方法可以有效降低重建误差。

3 结束语

本文提出的SR算法利用声慢度信号在DFT域的稀疏性,在稀疏重构的框架中求解声学CT逆问题,允许N>M的声学CT温度场重建。通过温度场重建实验可以得出以下结论:

(1)SR算法重建温度场的误差低于LSA算法和SIRT算法,证实了稀疏重构算法在温度场重建应用中的可行性与优越性。

(2)SR-b算法重建误差低于SR-a算法,说明通过优化收发器布局,可以有效提高声学CT采样随机性,通过满足压缩感知理论中对采样随机性的要求,可以有效提高基于稀疏重构的声学CT温度场重建质量。

(3)SR算法中采用的噪声水平相关的OMP迭代终止条件在温度场重建实验中,表现出了良好的稳定性,表明SR算法在室内实际应用中可以稳定地重建温度场。