劉 欣，潘 磊，呂躍凱

（天津師范大學(xué) 物理與電子信息學(xué)院，天津 300387）

SFTSVD反演方法及其在材料光學(xué)參數(shù)分布重構(gòu)中的應(yīng)用

劉 欣，潘 磊，呂躍凱

（天津師范大學(xué) 物理與電子信息學(xué)院，天津 300387）

以頻域熱傳導(dǎo)方程為研究目標(biāo)，通過(guò)數(shù)學(xué)技巧得到模型格林函數(shù)的解析解.通過(guò)分析模型的本征性質(zhì)，建立了一種能夠根據(jù)信號(hào)的誤差程度自動(dòng)尋取截?cái)帱c(diǎn)的新方法，并尋求其實(shí)用化的在線反演新途徑.

光熱檢測(cè)；奇異值分解；自適應(yīng)截?cái)嗥娈愔捣纸夥椒?/p>

近年來(lái)，人們不斷探索新的檢測(cè)手段和方法，用以提高檢測(cè)帶寬和數(shù)據(jù)質(zhì)量，拓寬檢測(cè)和分析范圍，擴(kuò)大應(yīng)用領(lǐng)域，有些研究小組已開(kāi)發(fā)出商品化的掃描熱像顯微鏡和聲掃描顯微鏡［1］.由于光聲（Photoacoustic，PA）和光熱（Photothermal，PT）技術(shù)具有靈敏度高、頻響寬和適應(yīng)性強(qiáng)（適用于氣體、液體和固體樣品）等優(yōu)點(diǎn)，已成為一種重要的檢測(cè)與評(píng)價(jià)手段，PA／PT檢測(cè)技術(shù)在20世紀(jì)80年代已經(jīng)發(fā)展到相當(dāng)成熟的水平，氣體麥克風(fēng)檢測(cè)、壓電檢測(cè)、熱電檢測(cè)、光熱位移技術(shù)、調(diào)制光散射技術(shù)、熱透鏡與熱光柵技術(shù)［2－6］以及近年來(lái)出現(xiàn)的紅外輻射檢測(cè)技術(shù)等檢測(cè)手段被廣泛應(yīng)用于材料科學(xué)、無(wú)損檢測(cè)、醫(yī)療診斷和航天航空等諸多領(lǐng)域［3－7］.同時(shí)，借助應(yīng)用數(shù)學(xué)與控制理論的研究成果，研究人員已經(jīng)建立了多種針對(duì)具體實(shí)際問(wèn)題的反演算法，如奇異值分解方法、最佳攝動(dòng)法和迭代正則化方法等.近年來(lái)，一些現(xiàn)代優(yōu)化計(jì)算方法，如模擬退火算法、遺傳算法、禁忌搜索算法、人工神經(jīng)網(wǎng)絡(luò)識(shí)別方法和蒙特卡羅算法等也相繼被應(yīng)用于反演的研究.

近20年來(lái)，基于PA／PT的各種反演技術(shù)一直是反演領(lǐng)域的研究熱點(diǎn)之一.在材料某一局部檢測(cè)到的PA／PT調(diào)制信號(hào)攜帶著樣品此位置的許多物理信息，利用適當(dāng)?shù)姆椒纯稍谝欢ǚ秶鷥?nèi)對(duì)材料的彈性系數(shù)、熱傳導(dǎo)系數(shù)、熱擴(kuò)散率、光吸收系數(shù)和熱源分布等物理量予以重構(gòu)，繼而實(shí)現(xiàn)對(duì)受檢樣品的局部成像.目前，針對(duì)PA／PT反演技術(shù)尚無(wú)全面的評(píng)述性文章，建立起的反演方法大多局限于數(shù)值模擬和實(shí)驗(yàn)個(gè)例意義上的成功，很難對(duì)各種重構(gòu)技術(shù)給出全面客觀的評(píng)價(jià).反演技術(shù)的實(shí)用化主要存在3個(gè)關(guān)鍵問(wèn)題：一是反演過(guò)程算程較長(zhǎng)，如奇異值分解（Sigular Value Decomposition，SVD）方法；二是算法抗噪能力較弱，很多方法對(duì)于噪聲擾動(dòng)非常敏感，如果噪聲太大就很難從實(shí)測(cè)信號(hào)中解析出預(yù)想結(jié)果；三是難以實(shí)現(xiàn)實(shí)用化的在線反演，如截?cái)嗥娈愔捣纸猓═ruancted Sigular Value Decomposition，TSVD）反演方法.因此，結(jié)合實(shí)際問(wèn)題繼續(xù)探索新方法或?qū)⑿碌南嚓P(guān)數(shù)學(xué)理論與方法有效應(yīng)用于反演領(lǐng)域是今后有關(guān)物理工作者的主要課題.

本研究以PT檢測(cè)技術(shù)為物理背景，對(duì)頻域熱傳導(dǎo)反問(wèn)題數(shù)學(xué)模型的本征性質(zhì)進(jìn)行深入研究與分析.在此基礎(chǔ)上，改進(jìn)了截?cái)嗥娈愔捣纸夥囱莘椒ǎ⒘艘环N能夠根據(jù)信號(hào)誤差程度自動(dòng)尋取截?cái)帱c(diǎn)的新方法，從而尋求在線反演實(shí)用化的新途徑.

1 理論模型與分析

1.1 實(shí)驗(yàn)背景

圖1是光熱檢測(cè)實(shí)驗(yàn)示意圖.在光聲光熱檢測(cè)中，周期調(diào)制的激光投射于樣品表面，樣品吸收激光能量后，在樣品內(nèi)激勵(lì)聲波和熱波，波在樣品異常處引起反射，通過(guò)適當(dāng)?shù)臋z測(cè)手段，可以在樣品某局部檢測(cè)到反射信號(hào).反射信號(hào)包含著樣品的許多物理信息，通過(guò)適當(dāng)?shù)姆囱菟惴梢詮臋z測(cè)信號(hào)中解析出有效信息，實(shí)現(xiàn)對(duì)非均勻樣品熱學(xué)和光學(xué)等物理參數(shù)的分布重構(gòu)，本研究提取其中的光熱（PT）信號(hào).

圖1 光聲光熱檢測(cè)實(shí)驗(yàn)示意圖Figure 1 Experimental diagram of photoacoustic and photothermal detection

1.2 頻域熱傳導(dǎo)模型

在光熱檢測(cè)實(shí)驗(yàn)中，樣品表面PT信號(hào)與其表面溫度成比例，理論上表面溫度可通過(guò)求解熱傳導(dǎo)方程獲得.對(duì)于熱參數(shù)均勻但光學(xué)參數(shù)分布非均勻的樣品，通過(guò)對(duì)時(shí)域熱傳導(dǎo)方程的傅里葉變換［7－9］可分別寫(xiě)出樣品上方空氣層、樣品和襯底3個(gè)區(qū)域在頻域中的熱傳導(dǎo)方程.

利用式（1a）—式（1c）及其定解條件（式（2））可將樣品內(nèi)熱傳導(dǎo)問(wèn)題歸結(jié)為

1.3 熱傳導(dǎo)方程反問(wèn)題模型

根據(jù)格林函數(shù)理論，可將式（3a）的解寫(xiě)為第一型Fredholm積分

由式（4）可以得到與實(shí)測(cè)PA信號(hào)相關(guān)聯(lián)的樣品表面溫度的積分表達(dá)式

熱源函數(shù)Q（z）是待解的目標(biāo)函數(shù)；T（ω）為檢測(cè)信號(hào)，可由實(shí)驗(yàn)測(cè)得.反演的目的就是求得式（8）的反解，以確定目標(biāo)函數(shù)Q（z）.Q（z）解出后，可利用

進(jìn)一步確定樣品光吸收系數(shù)的分布函數(shù).為達(dá)到此目的，用一組歸一化的正交基｛φk，k＝1，2，…，N｝將熱源函數(shù)展開(kāi)為

將式（11）帶入式（8），可將其離散化為M×N線性方程組

式（12）中，Q為待解的目標(biāo)向量，當(dāng)ω變化時(shí)T為與檢測(cè)信號(hào)相關(guān)的M×1列向量.系數(shù)矩陣G的矩陣元為

1.4 SFTSVD方法

形式上，模型方程式（12）的精確解和最小二乘解可分別表示為

式（14）中的H為轉(zhuǎn)置共軛.實(shí)際應(yīng)用中，由于物理?xiàng)l件的制約，對(duì)于數(shù)據(jù)信息量和數(shù)據(jù)質(zhì)量2個(gè)方面，式（14）難以保證解的適定性和穩(wěn)定性，即從實(shí)際問(wèn)題中抽象而來(lái)的物理模型從本質(zhì)上具有潛在的病態(tài)，從數(shù)學(xué)的角度說(shuō)，算子G－1和GH都是不連續(xù)算子，因此，上述定義的解在實(shí)際計(jì)算中不可避免地因發(fā)散而無(wú)意義.針對(duì)上述問(wèn)題，有2種具有代表性且廣泛使用的方法：一種是截?cái)嗥娈愔捣纸夥ǎ?－13］，另一種是正則法（Regularization Method，RM）［10］.TSVD 方 法是對(duì)矩陣G（核函數(shù)）進(jìn)行奇異值分解（Sigular value decomposition，SVD）［11］

則可將式（11）的解表達(dá)為級(jí)數(shù)形式：

容易證明，向量Vi是矩陣GHG對(duì)應(yīng)本征值λi的本征向量.由式（17）可以看出，在已知數(shù)據(jù)T不含任何誤差的理想情況下，只要取足夠多的級(jí)次，以本征向量｛Vi，i＝1，2，…，N｝作為展開(kāi)基就可以充分逼近目標(biāo)函數(shù)Q.而實(shí)際應(yīng)用中，由于模型的固有病態(tài)以及測(cè)量數(shù)據(jù)不可避免地?cái)y帶著某種程度的噪聲等原因，在逼近過(guò)程中，只有前若干個(gè)低級(jí)次項(xiàng)起有效作用，而高級(jí)次項(xiàng)不僅無(wú)助于提高逼近精度，反而會(huì)引起解的劇烈震蕩甚至發(fā)散.因此，實(shí)際應(yīng)用中可以將起不穩(wěn)定作用的高階項(xiàng)截?cái)?，從而將式?）改寫(xiě)為

以上推導(dǎo)即為T(mén)SVD的原理，不難看出截?cái)嗟哪康氖且跃鹊膿p失換取解的穩(wěn)定.運(yùn)用TSVD的關(guān)鍵和難點(diǎn)是如何選取截?cái)辔恢?，使解在穩(wěn)定的同時(shí)又具有盡可能高的精度.理論上，截?cái)帱c(diǎn)位置的選取取決于已知數(shù)據(jù)（信號(hào)）的誤差程度，而實(shí)際應(yīng)用中往往難以定量估計(jì)信噪比，因此很難做到截?cái)帱c(diǎn)的準(zhǔn)確定位.

從式（17）可以看出，級(jí)數(shù)是否收斂取決于內(nèi)積（U＊i，T）和奇異值λi的比值.圖2為不同噪聲下的內(nèi)積曲線.由圖2可以看出，方程的病態(tài)導(dǎo)致奇異值曲線迅速衰減為零，而內(nèi)積曲線在無(wú)噪聲擾動(dòng)的情況下趨于零的速率大于奇異值曲線，然而，在噪聲干擾下內(nèi)積曲線到達(dá)一個(gè)拐點(diǎn)后開(kāi)始上揚(yáng)，而且噪聲強(qiáng)度越強(qiáng)，拐點(diǎn)位置越靠前.

圖2 不同噪聲下的內(nèi)積曲線Figure 2 Inner product matrix of different noise intensity

圖3 不同噪聲強(qiáng)度下的系數(shù)ln曲線Figure 3 Coefficient matrix lnof different noise intensity

2 數(shù)值模擬與討論

2.1 反演結(jié)果

為檢驗(yàn)SFTSVD算法的實(shí)效，本研究以光熱檢測(cè)信號(hào)為已知信息T（ω），以光學(xué)性質(zhì)不均勻材料的熱源分布函數(shù)Q（z）和光吸收系數(shù)分布函數(shù)β（z）為重構(gòu)目標(biāo)，研究了不同樣品的反演效果.操作過(guò)程中，先利用SFTSVD算法得到Q（z）的重構(gòu)數(shù)據(jù)，再利用式（10）給出β（z），重構(gòu)結(jié)果如圖4所示，圖中實(shí)線代表原始數(shù)據(jù)，不同符號(hào)曲線代表不同的重建數(shù)據(jù).圖4（a）展示了在信號(hào)無(wú)誤差的理想條件下，對(duì)幾個(gè)單調(diào)分布的熱源函數(shù)Q（z）和光吸收系數(shù)β（z）的反演結(jié)果，由圖4（a）可以看出，重構(gòu)數(shù)據(jù)與原始數(shù)據(jù)重合良好.圖4（b）為非單調(diào)函數(shù)的反演結(jié)果，數(shù)據(jù)良好的符合性說(shuō)明反演方法對(duì)個(gè)別復(fù)雜函數(shù)的適應(yīng)性.在實(shí)際應(yīng)用中，信號(hào)不可避免地帶有一定強(qiáng)度的噪聲，為了檢驗(yàn)算法的實(shí)用性，我們?cè)谛盘?hào)中慘入1.00%的噪聲，反演效果如圖5所示.

圖4 在信號(hào)無(wú)噪聲的理想條件下，利用SFTSVD對(duì)函數(shù)Q（z）和β（z）的重構(gòu)效果Figure 4 Reconstructed results for Q（z）andβ（z）by SFTSVD under ideal conditions

圖5展示了在信號(hào)摻有1.00%噪聲強(qiáng)度的情況下，熱源函數(shù)Q（z）和光吸收系數(shù)β（z）的重構(gòu)結(jié)果，良好的符合性說(shuō)明SFTSVD算法具有較好的實(shí)用性.從圖5中可以看出，熱源函數(shù)Q（z）的重構(gòu)效果要優(yōu)于光吸收系數(shù)β（z），這是因?yàn)樗惴ㄖ苯訉?duì)熱源函數(shù)Q（z）進(jìn)行反演，而光吸收系數(shù)β（z）是通過(guò)式（10）解出的，而函數(shù)Q（z）和β（z）的關(guān)系是非線性的，會(huì)放大噪聲.同時(shí)，圖5中重構(gòu)效果的誤差會(huì)隨著深度的增加而變大，這是因?yàn)闃悠穼?duì)光熱信號(hào)具有吸收作用，隨著深度的增大，反射回的信號(hào)強(qiáng)度會(huì)有所減弱，所攜帶的信息也會(huì)減少，因此導(dǎo)致誤差增大.

為了進(jìn)一步檢驗(yàn)算法對(duì)于噪聲的自適應(yīng)性，本研究對(duì)不同噪聲強(qiáng)度下的反演結(jié)果進(jìn)行比較，結(jié)果如圖6所示.由圖6可以看出，當(dāng)噪聲強(qiáng)度增大到5.00%時(shí)，重構(gòu)數(shù)據(jù)與原數(shù)據(jù)仍具有較好的符合，這說(shuō)明SFTSVD算法有利于實(shí)用化反演.

圖5 在信號(hào)摻有1.00%噪聲的條件下，利用SFTSVD對(duì)函數(shù)Q（z）和β（z）的重構(gòu)效果Figure 5 Reconstructed results for Q（z）andβ（z）by SFTSVD under 1.00%noise in the signal conditions

圖6 在不同噪聲強(qiáng)度下，利用SFTSVD對(duì)Q（z）與β（z）的反演結(jié)果Figure 6 Reconstructed results for Q（z）andβ（z）by SFTSVD under different noise intensity

2.2 SFTSVD與正則方法的比較

正則方法的關(guān)鍵在于尋找正則參數(shù)，但尋找適當(dāng)?shù)恼齽t參數(shù)非常困難，正則參數(shù)過(guò)大會(huì)導(dǎo)致精度下降，過(guò)小則導(dǎo)致解發(fā)散.圖7為正則參數(shù)α對(duì)奇異值倒數(shù)發(fā)散的抑制效果.

圖7 正則參數(shù)α對(duì)奇異值倒數(shù)發(fā)散的抑制效果Figure 7 Inhibitory effect of the regularization parameterα on divergence of reciprocal of the singular value

3 結(jié)論

通過(guò)數(shù)學(xué)技巧得到正問(wèn)題格林函數(shù)的解析解，在此基礎(chǔ)上，通過(guò)對(duì)模型的分析建立了SFTSVD算法，利用此算法檢驗(yàn)了具有不同單調(diào)性的物理參數(shù)分布樣品的重構(gòu)效果，并進(jìn)一步檢驗(yàn)了算法對(duì)噪聲的適應(yīng)性.通過(guò)研究不同樣品的反演效果，證明SFTSVD方法可以高效準(zhǔn)確地對(duì)一維情況下的熱源函數(shù)Q（z）和光吸收系數(shù)β（z）進(jìn)行反演，并且能夠抵抗一定強(qiáng)度的噪聲且穩(wěn)定性高，有利于在實(shí)際應(yīng)用中進(jìn)行實(shí)用化反演.

［1］ 李云紅，孫曉剛，原桂彬.紅外熱像儀測(cè)溫技術(shù)［J］.光學(xué)精密工程，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］ 張弛，汪源源.聲速不均勻介質(zhì)的光聲成像重建算法［J］.光學(xué)學(xué)報(bào)，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，F(xiàn)aubel 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—），男，教授，主要從事無(wú)損檢測(cè)和數(shù)學(xué)物理反問(wèn)題方面的研究.

（責(zé)任編校 紀(jì)翠榮）