基于PCA 降维和迭代正则化的温度场重建

白 云，颜 华，魏元焜

（沈阳工业大学 信息科学与工程学院，沈阳 110870）

声学CT 温度场重建技术[1-6]具有非接触、环境适应力强、维护方便、实时监测等优点[7-8]，备受工业界青睐。重建系统的性能很大程度上依赖于重建算法的精度。依正问题建模方式的不同重建算法主要有两大类。第一类基于网格内声线长度建模[9]，典型算法有最小二乘法（简称LSM 法）[10]，网格数N 通常小于声线数M。针对边缘信息缺失等问题，文献[11]利用均衡优化算法得到粗网格下的温度分布，再利用高斯过程回归算法预测细密网格下的温度分布。第二类采用径向基函数逼近声慢度分布的方式建模，允许N>M，用正则化等手段求解严重病态的逆问题[12]。例如文献[13]采用截断的广义奇异值分解法。

重建区域划分的网格数越多，对复杂温度场的描述能力越强，但重建矩阵的奇异性显著增加。针对此问题本文提出一种基于主成分分析（principal component analysis，PCA）降维和迭代正则化的重建算法。该算法采用第二类建模方式，用PCA 降维改善逆问题的病态性，用迭代正则化法求解逆问题，进而得到温度分布。

1 声学CT 温度场重建原理

声波在气体介质中的传播速度c 与气体介质的绝对温度间T 间的关系为

式中：z 是声音常数，由气体组成成分决定，对空气而言有z=20.045。

用声学CT 法重建温度场需事先将多组声波收发器布置在待测区域周围。系统工作时依据收发器采集到的声波数据，计算声波在收发器间的飞行时间，通过预设的重建算法，计算被测层面的声慢度（声速的倒数）分布，进而利用式（1）得到温度分布。

2 基于PCA 降维和迭代正则化的重建算法

2.1 正问题建模

利用布置在被测层面周围的收发器在被测层面形成M 条声线，并将被测区域离散成N 个网格。设被测区域内的声慢度分布为s（x，y），第i 个网格中心点的坐标为（xi，yi），则可用式（2）计算声波在第k 条声线lk上的飞行时间：

本文利用N 个径基函数的线性组合逼近声慢度分布，如式（3）所示：

式中：βi为待定系数；α＞0 为径向基函数的形状参数，根据所需测量范围以及收发器位置摆放，利用数值实验确定此参数。令：

将式（4）代入式（2），有：

考虑M 条路径并令A=（Aki）k=1，…，M，i=1，…，N，t=（t1，…，tM）T，β=（β1，…，βN）T，建立正问题模型如式（6）所示：

式中：β为N 维待定系数向量；A为M×N 维重建矩阵；t为M 维声波飞行时间向量。本文M=54，N=81。

2.2 逆问题求解

2.2.1 基于主成分分析的降维处理

为获得良好的重建质量，需要获得尽可能多的反映温（声慢）度分布信息的声波飞行时间数据。相应的，重建矩阵也成为典型的高维度数据，有很强的行相关性。为减轻冗余特征带来的多重共线性问题，本文通过对重建矩阵A的主成分分析，对正问题模型进行降维处理。

主成分分析[14-16]通过对原始的M 维变量的正交变换，得到被称为主成分的线性无关的新变量。第一个主成分具有最大的信息量（方差值），后续主成分包含的信息量依次降低。前K（K

（1）将矩阵A的每一个行向量零均值化处理，得到矩阵Y，而后求矩阵Y的协方差矩阵C：

（2）先求矩阵C的特征值和特征向量，而后由上到下排列前K 大的特征值所对应的特征向量，形成K×N 维转换矩阵P；

（3）令Adr=PA得到降维后的K×N 维重建矩阵Adr。

本文K 取47，P为47×54 维转换矩阵。去掉的7 个特征值累积贡献率仅为0.0001%，降维所带来的信息损失完全可以忽略。重建矩阵A条件数为1.68×108，降维重建矩阵Adr的条件数为8.77×107，重建矩阵的条件数降低47.8%。

用变换矩阵P左乘式（6）两端，并令tdr=Pt，可得到PCA 降维后的正问题模型如式（8）所示：

2.2.2 逆问题求解

降维处理可有效地改善重建矩阵的病态性，但为了避免信息损失所带来的重建失真，降维后Adr的仍具有严重的病态性。为此，本文提出的重建算法采用了式（9）所示的迭代正则化获得式（8）的稳定解：

式中：α 为松弛因子；μ 为正则化参数；I为单位矩阵；k为迭代次数。求出β后，用式（3）可以计算出N 个网格中心点处的声慢度。由于声慢度是声速的倒数，故利用式（1）可得到N 个网格中心点处的温度值。

3 重建误差评价

本文采用均方根误差Rrms对重建质量进行评估：

4 重建算法的仿真数据验证

本文的被测层面为1.3 m×1.3 m 的正方形，每条边长的四等分点上布置一个声波收发器，形成54 条穿越被测区域的声线，如图1 所示。采用本文提出的基于PCA 降维和迭代正则化的重建算法（简称pcaIR 法）对单峰、双峰和三峰模型温度场分别进行仿真重建。单峰模型温度场的热点位于（0，0），双峰模型温度场的热点位于（-0.25，-0.25）和（0.25，0.25），三峰模型温度场的热点位于（-0.35，-0.35），（-0.35，0.35）和（0.2，-0.1）。作为对比算法，基于奇异值分解的直接正则化算法（简称svdDR法）和最小二乘算法（LSM）的重建结果也在本文中给出。

图1 声波收发器及所形成的声波路径Fig.1 Acoustic transceivers and acoustic paths formed

LSM 法[9]用式（11）重建出N 维声慢度向量S：

式中：L=（Lki）k=1，…，M，i=1，…，N，Lki表示第k（k=1，2，…，M）条声波路径在第i（i=1，2，…，N）个网格内的路径长度，t=（t1，…，tM）T为M 维声波飞行时间向量。

svdDR 法[1，4]基于重建矩阵A的奇异值分解和Tikhonov 正则化获得式（6）的稳定解，如式（12）所示：

式中：uj、vj分别是A的左、右奇异值向量；σ1≥σ2≥…≥σp＞0 是A 的奇异值；p 是非零奇异值的数目；μ（μ＞0）是正则化参数。

用LSM 法重建时，被测区域均匀地划分为5×5=25 个网格。采用svdDR 法和本文提出的pcaIR 法时，被测区域划分成9×9=81 个网格，径向基函数形状参数为0.0001。svdDR 法的正则化参数为1×1010。pcaIR 法的正则化参数为1×109，松弛因子α 为0.9，迭代次数k 为13 次。为更好地表述复杂温度场，本文给出的重建图像均采用三次样条插值，细化后的像素为51×51=2601，即N′＝2601。

对声波飞行时间数据添加标准差为1×10-5的高斯白噪声。三种算法对应的重建误差如表1 所示。模型温度场以及3 种重建算法给出的重建图像如图2 所示。由表1 和图2 可以看出：①三幅重建图像都正确地体现了模型温度场的基本特征；②本文算法获得的重建图像与模型温度场最接近；③本文算法对应的重建误差最小，与LSM 法、svdDR 法相比，本文算法的重建误差最高可降低86.62%和29.1%；④LSM 法、svdDR 法和本文算法的重建时间分别为0.38 ms、0.82 ms、3.87 ms，虽本文算法的重建时间最长，但也只有3.87 ms，即3 种算法都是实时性很好的重建算法。

表1 重建误差比较Tab.1 Comparison of reconstruction errors

图2 模型温度场及重建温度场Fig.2 Model temperature fields and reconstructed temperature fields

5 重建算法的实测数据验证

采用自开发的声学CT 温度场重建系统对本文算法进行了实测数据验证。收发器（由传声器和扬声器组成）布局及形成的54 条有效穿越被测层面的声波路径，如图1 所示。在被测层面下方放置1～2个电加热器，模拟单热点、偏置的单热点和双热点温度场。在每个完整测量周期内，扬声器轮流发声；传声器将接收到的声波转换为电信号；计算机通过同步数据采集卡采集调理后的电信号，用基于互相关的时延估计算法，计算出声波在54 条声波路径上的传播时间，获得实测飞行时间数据。

采用LSM 法、svdDR 法和本文算法对实测数据进行了温度场重建。各算法参数设置同仿真数据验证。单热点、偏置的单热点和双热点温度场的布置照片以及用3 种重建算法获得的重建温度场，如图3 所示。考虑到照片形变，热点位置示意图也在图3中给出。用Fluke 54 II 和K 型热电偶实测的热点温度与各算法重建的热点温度对比，如表2 所示。由图3 和表2 可以看出，与LSM、svdDR 法相比，本文算法能更准确地重建出温度场特征，热点温度也更接近于热电偶的测量值。

表2 热点温度Tab.2 Hot-spot temperature

图3 实测温度场Fig.3 Measured temperature fields

6 结语

本文提出一种基于PCA 降维和迭代正则化的声学CT 温度场重建算法，即pcaIR 法。数值仿真和实验验证表明，与经典的最小二乘法、基于奇异值分解的直接正则化法相比，pcaIR 法的重建误差更低，重建的温度场更接近于真实分布。pcaIR 法虽是迭代算法，但由于计算简单，所需要的迭代次数也很少，所以pcaIR 法的重建时间不足4 ms，远小于获取一帧投影数据的数据采集时间。因此本文所提pcaIR 法有望应用于实际的温度场重建。