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—），男，教授，主要從事無損檢測和數學物理反問題方面的研究.

（責任編校 紀翠榮）