基于磁聲耦合效應的生物組織電特性功能成像方法在腫瘤的早期診斷方面具有重要研究意義,其中磁聲信號時頻特性是該檢測成像方法信號檢測處理和圖像重建中的重要研究內容。本文基于磁聲耦合格林函數解,采用基于時間延遲的波形求和方法,對矩形介質模型的準一維下聲傳播情況進行分析,計算磁聲耦合聲信號的時域波形,對不同厚度介質模型的磁聲信號進行正問題仿真計算,并通過實驗測量磁聲信號波形驗證仿真結果,證明了波形的疊加求和使不同電導率分布介質的磁聲信號具有不同的時頻特性。仿真結果表明,當介質厚度小于波形周期時,磁聲信號寬度大于激勵脈寬;隨著介質厚度的增加,磁聲信號脈沖寬度逐漸增加,峰值頻率逐漸降低;當介質厚度大于激勵周期時,磁聲信號脈沖頻率特性與激勵頻率近似相同。不同厚度介質樣本的實驗驗證了仿真結果。該研究為磁聲耦合成像中基于磁聲信號時頻特性對介質聲源電導率的重建打下基礎。
引用本文： 張順起, 殷濤, 馬任, 劉志朋. 基于波形求和的不同厚度介質磁聲信號時頻特性研究. 生物醫學工程學雜志, 2015, 32(4): 717-724. doi: 10.7507/1001-5515.20150131 復制
引言
生物組織電特性反映了組織的生理病理特征[1]，研究表明，正常組織與腫瘤病變組織電導率具有明顯差異[2-4]，因此組織電特性的檢測和成像將為相關疾病的的早期診斷提供重要依據。磁聲耦合效應的生物組織電特性功能成像方法[5-9]，由于具有無創性，造價較低，能夠反映生物組織早期病變，同時兼具電阻抗成像高對比度和超聲成像高空間分辨率的優點，成為了近些年研究的熱點。
目前國內外有多篇文獻對磁聲成像方法進行研究報道[10-19]，由于該方法涉及電、磁、聲、信號與系統等多學科交叉，目前還處于實驗室研究階段，相關機制以及檢測成像技術研究以仿體樣本為主。其中磁聲信號正問題求解采用有限元方法(finite elements method,FEM)，設置激勵為沖激函數[11]，則正問題磁聲信號脈沖即對應空間電導率的邊界。在磁聲信號處理和逆問題電導率重建過程中[12-13]，信號經放大濾波和數字平均處理去除隨機噪聲[14]后，進行聲源分布重建，然后在已知介質為實心體的情況下，根據介質邊界聲源的幅值推導介質的相對電導率圖像。然而實際測量中，由于儀器設備帶寬限制，難以實現無限窄脈沖的激勵，使得實測磁聲信號脈沖特性復雜。因此成像分辨率也受到影響，即受限于激勵脈沖寬度，無法分辨小于脈沖寬度的成像目標，如比較薄的結締組織、隔膜、血管壁等，其前后邊界聲源發出的聲信號發生混疊，即兩個或多個邊界聲源脈沖信號疊加，使得逆問題研究中，在未知介質電導率分布的情況下，無法分辨介質為實心電導組織，還是薄層膜狀組織，從而影響圖像重建的準確性。
因此，本文利用波形求和的正問題求解方法，對不同厚度的電導率介質的磁聲信號時頻特性進行研究，從磁聲信號中提取并分析比較薄層介質的電導率聲源和厚層介質的邊界聲源信息的異同點，為在未知介質電導率分布的情況下由磁聲信號分辨和重建內部聲源電導率圖像提供依據。
1 原理
1.1 磁聲耦合原理
磁聲成像方法通過對介質施加電磁激勵，介質內部質點受到洛倫茲力作用發生振動，形成聲波，在介質外部可由聲換能器接收聲信號響應，如圖 1。磁聲成像滿足波動方程[15]：
| ${\nabla ^2}p\left({\boldsymbol{r},t} \right)- \frac{1}{{c_x^2}}\frac{{\partial p\left({\boldsymbol{r},t} \right)}}{{\partial {t^2}}}=\nabla \left({\boldsymbol{r},t} \right) \cdot \left[{\boldsymbol{J}\left({\boldsymbol{r},t} \right) \times {\boldsymbol{B}_0}\left(\boldsymbol{r} \right)} \right]$ |
其中p(r,t)為r處聲信號，J為電流密度,B0為靜磁場的磁感應強度，R為r到傳感器的距離,cs介質內部聲傳播速度。
由格林函數求解波動方程，介質外部聲信號響應[15]
| $\begin{array}{l} p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left[{\sigma \left({\boldsymbol{r}'} \right.\boldsymbol{E}\left({\boldsymbol{r}'} \right) \times {\boldsymbol{B}_0}\left({\boldsymbol{r}'} \right)} \right]} \\ \frac{{\delta \left({t- R/{c_s}} \right)}}{R} \end{array}$ |
其中▽r′·[σ(r′)E(r′)×B0(r′)]為聲源，σ為介質電導率，E為介質內部電場，洛倫茲力為f=σ(r′)E(r′)×B0(r′)。其中延遲項包含了聲振動源到傳感器位置r的傳播距離引起的傳播延遲[16-17]。
因此，r處檢測的磁聲信號p(r,t)可表示為介質內部聲源在不同傳輸延遲下的積分。
圖1
磁聲耦合裝置原理圖
Figure1.
Scheme of magneto-acoustic effect experimental set-up
1.2 波形求和算法
由前所述，電流密度J可分離為空間項和時間項的乘積[15]，為了計算磁聲信號的時域波形，我們將計算過程分為2步：空間項求解和時間項求解，即計算介質內電流密度的空間分布，再對各剖分單元分別考慮其隨時間變化的信號波形，根據各剖分單元電流密度模和方向，利用磁聲信號格林函數解，求解磁聲信號。
考慮簡化的矩形介質模型，設電流沿y方向，靜磁場沿－z方向，則洛倫茲力方向沿x方向，這時聲信號可視為沿一維方向(x軸)傳播，如圖 1所示。設檢測器端面與樣本前端面距離為d，樣本長度l，設樣本尺寸位于0≤x≤l區間內，則接收到的磁聲信號為整個介質各剖分單元聲信號求和的結果。
| $p=\sum\limits_{i=1}^n {{p_i}} $ |
其中pi為第i個剖分單元的聲壓。
由于聲源項▽r′·(J×B0)中電流密度J可以分離為空間量和時間量乘積的形式[15],則聲源表達式(2)可變形為
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left({\boldsymbol{J} \times {\boldsymbol{B}_0}} \right)\frac{{\delta \left({t- R/{c_s}} \right)}}{R}} \otimes s\left(t \right)$ |
其中s(t)為電流激勵信號,應用卷積性質和沖激函數的抽樣性質，式(4)可寫為
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left({\boldsymbol{J} \times {\boldsymbol{B}_0}} \right)\frac{{s\left({t- R/{c_s}} \right)}}{R}} $ |
將積分過程離散為求和形式
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\sum\limits_i {\frac{{\rm{d}}}{{{\rm{d}}y}}} \left({{\boldsymbol{J}_i}{\boldsymbol{B}_0}} \right)\cos {\theta _i}\frac{{s\left({t- dela{y_i}} \right)}}{{dela{y_i} \cdot {c_s}}}$ |
其中θi為剖分單元i受到洛倫茲力方向與傳感器軸線方向的夾角，式(6)可表示為以下形式
| $p\left({\boldsymbol{r},t} \right)=\frac{{{\boldsymbol{B}_0}}}{{4\pi }}\sum\limits_{i=1}^n {\frac{{\rm{d}}}{{{\rm{d}}y}}} {\boldsymbol{J}_i}\frac{{s\left({t- dela{y_i}} \right)}}{{dela{y_i} \cdot {c_s}}}\cos {\theta _i}$ |
應用式(7),則可計算出由介質全部剖分單元形成的磁聲信號，通過設置不同的介質厚度l，則可獲得其對應的磁聲信號波形。
2 方法
2.1 仿真方法
仿真研究中設電流方向沿y軸，靜磁場沿z軸，洛倫茲力即聲振動方向沿x軸，如圖 1所示。本研究利用有限元工具Comsol Multiphysics，將介質聲源剖分為若干單元，計算空間電流密度，再將結果導出，并由MATLAB計算，獲得檢測位置r處的磁聲信號，仿真計算過程如圖 2。
考慮矩形介質的準一維的聲傳播情況，本研究中選擇樣本厚度分別為l<λ/4、l=λ/4、λ/4<l<λ/2、l=λ/2、λ/2<l<λ、l=λ以及l>λ幾種情況進行仿真：仿真中樣本厚度分別設為0.6、0.75、1.125、1.5、2.25、3和20 mm；電導率均為1 S/m，聲學參數均勻；通過電極對樣本施加電流激勵，激勵為一個周期正弦信號。
| $\boldsymbol{J}\left(t \right)=A\sin 2\pi ft{a_x}\left({0 \le 2\pi ft \le 2\pi } \right)$ |
其中A為信號幅值，頻率f=500 kHz。
聲檢測器與樣本前端面距離為d，檢測器中心頻率為500 kHz；圖 2中顯示了計算步驟。有限元剖分網格尺寸小于0.12 mm。
圖2
磁聲信號正問題計算過程
Figure2.
Computing procedure of magneto-acoustic signal
2.2 實驗方法
為驗證不同厚度磁聲信號時頻特性，本研究采用磁聲成像實驗系統[19-20]，和明膠仿體樣本進行實驗。裝置如圖 1，由函數發生器(AFG3252，美國泰克)產生激勵信號，通過雙極性放大器(HSA4104，日本NF)放大，將激勵信號加載到實驗樣本上。由電磁鐵產生恒穩磁場，磁場在半徑50 mm、高度100 mm的圓柱形空間范圍內均勻。由超聲傳感器(V303，美國Olympus，中心頻率500 kHz)接收聲信號。超聲傳感器接收的聲信號經由放大器(5307，NF)100倍放大后，由示波器(MSO4000，美國泰克)進行波形顯示。采用單周期正弦脈沖信號為激勵，頻率500 kHz，靜磁場0.4 T。檢測距離為60 mm，檢測波形平均512次，以去除隨機噪聲。樣本為長方體，采用明膠制作，混入5% NaCl，樣本截面為10mm×5mm，厚度分別為l=0.6 mm和l=22.5 mm。
3 結果
3.1 仿真結果
介質內部電流密度仿真結果如圖 3所示，其中的樣本厚度l分別為0.6、0.75、1.125、1.5、2.25、3和20 mm，寬度均為10 mm，各樣本電流沿y方向均勻分布。其對應檢測位置的波形和短時傅里葉變換頻譜如圖 4所示。
圖3
介質內電流密度(y方向)有限元仿真結果
樣本寬度均為10 mm，由圖可見，各樣本電流沿
the width of each sample is 10 mm. the result shows that the current in the samples with different thickness distributes along
為進一步比較分析不同厚度介質磁聲信號時頻特性，對各厚度樣本的脈沖半寬度、峰值頻率進行了比較，如表 1所示。
表1
不同樣本厚度的時頻特性
Table1.
Time-frequency characteristics of signals from samples with different thickness
時域上，由圖 4(a)(c)(e)(g)(i)(k)(m)可見，激勵信號為單周期正弦脈沖時，對于厚度l小于聲信號波長λ的薄層樣本，當厚度小于1/4聲信號波長時兩邊界對應脈沖信號發生混疊，形成兩個正峰一個負峰的波形，且隨著樣本厚度增加，脈沖逐漸展寬。對于厚度l大于聲信號波長λ的厚層樣本，當厚度大于波長2倍時，信號波形中兩脈沖分離，兩脈沖極性相反，說明電導率變化趨勢相反。由表 1第4列可見，激勵頻率500 kHz，其脈沖波長λ=v/f=1 500 ms-1/500 kHz=3 mm，當樣本厚度分別為0.6、0.75、1.125、1.5、2.25、3 mm時，對應脈沖的半寬度(大于波形峰值1/的寬度)，分別由激勵信號的0.66 μs變為0.67、0.68、0.73、0.83、0.99和0.99 μs。
頻域上，由短時傅里葉變換，即圖 4(b)(d)(f)(h)(j)(l)(n)和表 1最右列可見，各厚度的薄層樣本聲信號脈沖對應的峰值頻率由423 kHz變為406、397、369、336、273和223 kHz，說明薄層樣本峰值頻率隨樣本厚度增加逐漸降低；而厚層樣本(l=20 mm)的脈沖峰值頻率與激勵信號近似相同，均為423 kHz。
圖4
介質磁聲信號波形及頻譜仿真結果
(a)為樣本厚度
(a) is waveform of sample with thickness
仿真結果表明，樣本厚度小于激勵波長時，隨著厚度增加，脈沖展寬，峰值頻率對應下降；樣本厚度大于波長時，樣本聲信號與激勵波形和頻譜近似相同。此結果表明磁聲信號時頻特性中包含了介質電導率分布信息，尤其對于厚度大于波長的厚層樣本，與薄層樣本的時頻特性具有顯著區別。
3.2 實驗結果
實驗測量的薄層樣本(0.6 mm厚)和厚層樣本(22.5 mm厚)的磁聲信號及短時傅里葉變換頻譜如圖 5、圖 6所示。
圖5
實測樣本的磁聲信號及頻譜
(a)磁聲信號波形；(b)磁聲信號頻譜。樣本厚度
(a) measured waveform; (b) short-time Fourier transform spectrum of the waveform. The rectangular sample’s thickness
圖6
實測樣本的磁聲信號及頻譜
(a)磁聲信號波形；(b)(c)分別為圖中樣本前后邊界對應的脈沖w1、w2的短時傅里葉變換頻譜。樣本厚度
(a) measured waveform; (b) the short-time Fourier transform spectrum of the pulse w1 in (a); (c) the short-time Fourier transform spectrum of the pulse w2 in (a). The rectangular sample’s thickness
由圖 5、圖 6可見，對于相同的激勵信號，在相同的檢測距離下，厚層樣本和薄層樣本的磁聲信號具有顯著不同，對于薄層樣本，其波形呈現W型，其主要原因在于正弦波激勵信號經過聲換能器特性以及成像系統各儀器系統函數的卷積后形成[16-17]；而厚層介質呈現具有正負峰的正弦脈沖波形，如圖 5(a)，其脈沖波形展寬為原來的2倍左右(1.29 μs)，該現象與理論分析一致。
此外，比較圖 4(m)與圖 5(a)可見，仿真結果與實驗具有一定差異，這是由于實驗中介質內部的反射散射以及噪聲的影響，檢測通路的頻率響應限制，使得波形產生一定失真[21]。為此，對實驗數據進一步分析處理。首先對實驗信號進行小波去噪，采用1維dmey小波5層分解，去除雜波成分，并將薄層樣本的測量信號作為一個剖分單元，則一個厚層介質可視為若干個薄層樣本的組合，因此利用一定數量薄層介質的求和，即可獲得對應樣本厚度的厚層介質的磁聲信號。由于薄層樣本的實驗測量信號中包含了成像系統響應函數的卷積，根據卷積的線性性質，采用此處理方法求和計算的厚層樣本磁聲信號同樣包含了成像系統響應函數，因此去掉系統響應函數影響，便于對薄、厚樣本的時頻特性進行分析驗證。結果如圖 7所示。由實驗測量的薄層介質聲信號計算的厚層介質聲信號如圖 7(a)，對信號波形左右兩個脈沖w4、w5進行短時傅里葉變換如圖 7(b)(c)，計算波形短時的峰值頻率可知，薄層仿體信號峰值頻率為459 kHz。而由薄層疊加計算的左右脈沖峰值頻率分別為305 kHz和386 kHz，說明由于波形疊加導致信號的峰值頻率向低頻移動，與仿真結果一致。而對于實測厚層介質峰值頻率為310 kHz和393 kHz，其頻率誤差為1.6%和1.8%，說明疊加處理方法計算的結果與實測結果近似，驗證了波形求和求解及薄厚層時頻特性分析方法的可行性。
圖7
由薄層樣本磁聲信號計算得到的厚層樣本的磁聲信號
(a)為模擬計算獲得的厚層樣本磁聲信號波形；(b)(c)為w4、w5兩個信號的短時傅里葉變換
Figure7. Magneto-acoustic signal of thick sample calculated from thin sample(a) the calculated signal using the measured pulse in Fig. 5; (b) the short-time Fourier transform spectrum of the pulse w4 in (a); (c) the short-time Fourier transform spectrum of the pulse w5 in (a)
為進一步分析時域波形的相似程度，采用相關系數進行計算。相關系數通常用來表征波形的相似程度[22]。本研究選擇實驗波形的前后邊界對應的兩個脈沖46.61~51.76 μs(令其為w1)，60.96~65.75 μs(令其為w2)，同時選擇兩個脈沖中間部分的55.43~59.53 μs(令其為w3)的干擾波形，如圖 6(a)所示，與處理后波形中的脈沖，即圖 7(a)中的w4、w5，分別進行相關系數計算。結果顯示，實驗測量的仿體前后邊界對應波形w1、w2與由薄層波形疊加計算w4、w5的相關系數，分別為97.96%和76.04%，可見左側波形w1相似度較高，而右側仿體邊界w2相似度降低，這是由于聲信號在介質中傳播同時與散射等干擾的疊加結果，導致波形產生變化。而對于干擾脈沖，相關系數計算結果表明，介于仿體兩個邊界中間的內部干擾信號w3，其峰值頻率131 kHz，與計算得到的仿體左右側峰w1、w2波形的相關系數為29.88%和23.85%。可見，對于非介質邊界產生的聲信號，其時頻特性與介質邊界具有差異。
4 討論
本研究提出了一種基于波形求和的磁聲信號求解方法。在準一維的聲傳播情況下，對具有不同厚度樣本的磁聲信號時頻特性進行了仿真計算和分析。仿真結果表明，磁聲信號的時頻特性包含樣品厚度的信息。在單周期正弦脈沖激勵下，厚層樣本信號波形半寬度大于薄層樣本。
實驗結果驗證了波形求和方法的仿真和分析結果。薄層樣本和厚層樣本的電導率邊界對應的波形和頻譜具有明顯不同的特征。由于常見窄帶超聲波換能器脈沖響應為“W”型波形，薄樣品的聲信號呈現“W”型，而厚樣品的聲信號為具有正負峰的正弦形脈沖。厚層樣本的邊界聲源可以通過薄樣品聲源的求和計算得到。
采用波形求和方法分析磁聲波信號，可以根據傳統的多物理場有限元法得到聲源空間分布，進而獲得時域解。該方法可以求解任意幾何形狀樣本任意激勵波形的磁聲耦合正問題，能更準確地反映聲信號的時頻特性。為后續根據磁聲信號時頻特性，反演重建聲源分布特征和電導率圖像重建提供重要的依據。
然而，本文研究工作中仍存在需要解決的問題。一方面如上所述，為了獲得反映兩個邊界和介質內部信息的磁聲信號，應開發先進的寬帶測量系統。同時，考慮到商用超聲換能器是窄帶換能器，因此，可通過設計寬帶換能器或采用多頻帶傳感器檢測技術，從而獲得包含寬頻帶信息的磁聲信號[23]。另一方面，聲信號幅值微弱，容易受到噪聲的干擾，會導致檢測信號的失真。因此，應考慮具有高信噪比的檢測技術，如鎖相放大方法和技術的利用等。實驗裝置還需要進一步完善，以抑制噪音的干擾。另外，本文僅對一維簡單模型進行了聲源信號研究，對于多層復雜的二維、三維介質電導率分布模型，還需要進一步分析和研究，從而為二維、三維圖像重建的應用提供依據。最后，在上述聲源厚度與時頻特性研究的基礎上，在圖像重建過程中，根據聲源時頻特性，對介質厚度反演重建，將對獲取介質內部聲源的電導率圖像具有重要意義。
總之，本文對磁聲信號的時頻特性研究，為區分不同厚度磁聲成像樣本邊界提供了一種新的思路和方法，為進行小于一個波長精細結構的成像分析提供依據。
引言
生物組織電特性反映了組織的生理病理特征[1]，研究表明，正常組織與腫瘤病變組織電導率具有明顯差異[2-4]，因此組織電特性的檢測和成像將為相關疾病的的早期診斷提供重要依據。磁聲耦合效應的生物組織電特性功能成像方法[5-9]，由于具有無創性，造價較低，能夠反映生物組織早期病變，同時兼具電阻抗成像高對比度和超聲成像高空間分辨率的優點，成為了近些年研究的熱點。
目前國內外有多篇文獻對磁聲成像方法進行研究報道[10-19]，由于該方法涉及電、磁、聲、信號與系統等多學科交叉，目前還處于實驗室研究階段，相關機制以及檢測成像技術研究以仿體樣本為主。其中磁聲信號正問題求解采用有限元方法(finite elements method,FEM)，設置激勵為沖激函數[11]，則正問題磁聲信號脈沖即對應空間電導率的邊界。在磁聲信號處理和逆問題電導率重建過程中[12-13]，信號經放大濾波和數字平均處理去除隨機噪聲[14]后，進行聲源分布重建，然后在已知介質為實心體的情況下，根據介質邊界聲源的幅值推導介質的相對電導率圖像。然而實際測量中，由于儀器設備帶寬限制，難以實現無限窄脈沖的激勵，使得實測磁聲信號脈沖特性復雜。因此成像分辨率也受到影響，即受限于激勵脈沖寬度，無法分辨小于脈沖寬度的成像目標，如比較薄的結締組織、隔膜、血管壁等，其前后邊界聲源發出的聲信號發生混疊，即兩個或多個邊界聲源脈沖信號疊加，使得逆問題研究中，在未知介質電導率分布的情況下，無法分辨介質為實心電導組織，還是薄層膜狀組織，從而影響圖像重建的準確性。
因此，本文利用波形求和的正問題求解方法，對不同厚度的電導率介質的磁聲信號時頻特性進行研究，從磁聲信號中提取并分析比較薄層介質的電導率聲源和厚層介質的邊界聲源信息的異同點，為在未知介質電導率分布的情況下由磁聲信號分辨和重建內部聲源電導率圖像提供依據。
1 原理
1.1 磁聲耦合原理
磁聲成像方法通過對介質施加電磁激勵，介質內部質點受到洛倫茲力作用發生振動，形成聲波，在介質外部可由聲換能器接收聲信號響應，如圖 1。磁聲成像滿足波動方程[15]：
| ${\nabla ^2}p\left({\boldsymbol{r},t} \right)- \frac{1}{{c_x^2}}\frac{{\partial p\left({\boldsymbol{r},t} \right)}}{{\partial {t^2}}}=\nabla \left({\boldsymbol{r},t} \right) \cdot \left[{\boldsymbol{J}\left({\boldsymbol{r},t} \right) \times {\boldsymbol{B}_0}\left(\boldsymbol{r} \right)} \right]$ |
其中p(r,t)為r處聲信號，J為電流密度,B0為靜磁場的磁感應強度，R為r到傳感器的距離,cs介質內部聲傳播速度。
由格林函數求解波動方程，介質外部聲信號響應[15]
| $\begin{array}{l} p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left[{\sigma \left({\boldsymbol{r}'} \right.\boldsymbol{E}\left({\boldsymbol{r}'} \right) \times {\boldsymbol{B}_0}\left({\boldsymbol{r}'} \right)} \right]} \\ \frac{{\delta \left({t- R/{c_s}} \right)}}{R} \end{array}$ |
其中▽r′·[σ(r′)E(r′)×B0(r′)]為聲源，σ為介質電導率，E為介質內部電場，洛倫茲力為f=σ(r′)E(r′)×B0(r′)。其中延遲項包含了聲振動源到傳感器位置r的傳播距離引起的傳播延遲[16-17]。
因此，r處檢測的磁聲信號p(r,t)可表示為介質內部聲源在不同傳輸延遲下的積分。
圖1
磁聲耦合裝置原理圖
Figure1.
Scheme of magneto-acoustic effect experimental set-up
1.2 波形求和算法
由前所述，電流密度J可分離為空間項和時間項的乘積[15]，為了計算磁聲信號的時域波形，我們將計算過程分為2步：空間項求解和時間項求解，即計算介質內電流密度的空間分布，再對各剖分單元分別考慮其隨時間變化的信號波形，根據各剖分單元電流密度模和方向，利用磁聲信號格林函數解，求解磁聲信號。
考慮簡化的矩形介質模型，設電流沿y方向，靜磁場沿－z方向，則洛倫茲力方向沿x方向，這時聲信號可視為沿一維方向(x軸)傳播，如圖 1所示。設檢測器端面與樣本前端面距離為d，樣本長度l，設樣本尺寸位于0≤x≤l區間內，則接收到的磁聲信號為整個介質各剖分單元聲信號求和的結果。
| $p=\sum\limits_{i=1}^n {{p_i}} $ |
其中pi為第i個剖分單元的聲壓。
由于聲源項▽r′·(J×B0)中電流密度J可以分離為空間量和時間量乘積的形式[15],則聲源表達式(2)可變形為
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left({\boldsymbol{J} \times {\boldsymbol{B}_0}} \right)\frac{{\delta \left({t- R/{c_s}} \right)}}{R}} \otimes s\left(t \right)$ |
其中s(t)為電流激勵信號,應用卷積性質和沖激函數的抽樣性質，式(4)可寫為
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\int\limits_V {{\rm{d}}\boldsymbol{r}'\nabla r' \cdot \left({\boldsymbol{J} \times {\boldsymbol{B}_0}} \right)\frac{{s\left({t- R/{c_s}} \right)}}{R}} $ |
將積分過程離散為求和形式
| $p\left({\boldsymbol{r},t} \right)=\frac{{- 1}}{{4\pi }}\sum\limits_i {\frac{{\rm{d}}}{{{\rm{d}}y}}} \left({{\boldsymbol{J}_i}{\boldsymbol{B}_0}} \right)\cos {\theta _i}\frac{{s\left({t- dela{y_i}} \right)}}{{dela{y_i} \cdot {c_s}}}$ |
其中θi為剖分單元i受到洛倫茲力方向與傳感器軸線方向的夾角，式(6)可表示為以下形式
| $p\left({\boldsymbol{r},t} \right)=\frac{{{\boldsymbol{B}_0}}}{{4\pi }}\sum\limits_{i=1}^n {\frac{{\rm{d}}}{{{\rm{d}}y}}} {\boldsymbol{J}_i}\frac{{s\left({t- dela{y_i}} \right)}}{{dela{y_i} \cdot {c_s}}}\cos {\theta _i}$ |
應用式(7),則可計算出由介質全部剖分單元形成的磁聲信號，通過設置不同的介質厚度l，則可獲得其對應的磁聲信號波形。
2 方法
2.1 仿真方法
仿真研究中設電流方向沿y軸，靜磁場沿z軸，洛倫茲力即聲振動方向沿x軸，如圖 1所示。本研究利用有限元工具Comsol Multiphysics，將介質聲源剖分為若干單元，計算空間電流密度，再將結果導出，并由MATLAB計算，獲得檢測位置r處的磁聲信號，仿真計算過程如圖 2。
考慮矩形介質的準一維的聲傳播情況，本研究中選擇樣本厚度分別為l<λ/4、l=λ/4、λ/4<l<λ/2、l=λ/2、λ/2<l<λ、l=λ以及l>λ幾種情況進行仿真：仿真中樣本厚度分別設為0.6、0.75、1.125、1.5、2.25、3和20 mm；電導率均為1 S/m，聲學參數均勻；通過電極對樣本施加電流激勵，激勵為一個周期正弦信號。
| $\boldsymbol{J}\left(t \right)=A\sin 2\pi ft{a_x}\left({0 \le 2\pi ft \le 2\pi } \right)$ |
其中A為信號幅值，頻率f=500 kHz。
聲檢測器與樣本前端面距離為d，檢測器中心頻率為500 kHz；圖 2中顯示了計算步驟。有限元剖分網格尺寸小于0.12 mm。
圖2
磁聲信號正問題計算過程
Figure2.
Computing procedure of magneto-acoustic signal
2.2 實驗方法
為驗證不同厚度磁聲信號時頻特性，本研究采用磁聲成像實驗系統[19-20]，和明膠仿體樣本進行實驗。裝置如圖 1，由函數發生器(AFG3252，美國泰克)產生激勵信號，通過雙極性放大器(HSA4104，日本NF)放大，將激勵信號加載到實驗樣本上。由電磁鐵產生恒穩磁場，磁場在半徑50 mm、高度100 mm的圓柱形空間范圍內均勻。由超聲傳感器(V303，美國Olympus，中心頻率500 kHz)接收聲信號。超聲傳感器接收的聲信號經由放大器(5307，NF)100倍放大后，由示波器(MSO4000，美國泰克)進行波形顯示。采用單周期正弦脈沖信號為激勵，頻率500 kHz，靜磁場0.4 T。檢測距離為60 mm，檢測波形平均512次，以去除隨機噪聲。樣本為長方體，采用明膠制作，混入5% NaCl，樣本截面為10mm×5mm，厚度分別為l=0.6 mm和l=22.5 mm。
3 結果
3.1 仿真結果
介質內部電流密度仿真結果如圖 3所示，其中的樣本厚度l分別為0.6、0.75、1.125、1.5、2.25、3和20 mm，寬度均為10 mm，各樣本電流沿y方向均勻分布。其對應檢測位置的波形和短時傅里葉變換頻譜如圖 4所示。
圖3
介質內電流密度(y方向)有限元仿真結果
樣本寬度均為10 mm，由圖可見，各樣本電流沿
the width of each sample is 10 mm. the result shows that the current in the samples with different thickness distributes along
為進一步比較分析不同厚度介質磁聲信號時頻特性，對各厚度樣本的脈沖半寬度、峰值頻率進行了比較，如表 1所示。
表1
不同樣本厚度的時頻特性
Table1.
Time-frequency characteristics of signals from samples with different thickness
時域上，由圖 4(a)(c)(e)(g)(i)(k)(m)可見，激勵信號為單周期正弦脈沖時，對于厚度l小于聲信號波長λ的薄層樣本，當厚度小于1/4聲信號波長時兩邊界對應脈沖信號發生混疊，形成兩個正峰一個負峰的波形，且隨著樣本厚度增加，脈沖逐漸展寬。對于厚度l大于聲信號波長λ的厚層樣本，當厚度大于波長2倍時，信號波形中兩脈沖分離，兩脈沖極性相反，說明電導率變化趨勢相反。由表 1第4列可見，激勵頻率500 kHz，其脈沖波長λ=v/f=1 500 ms-1/500 kHz=3 mm，當樣本厚度分別為0.6、0.75、1.125、1.5、2.25、3 mm時，對應脈沖的半寬度(大于波形峰值1/的寬度)，分別由激勵信號的0.66 μs變為0.67、0.68、0.73、0.83、0.99和0.99 μs。
頻域上，由短時傅里葉變換，即圖 4(b)(d)(f)(h)(j)(l)(n)和表 1最右列可見，各厚度的薄層樣本聲信號脈沖對應的峰值頻率由423 kHz變為406、397、369、336、273和223 kHz，說明薄層樣本峰值頻率隨樣本厚度增加逐漸降低；而厚層樣本(l=20 mm)的脈沖峰值頻率與激勵信號近似相同，均為423 kHz。
圖4
介質磁聲信號波形及頻譜仿真結果
(a)為樣本厚度
(a) is waveform of sample with thickness
仿真結果表明，樣本厚度小于激勵波長時，隨著厚度增加，脈沖展寬，峰值頻率對應下降；樣本厚度大于波長時，樣本聲信號與激勵波形和頻譜近似相同。此結果表明磁聲信號時頻特性中包含了介質電導率分布信息，尤其對于厚度大于波長的厚層樣本，與薄層樣本的時頻特性具有顯著區別。
3.2 實驗結果
實驗測量的薄層樣本(0.6 mm厚)和厚層樣本(22.5 mm厚)的磁聲信號及短時傅里葉變換頻譜如圖 5、圖 6所示。
圖5
實測樣本的磁聲信號及頻譜
(a)磁聲信號波形；(b)磁聲信號頻譜。樣本厚度
(a) measured waveform; (b) short-time Fourier transform spectrum of the waveform. The rectangular sample’s thickness
圖6
實測樣本的磁聲信號及頻譜
(a)磁聲信號波形；(b)(c)分別為圖中樣本前后邊界對應的脈沖w1、w2的短時傅里葉變換頻譜。樣本厚度
(a) measured waveform; (b) the short-time Fourier transform spectrum of the pulse w1 in (a); (c) the short-time Fourier transform spectrum of the pulse w2 in (a). The rectangular sample’s thickness
由圖 5、圖 6可見，對于相同的激勵信號，在相同的檢測距離下，厚層樣本和薄層樣本的磁聲信號具有顯著不同，對于薄層樣本，其波形呈現W型，其主要原因在于正弦波激勵信號經過聲換能器特性以及成像系統各儀器系統函數的卷積后形成[16-17]；而厚層介質呈現具有正負峰的正弦脈沖波形，如圖 5(a)，其脈沖波形展寬為原來的2倍左右(1.29 μs)，該現象與理論分析一致。
此外，比較圖 4(m)與圖 5(a)可見，仿真結果與實驗具有一定差異，這是由于實驗中介質內部的反射散射以及噪聲的影響，檢測通路的頻率響應限制，使得波形產生一定失真[21]。為此，對實驗數據進一步分析處理。首先對實驗信號進行小波去噪，采用1維dmey小波5層分解，去除雜波成分，并將薄層樣本的測量信號作為一個剖分單元，則一個厚層介質可視為若干個薄層樣本的組合，因此利用一定數量薄層介質的求和，即可獲得對應樣本厚度的厚層介質的磁聲信號。由于薄層樣本的實驗測量信號中包含了成像系統響應函數的卷積，根據卷積的線性性質，采用此處理方法求和計算的厚層樣本磁聲信號同樣包含了成像系統響應函數，因此去掉系統響應函數影響，便于對薄、厚樣本的時頻特性進行分析驗證。結果如圖 7所示。由實驗測量的薄層介質聲信號計算的厚層介質聲信號如圖 7(a)，對信號波形左右兩個脈沖w4、w5進行短時傅里葉變換如圖 7(b)(c)，計算波形短時的峰值頻率可知，薄層仿體信號峰值頻率為459 kHz。而由薄層疊加計算的左右脈沖峰值頻率分別為305 kHz和386 kHz，說明由于波形疊加導致信號的峰值頻率向低頻移動，與仿真結果一致。而對于實測厚層介質峰值頻率為310 kHz和393 kHz，其頻率誤差為1.6%和1.8%，說明疊加處理方法計算的結果與實測結果近似，驗證了波形求和求解及薄厚層時頻特性分析方法的可行性。
圖7
由薄層樣本磁聲信號計算得到的厚層樣本的磁聲信號
(a)為模擬計算獲得的厚層樣本磁聲信號波形；(b)(c)為w4、w5兩個信號的短時傅里葉變換
Figure7. Magneto-acoustic signal of thick sample calculated from thin sample(a) the calculated signal using the measured pulse in Fig. 5; (b) the short-time Fourier transform spectrum of the pulse w4 in (a); (c) the short-time Fourier transform spectrum of the pulse w5 in (a)
為進一步分析時域波形的相似程度，采用相關系數進行計算。相關系數通常用來表征波形的相似程度[22]。本研究選擇實驗波形的前后邊界對應的兩個脈沖46.61~51.76 μs(令其為w1)，60.96~65.75 μs(令其為w2)，同時選擇兩個脈沖中間部分的55.43~59.53 μs(令其為w3)的干擾波形，如圖 6(a)所示，與處理后波形中的脈沖，即圖 7(a)中的w4、w5，分別進行相關系數計算。結果顯示，實驗測量的仿體前后邊界對應波形w1、w2與由薄層波形疊加計算w4、w5的相關系數，分別為97.96%和76.04%，可見左側波形w1相似度較高，而右側仿體邊界w2相似度降低，這是由于聲信號在介質中傳播同時與散射等干擾的疊加結果，導致波形產生變化。而對于干擾脈沖，相關系數計算結果表明，介于仿體兩個邊界中間的內部干擾信號w3，其峰值頻率131 kHz，與計算得到的仿體左右側峰w1、w2波形的相關系數為29.88%和23.85%。可見，對于非介質邊界產生的聲信號，其時頻特性與介質邊界具有差異。
4 討論
本研究提出了一種基于波形求和的磁聲信號求解方法。在準一維的聲傳播情況下，對具有不同厚度樣本的磁聲信號時頻特性進行了仿真計算和分析。仿真結果表明，磁聲信號的時頻特性包含樣品厚度的信息。在單周期正弦脈沖激勵下，厚層樣本信號波形半寬度大于薄層樣本。
實驗結果驗證了波形求和方法的仿真和分析結果。薄層樣本和厚層樣本的電導率邊界對應的波形和頻譜具有明顯不同的特征。由于常見窄帶超聲波換能器脈沖響應為“W”型波形，薄樣品的聲信號呈現“W”型，而厚樣品的聲信號為具有正負峰的正弦形脈沖。厚層樣本的邊界聲源可以通過薄樣品聲源的求和計算得到。
采用波形求和方法分析磁聲波信號，可以根據傳統的多物理場有限元法得到聲源空間分布，進而獲得時域解。該方法可以求解任意幾何形狀樣本任意激勵波形的磁聲耦合正問題，能更準確地反映聲信號的時頻特性。為后續根據磁聲信號時頻特性，反演重建聲源分布特征和電導率圖像重建提供重要的依據。
然而，本文研究工作中仍存在需要解決的問題。一方面如上所述，為了獲得反映兩個邊界和介質內部信息的磁聲信號，應開發先進的寬帶測量系統。同時，考慮到商用超聲換能器是窄帶換能器，因此，可通過設計寬帶換能器或采用多頻帶傳感器檢測技術，從而獲得包含寬頻帶信息的磁聲信號[23]。另一方面，聲信號幅值微弱，容易受到噪聲的干擾，會導致檢測信號的失真。因此，應考慮具有高信噪比的檢測技術，如鎖相放大方法和技術的利用等。實驗裝置還需要進一步完善，以抑制噪音的干擾。另外，本文僅對一維簡單模型進行了聲源信號研究，對于多層復雜的二維、三維介質電導率分布模型，還需要進一步分析和研究，從而為二維、三維圖像重建的應用提供依據。最后，在上述聲源厚度與時頻特性研究的基礎上，在圖像重建過程中，根據聲源時頻特性，對介質厚度反演重建，將對獲取介質內部聲源的電導率圖像具有重要意義。
總之，本文對磁聲信號的時頻特性研究，為區分不同厚度磁聲成像樣本邊界提供了一種新的思路和方法，為進行小于一個波長精細結構的成像分析提供依據。
