SFTSVD反演方法及其在材料光学参数分布重构中的应用

刘 欣，潘 磊，吕跃凯

（天津师范大学 物理与电子信息学院，天津 300387）

SFTSVD反演方法及其在材料光学参数分布重构中的应用

刘 欣，潘 磊，吕跃凯

（天津师范大学 物理与电子信息学院，天津 300387）

以频域热传导方程为研究目标，通过数学技巧得到模型格林函数的解析解.通过分析模型的本征性质，建立了一种能够根据信号的误差程度自动寻取截断点的新方法，并寻求其实用化的在线反演新途径.

光热检测；奇异值分解；自适应截断奇异值分解方法

近年来，人们不断探索新的检测手段和方法，用以提高检测带宽和数据质量，拓宽检测和分析范围，扩大应用领域，有些研究小组已开发出商品化的扫描热像显微镜和声扫描显微镜［1］.由于光声（Photoacoustic，PA）和光热（Photothermal，PT）技术具有灵敏度高、频响宽和适应性强（适用于气体、液体和固体样品）等优点，已成为一种重要的检测与评价手段，PA／PT检测技术在20世纪80年代已经发展到相当成熟的水平，气体麦克风检测、压电检测、热电检测、光热位移技术、调制光散射技术、热透镜与热光栅技术［2－6］以及近年来出现的红外辐射检测技术等检测手段被广泛应用于材料科学、无损检测、医疗诊断和航天航空等诸多领域［3－7］.同时，借助应用数学与控制理论的研究成果，研究人员已经建立了多种针对具体实际问题的反演算法，如奇异值分解方法、最佳摄动法和迭代正则化方法等.近年来，一些现代优化计算方法，如模拟退火算法、遗传算法、禁忌搜索算法、人工神经网络识别方法和蒙特卡罗算法等也相继被应用于反演的研究.

近20年来，基于PA／PT的各种反演技术一直是反演领域的研究热点之一.在材料某一局部检测到的PA／PT调制信号携带着样品此位置的许多物理信息，利用适当的方法即可在一定范围内对材料的弹性系数、热传导系数、热扩散率、光吸收系数和热源分布等物理量予以重构，继而实现对受检样品的局部成像.目前，针对PA／PT反演技术尚无全面的评述性文章，建立起的反演方法大多局限于数值模拟和实验个例意义上的成功，很难对各种重构技术给出全面客观的评价.反演技术的实用化主要存在3个关键问题：一是反演过程算程较长，如奇异值分解（Sigular Value Decomposition，SVD）方法；二是算法抗噪能力较弱，很多方法对于噪声扰动非常敏感，如果噪声太大就很难从实测信号中解析出预想结果；三是难以实现实用化的在线反演，如截断奇异值分解（Truancted Sigular Value Decomposition，TSVD）反演方法.因此，结合实际问题继续探索新方法或将新的相关数学理论与方法有效应用于反演领域是今后有关物理工作者的主要课题.

本研究以PT检测技术为物理背景，对频域热传导反问题数学模型的本征性质进行深入研究与分析.在此基础上，改进了截断奇异值分解反演方法，建立了一种能够根据信号误差程度自动寻取截断点的新方法，从而寻求在线反演实用化的新途径.

1 理论模型与分析

1.1 实验背景

图1是光热检测实验示意图.在光声光热检测中，周期调制的激光投射于样品表面，样品吸收激光能量后，在样品内激励声波和热波，波在样品异常处引起反射，通过适当的检测手段，可以在样品某局部检测到反射信号.反射信号包含着样品的许多物理信息，通过适当的反演算法可以从检测信号中解析出有效信息，实现对非均匀样品热学和光学等物理参数的分布重构，本研究提取其中的光热（PT）信号.

图1 光声光热检测实验示意图Figure 1 Experimental diagram of photoacoustic and photothermal detection

1.2 频域热传导模型

在光热检测实验中，样品表面PT信号与其表面温度成比例，理论上表面温度可通过求解热传导方程获得.对于热参数均匀但光学参数分布非均匀的样品，通过对时域热传导方程的傅里叶变换［7－9］可分别写出样品上方空气层、样品和衬底3个区域在频域中的热传导方程.

利用式（1a）—式（1c）及其定解条件（式（2））可将样品内热传导问题归结为

1.3 热传导方程反问题模型

根据格林函数理论，可将式（3a）的解写为第一型Fredholm积分

由式（4）可以得到与实测PA信号相关联的样品表面温度的积分表达式

热源函数Q（z）是待解的目标函数；T（ω）为检测信号，可由实验测得.反演的目的就是求得式（8）的反解，以确定目标函数Q（z）.Q（z）解出后，可利用

进一步确定样品光吸收系数的分布函数.为达到此目的，用一组归一化的正交基｛φk，k＝1，2，…，N｝将热源函数展开为

将式（11）带入式（8），可将其离散化为M×N线性方程组

式（12）中，Q为待解的目标向量，当ω变化时T为与检测信号相关的M×1列向量.系数矩阵G的矩阵元为

1.4 SFTSVD方法

形式上，模型方程式（12）的精确解和最小二乘解可分别表示为

式（14）中的H为转置共轭.实际应用中，由于物理条件的制约，对于数据信息量和数据质量2个方面，式（14）难以保证解的适定性和稳定性，即从实际问题中抽象而来的物理模型从本质上具有潜在的病态，从数学的角度说，算子G－1和GH都是不连续算子，因此，上述定义的解在实际计算中不可避免地因发散而无意义.针对上述问题，有2种具有代表性且广泛使用的方法：一种是截断奇异值分解法［7－13］，另一种是正则法（Regularization Method，RM）［10］.TSVD 方 法是对矩阵G（核函数）进行奇异值分解（Sigular value decomposition，SVD）［11］

则可将式（11）的解表达为级数形式：

容易证明，向量Vi是矩阵GHG对应本征值λi的本征向量.由式（17）可以看出，在已知数据T不含任何误差的理想情况下，只要取足够多的级次，以本征向量｛Vi，i＝1，2，…，N｝作为展开基就可以充分逼近目标函数Q.而实际应用中，由于模型的固有病态以及测量数据不可避免地携带着某种程度的噪声等原因，在逼近过程中，只有前若干个低级次项起有效作用，而高级次项不仅无助于提高逼近精度，反而会引起解的剧烈震荡甚至发散.因此，实际应用中可以将起不稳定作用的高阶项截断，从而将式（4）改写为

以上推导即为TSVD的原理，不难看出截断的目的是以精度的损失换取解的稳定.运用TSVD的关键和难点是如何选取截断位置，使解在稳定的同时又具有尽可能高的精度.理论上，截断点位置的选取取决于已知数据（信号）的误差程度，而实际应用中往往难以定量估计信噪比，因此很难做到截断点的准确定位.

从式（17）可以看出，级数是否收敛取决于内积（U＊i，T）和奇异值λi的比值.图2为不同噪声下的内积曲线.由图2可以看出，方程的病态导致奇异值曲线迅速衰减为零，而内积曲线在无噪声扰动的情况下趋于零的速率大于奇异值曲线，然而，在噪声干扰下内积曲线到达一个拐点后开始上扬，而且噪声强度越强，拐点位置越靠前.

图2 不同噪声下的内积曲线Figure 2 Inner product matrix of different noise intensity

图3 不同噪声强度下的系数ln曲线Figure 3 Coefficient matrix lnof different noise intensity

2 数值模拟与讨论

2.1 反演结果

为检验SFTSVD算法的实效，本研究以光热检测信号为已知信息T（ω），以光学性质不均匀材料的热源分布函数Q（z）和光吸收系数分布函数β（z）为重构目标，研究了不同样品的反演效果.操作过程中，先利用SFTSVD算法得到Q（z）的重构数据，再利用式（10）给出β（z），重构结果如图4所示，图中实线代表原始数据，不同符号曲线代表不同的重建数据.图4（a）展示了在信号无误差的理想条件下，对几个单调分布的热源函数Q（z）和光吸收系数β（z）的反演结果，由图4（a）可以看出，重构数据与原始数据重合良好.图4（b）为非单调函数的反演结果，数据良好的符合性说明反演方法对个别复杂函数的适应性.在实际应用中，信号不可避免地带有一定强度的噪声，为了检验算法的实用性，我们在信号中惨入1.00%的噪声，反演效果如图5所示.

图4 在信号无噪声的理想条件下，利用SFTSVD对函数Q（z）和β（z）的重构效果Figure 4 Reconstructed results for Q（z）andβ（z）by SFTSVD under ideal conditions

图5展示了在信号掺有1.00%噪声强度的情况下，热源函数Q（z）和光吸收系数β（z）的重构结果，良好的符合性说明SFTSVD算法具有较好的实用性.从图5中可以看出，热源函数Q（z）的重构效果要优于光吸收系数β（z），这是因为算法直接对热源函数Q（z）进行反演，而光吸收系数β（z）是通过式（10）解出的，而函数Q（z）和β（z）的关系是非线性的，会放大噪声.同时，图5中重构效果的误差会随着深度的增加而变大，这是因为样品对光热信号具有吸收作用，随着深度的增大，反射回的信号强度会有所减弱，所携带的信息也会减少，因此导致误差增大.

为了进一步检验算法对于噪声的自适应性，本研究对不同噪声强度下的反演结果进行比较，结果如图6所示.由图6可以看出，当噪声强度增大到5.00%时，重构数据与原数据仍具有较好的符合，这说明SFTSVD算法有利于实用化反演.

图5 在信号掺有1.00%噪声的条件下，利用SFTSVD对函数Q（z）和β（z）的重构效果Figure 5 Reconstructed results for Q（z）andβ（z）by SFTSVD under 1.00%noise in the signal conditions

图6 在不同噪声强度下，利用SFTSVD对Q（z）与β（z）的反演结果Figure 6 Reconstructed results for Q（z）andβ（z）by SFTSVD under different noise intensity

2.2 SFTSVD与正则方法的比较

正则方法的关键在于寻找正则参数，但寻找适当的正则参数非常困难，正则参数过大会导致精度下降，过小则导致解发散.图7为正则参数α对奇异值倒数发散的抑制效果.

图7 正则参数α对奇异值倒数发散的抑制效果Figure 7 Inhibitory effect of the regularization parameterα on divergence of reciprocal of the singular value

3 结论

通过数学技巧得到正问题格林函数的解析解，在此基础上，通过对模型的分析建立了SFTSVD算法，利用此算法检验了具有不同单调性的物理参数分布样品的重构效果，并进一步检验了算法对噪声的适应性.通过研究不同样品的反演效果，证明SFTSVD方法可以高效准确地对一维情况下的热源函数Q（z）和光吸收系数β（z）进行反演，并且能够抵抗一定强度的噪声且稳定性高，有利于在实际应用中进行实用化反演.

［1］ 李云红，孙晓刚，原桂彬.红外热像仪测温技术［J］.光学精密工程，2007，15（9）：1336－1341.

［2］ Andreas M，Peter H.Progress in Photoacoustic and Photothermal Science and Technology［M］.New York：Elsevier，1997：1－10.

［3］ Venkata S C，Michael S S，Daumantas C.Recent advances in application of acoustic，acousto－optic and photoacoustic methods in biology and medicine［J］.Physica Stautus Solidi，2007，204（10）：3209－3236.

［4］ 张弛，汪源源.声速不均匀介质的光声成像重建算法［J］.光学学报，2008，28（12）：2296－2310.

［5］ Martin M，Majaron B，Joseph S N.Pulsed photothermal temperature profiling of agar tissue phantoms［J］.Lasers in Medical Science，2007，22（4）：279－284.

［6］ Hanh B D，Faubel W，Heissler S，et al.Pharmaceutical applications of photothermal beam deflection［J］.Laser Physics，2006，16（5）：794－798.

［7］ Qu H，Qin H W，Xin X G，et al.Reconstruction of depth profiles of thermal conductivity of case hardened steels using a three－dimensional photothermal technique［J］.Journal of Applied Physics，2008，104（11）：113518－113527.

［8］ Yotsna R，Lu Y，Longuemart S，et al.Optothermal depth profiling by neural network infrared radiometry signal recognition［J］.Appl Phys，2005，97（1）：4701－4706.

［9］ Voti R L，Sibilia C，Bertolotti M.Photothermal depth profiling by thermal wave backscattring and genetic algorithms［J］.International Journal of Thermophysics，2005，26（6）：1833－1848.

［10］ Wang W，Ye C H，Li P.Multi－exponential inversions of nuclear magnetic resonance relaxation signal［J］.Science in China，2001，44（11）：1477－1484.

［11］ Gérard K，George K，Erwan L P，et al.Dominique picard inversion of noisy radon transform by SVD based needlets［J］.Applied and Computational Harmonic Analysis，2010，28（1）：24－25.

［12］ Barriere P，Idier J，Goussard Y，et al.Fast solutions of the 2D inverse scattering problem based on a TSVD approximation of the internal field for the forward model［J］.Antennas and Propagation，2010，58（12）：4015－4024.

［13］ Zhi L.The canonical TSVD method for L－generalized solution of linear ill－posed equations［J］.Advanced Materials Research，2010，211（212）：177－181.

SFTSVD inversion method and its applications on reconstruction of optical parameters of materials

LIUXin，PANLei，LÜYuekai
（College of Physics and Electronic Information Science，Tianjin Normal University，Tianjin 300387，China）

Heat conduction equation of frequency domain is the research target.The analytical solution of Green's function of the model is obtained by mathematical techniques.Through the analysis of intrinsic properties of inverse model，a new way which can choose cut－off point automatically according to the error of signal is established to seek new ways of practical online inversion.

photo－thermal detection；singular value decomposition；self－fitted truncated singular decomposition

O175.1

A

1671－1114（2011）04－0027－06

投稿日期：2011－05－12

刘 欣（1986—），男，硕士研究生.

吕跃凯（1958—），男，教授，主要从事无损检测和数学物理反问题方面的研究.

（责任编校 纪翠荣）