[發明專利]一種基于非均勻背景介質的彈性波成像方法

有效

• 申請號: 201810906598.9
• 申請日: 2018-08-10
• 公開（公告）號: CN109239771B
• 公開（公告）日: 2020-01-31
• 發明（設計）人: 徐魁文;楚彥青;趙文生;陳世昌;趙鵬;王高峰
• 申請（專利權）人: 杭州電子科技大學
• 主分類號: G01V1/28
• 分類號: G01V1/28
• 代理公司: 33240 杭州君度專利代理事務所(特殊普通合伙)
• 代理人: 黃前澤

1.一種基于非均勻背景介質彈性波成像方法，彈性波為電磁波，其特征在于包括以下：

步驟(1)、根據離散的網格位置和發射裝置、接收裝置的位置，計算和并根據非均勻背景的對比度和場積分方程(4)-(6)計算出相應的非均勻背景總場場強散射場場強和對比源

其中為離散后的格林函數G(r_s，r′)的積分算子；是離散后的格林函數G(r，r′)的積分算子；表示一個位于空間r_s處的點源對其周圍空間某一點r′所產生的場；表示一個位于空間r處的點源對其周圍空間某一點r′所產生的場；為第一類零階漢克爾函數，i表示虛數，k₀是彈性波的波數；具體是：

總場積分方程：

其中表示位于r處的入射場場強；χ(r′)＝(∈(r′)-∈₀)/∈₀，為∈_r的對比度函數，∈₀表示彈性波穿過的介質的某種物理特性；L為發射裝置的個數；

散射場積分方程：

表示的是位于r_s處的接收裝置所接收到的散射體產生的散射場數據，M為接收裝置的個數；

對比源為對比度和總場的乘積，定義為：

I_l(r)＝χ(r)E_l(r) (3)

將公式(1)-(3)離散化得到：

第(m，n)個離散網格的總場場強：

散射場場強：

對比源：

其中(m，n)代表離散網格的中心坐標，為第(m，n)個離散網格的感應電流，為各個離散網格感應電流的集合，是離散后的格林函數G(r，r′)的積分算子，為離散后的格林函數G(r_s，r′)的積分算子，是χ(r′)的對角矩陣形式；

由于非均勻背景介質的對比度已知，因此進一步可以得到：

非均勻背景總場場強

非均勻背景散射場場強

非均勻背景對比源

步驟(2)、初始化未知散射體對比度和正則化參數α，并對它們賦值初值為0，同時設定迭代次數p＝0；

步驟(3)、將未知散射體對比度代入到公式(10)和(11)中，得到未知散射體的對比源和未知散射場理論值F(Δχ)，并求取此次迭代獲得的未知散射場理論值F(Δχ)和未知散射體對比度之間的雅可比矩陣D，并進行奇異值分解；具體是：

由于非均勻背景介質的對比度已知，因此可將探測區域的對比度和總場以及對比源分為已知的非均勻背景介質和背景介質里面的未知散射體這兩部分：

其中分別表示非均勻背景的對比度、總場、對比源，分別表示未知散射體的對比度、總場、對比源；

將公式(7)、(8)、(9)代入公式(4)、(6)則得到關于未知散射體對比源的對比源積分方程：

將公式(10)代入到公式(5)中，則可以得到關于未知散射體的散射場場強理論值，也即目標函數：

其中⊙的定義為兩個矩陣的乘法，vec{}定義為向量化張量的操作，公式(11)表明目標函數僅僅只是關于未知散射體的對比度的函數；

根據e^s表示已知非均勻背景的散射場場強數據和未知散射體的散射場測量數據之和，即接收裝置直接接收到的數據，表示已知非均勻背景的散射場場強數據，由于也是已知，故而得到背景里面的未知散射體的散射場實測數據：

根據上述構建未知散射體的對比度的成本函數：

f(Δχ)＝||F(Δχ)-Δe^s||²_min (12)

由于公式(12)的病態性，利用改進增強的Levenberg-Marquardt方法進行優化，可以得到以下方程：

[D*D+αI]Δ(Δχ_α)＝D*δe^s (13)

其中α表示正則化參數，I是單位矩陣；δe^s＝Δe^s-F(Δχ_c)表示測量的未知散射體的散射場測量數據與利用公式(11)計算的未知散射體的散射場場強理論值之間的差值；其中Δχ_c表示本次迭代的數值結果；

進而可以得到：

Δ(Δχ_α)＝[D*D+αI]^-1D*δe^s (14)

由于公式(14)計算量比較大，并且不容易濾除數據噪聲，故而采用奇異值分解的方法來求Δ(Δχ_α)；所述的奇異值分解具體是：

采用雅可比矩陣D＝U∑V^*求解方程(14)，得到前后相鄰兩次迭代未知散射體的對比度誤差：

進一步獲取下一次迭代的對比度參數：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (16)

其中V是大小為N×N的單位矩陣，表示雅可比矩陣D的右奇異矩陣；∑為LM×N的對角矩陣，[∑]_k＝σ_k，k＝1，2，...，min(LM，N)，σ_k表示的是∑第k行k列上的對角矩陣元素；

步驟(4)、利用公式(11)重新計算下一次迭次的未知散射體的散射場數據F(Δχ_p+1)，并求出未知散射體的散射場場強的理論值和測量值之間的差值δe^s，然后判斷是否滿足迭代停止條件，如果滿足條件則結束并輸出未知散射體最佳的對比度值進而重建相應的圖像，如果不滿足則進行步驟(5)，繼續進行迭代優化；

步驟(5)、更新迭代次數p＝p+1，利用公式(17)函數V(α)最小值時所對應的α值計算新的正則化參數α_p，返回步驟(3)繼續進行優化；

改進的廣義交叉函數如下所示：

其中，U是大小為LM×LM的單位矩陣，表示雅可比矩陣D的左奇異矩陣；n是矩陣U的行數，q是矩陣U的列數；K為雅可比矩陣D奇異值截取的個數。

2.一種基于非均勻背景介質彈性波成像方法，彈性波為聲波，其特征在于包括以下：

步驟(1)、根據離散的網格位置和發射裝置、接收裝置的位置，計算和并根據非均勻背景的對比度和場積分方程(20)-(21)計算出相應的非均勻背景總場場強P^b，散射場場強P^b，sca；

其中為離散后的格林函數g(r_s，r′)的積分算子，是離散后的格林函數g(r，r′)的積分算子，表示一個位于空間r_s處的點源對其周圍空間某一點r′所產生的場；表示一個位于空間r處的點源對其周圍空間某一點r′所產生的場；為第一類零階漢克爾函數，i表示虛數，k₀是彈性波的波數；具體是：

總場積分方程：

散射場積分方程：

將以上兩個公式離散化之后可得：

其中χ₂均為密度的對比度函數值，以的形式表示χ的離散化形式，其它公式字母也均以此方式表示，P^inc、P、P^sca分別表示入射場、總場和散射場，由于非均勻背景介質的對比度已知，因此進一步可以得到：

非均勻背景總場場強P^b：

非均勻背景散射場場強：P^b，sca：

步驟(2)、初始化未知散射體對比度Δχ₁，Δχ₂和正則化參數α，并對它們賦值初值為0，同時設定迭代次數p＝0；

步驟(3)、將Δχ₁，Δχ₂代入到公式(27)和(28)中，得到F(Δχ₁，Δχ₂)，并求取此次迭代獲得的未知散射場數據和未知散射體對比度之間的雅可比矩陣D，并進行奇異值分解；具體是：

由于非均勻背景介質的對比度已知，因此可將探測區域的對比度和總場P分為已知的非均勻背景介質和背景里面的未知散射體這兩部分：

P＝P^b+ΔP (26)

其中表示非均勻背景的對比度，P^b表示非均勻背景的總場；表示非均勻背景里面的未知散射體的對比度，ΔP表示非均勻背景里面的未知散射體的總場；將公式(24)，(25)，(26)代入公式(20)，(22)中，并相減，得到關于的場差分方程：

消去中間變量ΔP，并代入公式(21)得到關于未知散射體散射場的理論值，也即目標函數：

公式(28)表明目標函數僅僅只是關于未知散射體的對比度的函數，

根據Δe^s＝e^s-P^b，sca，e^s表示已知的非均勻背景的散射場數據和未知散射體的散射場的測量數據之和，也即接收裝置直接接收到的數據，P^b，sca表示已知的非均勻背景散射體的散射場測量數據，由于P^b，sca也是已知，故而得到背景里面的未知散射體的實測數據：

根據上述構建未知散射體的對比度的成本函數，令Δχ為Δχ₁，Δχ₂的集合，則：Δχ＝[Δχ₁，Δχ₂]，則F(Δχ₁，Δχ₂)＝F(Δχ)，構建成本函數為：

f(Δχ)＝||F(Δχ)-Δe^s||²_min (29)

由于公式(29)的病態性，因此利用改進增強的Levenberg-Marquardt方法來優化，得到以下方程：

[D*D+αI]Δ(Δχ_α)＝D*δe^s (30)

進而可以得到：

Δ(Δχ_α)＝[D*D+αI]^-1D*δe^s (31)

由于公式(31)計算量比較大，并且不容易濾除數據噪聲，故而采用奇異值分解的方法來求Δ(Δχ_α)，所述的奇異值分解具體是：

采用雅可比矩陣D＝U∑V^*求解方程(30)，得到前后相鄰兩次迭代未知散射體的對比度誤差：

進一步獲取下一次迭代的對比度參數：

Δχ_p+1＝Δχ_p+Δ(Δχ_α)_p (33)

其中U是大小為LM×LM的單位矩陣，它表示雅可比矩陣D的左奇異矩陣；V是大小為N×N的單位矩陣，它表示雅可比矩陣D的右奇異矩陣；∑為LM×N的對角矩陣，[∑]_k＝σ_k，k＝1，2，...，min(LM，N)，σ_k表示的是∑第k行k列上的對角矩陣元素，

步驟(4)、利用公式(28)重新計算下一次迭次的未知散射體的散射場數據F(Δχ_p+1)，并求出散射場的理論值和測量值之間的差值δe^s，判斷是否滿足迭代停止條件，如果滿足條件則結束并輸出未知散射體最佳的對比度值進而重建相應的圖像，如果不滿足，則進行步驟(5)，繼續進行迭代優化；

步驟(5)、更新迭代次數p＝p+1，利用公式(34)計算新的正則化參數α_p，返回步驟(3)繼續進行優化；

改進的廣義交叉函數如下所示：

其中，n是矩陣U的行數，q是矩陣U的列數；K為雅可比矩陣D奇異值截取的個數；當函數V(α)取得最小值時的α值，即為所要選取的正則化參數。