一種優(yōu)化全波反演法重建聲速圖像的方法
【專利摘要】本發(fā)明公開了一種優(yōu)化全波反演法重建聲速圖像的方法，包括以下步驟：利用超聲數(shù)據(jù)采集設備，得到傳感器所在位置的實測聲壓信號；在計算機上完全模擬發(fā)射、接收過程作為正向函數(shù)，該函數(shù)輸入變量為聲速圖矩陣，輸出變量為傳感器對應位置所記錄的聲壓信號與的實測聲壓信號的差的范數(shù)；進入迭代過程：給出一個聲速圖像作為初始圖像，使其作為自變量帶入到上述正向函數(shù)中，得到輸出范數(shù)；計算出正向函數(shù)在當前聲速圖位置的一階梯度圖；通過擬牛頓法解決最優(yōu)化問題，更新原有聲速圖，得到新的聲速圖；將新的聲速圖作為初始圖像重復迭代過程，直到得到滿足條件退出迭代過程，并得到當前的聲速圖像。
【專利說明】
一種優(yōu)化全波反演法重建聲速圖像的方法
技術領域
[0001] 本發(fā)明屬于醫(yī)學超聲成像及處理領域，特別是一種針對聲速圖像重建的優(yōu)化方 法。
【背景技術】
[0002] 經(jīng)過幾十年的研究，全波反演重建技術在醫(yī)學領域的應用也日趨深入。例如，在廣 泛應用于乳腺癌檢測的超聲檢測層析成像的重建過程中，全波反演的重建方法得到的圖像 具有極具優(yōu)化的空間分辨率，這一點也強過傳統(tǒng)的Ray-Based重建技術，Ray-Based將聲波 傳播假設為具有光的折射特性，所以不能重建出全波形所具有的特征。然而由于全波反演 方法本身需要相當長的計算時間和巨大的計算機存儲資源，這項技術沒有被廣泛的應用于 醫(yī)學臨床應用。因此，需要一種優(yōu)化的方法來提高全波反演重建聲速圖像方法的效率。

【發(fā)明內容】

[0003] 發(fā)明目的：本發(fā)明所要解決的技術問題是針對現(xiàn)有全波反演重建聲速圖像速度 慢，效率低的不足，提供了一種從算法的角度提高重建速度和效率，并保證圖像質量的優(yōu)化 方法。
[0004] 為了解決上述技術問題，本發(fā)明公開了一種利用擬牛頓法降低運算量并提高效率 的全波反演重建聲速圖像方法，包括如下步驟：
[0005] 步驟一，利用超聲數(shù)據(jù)采集設備，依次完成Μ個信號源單獨發(fā)射、所有的Μ個傳感器 接收的過程，直到所有信號源都發(fā)射過一次，得到傳感器所在位置的實測聲壓信號g m(r， t)，其中r為傳感器坐標向量參數(shù)，t為時間參數(shù)，m表示發(fā)射次序，m=l，2,. . .，M，M同時也表 示總共的發(fā)射次數(shù)；
[0006] 步驟二，根據(jù)步驟一所用設備在計算機上完全模擬步驟一中的發(fā)射、接收過程作 為正向函數(shù)，使之輸入變量為聲速圖矩陣c，輸出變量為傳感器對應位置所記錄的聲壓信號 gm(r，t)與步驟一中的實測聲壓信號gm(r，t)的差的范數(shù)L;
[0007] 步驟三，開始迭代過程，首次迭代中給出聲速圖像CQ作為初始圖像，也即c = CQ，作 為步驟二中的正向函數(shù)的輸入量，得到輸出變量L;
[0008] 步驟四，計算出正向函數(shù)在聲速圖c位置的一階梯度圖J;
[0009] 步驟五，通過擬牛頓法解決最優(yōu)化問題，利用步驟四中的梯度圖J更新聲速圖c，得 到新的聲速圖；
[0010] 步驟六，將新的聲速圖作為初始圖像重復迭代過程，重復步驟三到六，直到得到滿 足條件退出迭代過程；
[0011] 步驟七，得到滿足步驟六的條件時對應的聲速圖像c即為優(yōu)化的全波反演法重建 的聲速圖像。
[0012] 本發(fā)明中，優(yōu)選地，步驟一中所用信號源用于超聲波的發(fā)射，傳感器用做聲壓信號 也即gm的記錄、保存，其中下標m表示發(fā)射次序所有的信號源和傳感器成對地位于同一位 置。
[0013] 本發(fā)明中，優(yōu)選地，步驟一中發(fā)射接收過程遵循以下原則:每次只有一個信號源發(fā) 射，其余傳感器(包括其本身)接收，傳感器共Μ個，那么發(fā)射次數(shù)也應為Μ次。
[0014] 本發(fā)明中，優(yōu)選地，步驟一中選用環(huán)形傳感器圍繞目標，或雙排線型傳感器置于目 標兩端，共有Μ個傳感器依次排列，其位置坐標可由r對應取值唯一表示。
[0015] 本發(fā)明中，優(yōu)選地，步驟二中的正向函數(shù)模擬了計算區(qū)域、傳播介質、發(fā)射設備、接 收設備以及聲波傳播過程。
[0016] 具體而言，正向函數(shù)具體包括定義了計算區(qū)域、傳播介質、發(fā)射設備、接收設備。計 算區(qū)域的定義是傳播介質的定義的基礎，二者可以表示為大小一致的矩陣。定義計算區(qū)域 大小 n = NXNXNz，表示將傳播介質劃分為n個區(qū)域，我們稱其每個區(qū)域為一個格點，其中N 和Νζ分別表示長或寬和厚度。傳播介質在其基礎上具體化了每個格點的參數(shù)，包括但不限 于聲速、密度、衰減系數(shù)等信息。發(fā)射設備、接收設備的兩種功能由信號源和傳感器來完成， 信號源和傳感器位置相同并位于計算區(qū)域中，計算區(qū)域的邊界經(jīng)過處理后，聲波會在邊界 被完全吸收。
[0017] 本發(fā)明中，優(yōu)選地，步驟二中的聲速圖為一個矩陣，也即一個數(shù)字化的傳播介質模 型，聲速只是傳播介質的一部分信息，其余信息待定為一般參數(shù)，一般是目標周圍介質的對 應參數(shù)，傳播介質中應包含目標信息，該聲速圖矩陣每一點的值為對應位置的聲速值。
[0018] 本發(fā)明中，優(yōu)選地，步驟二中的正向函數(shù)為模擬過程，正向函數(shù)輸入變量為聲速圖 矩陣，該函數(shù)可以通過給定的聲速圖直接得到模擬聲壓信息，進而可得范數(shù)L。
[0019] 本發(fā)明中，優(yōu)選地，步驟二中的正向函數(shù)，也即模擬的聲波傳播過程中，接收的聲 壓信號與波源、傳播介質滿足方程：
[0020]

[0021] 其中r為坐標向量參數(shù)，t為時間變量參數(shù)，pm(r，t)為聲壓，m表示發(fā)射的順序次 數(shù)，m=l，2，. ..，M，M表示總共的發(fā)射次數(shù)，sm(r，t)為超聲信號，c(r)為聲速圖矩陣，c(r)中 共有n = NXNXNz個格點。
[0022] 由于奈奎斯特定律限制，上述的超聲信號sm(r，t)的頻率不能超過 ，其 中co在表不聲速(m/s)，dx代表每一格點換算成實際的邊長(m)。
[0023] 本發(fā)明中，優(yōu)選地，步驟二中的正向函數(shù)，時間變量參數(shù)t在計算中需要離散化表 示。在傳感器點記錄聲壓信息的采樣時間間隔以dt(s)表示，T(s)表示時間總長度。
[0024]本發(fā)明中，優(yōu)選地，步驟二中的正向函數(shù)，也即模擬過程中，r可以為任意坐標，不 限于傳感器的位置。故而，在模擬過程中，不僅僅可以得到傳感器位置的聲壓信息，包括可 以得到整個計算區(qū)域所有位置的聲壓信息。
[0025] 本發(fā)明中，優(yōu)選地，步驟二中范數(shù)可以表示為：
[0026]
[0027]其中gm，gm*表表示步驟一得到的實測聲壓和在接下來的步驟三得到的模擬聲壓。 與以&，〇區(qū)^是，gm表示在測量位置的聲壓信號，也就是r等于傳感器位置坐標時的pm(r， t) ο
[0028] 本發(fā)明中，優(yōu)選地，步驟三中得到L，求解聲速圖的問題就轉化為在接下來的步驟 中解決如下表達式表述的問題：
[0029]
[0030] 該式表達了一個非線性的最優(yōu)化問題；即為最后要求得的矩陣。
[0031] 本發(fā)明中，優(yōu)選地，步驟三中除了得到輸出值L，函數(shù)產(chǎn)生的聲壓信息pm(r，t)也在 此步驟中得到并保留，并被用于接下來步驟中的運算。
[0032] 本發(fā)明中，優(yōu)選地，步驟四中c指的是當前聲速圖，如果為首次迭代，那么c = CQ，否 則應按上次迭代結果帶入。
[0033] 本發(fā)明中，優(yōu)選地，步驟四中的梯度的計算可以用伴隨狀態(tài)法(adjoint state method)計算并說明如下。首先，L計算式表述如下：
[0034]
[0035] 其中Ωη(ΝΧΝΧΝζ)表示整個包括目標在內的計算區(qū)域，T表示時間總長度。那么 根據(jù)伴隨狀態(tài)法和L的表述式，L的梯度計算式表述如下：
[0036]
[0037]其中qm(r，t)是下面波動方程的解。
[0038]
[0039] 其中1111(1'，1:)可以表述為1：111(1'，1:) = 8111(1'，1'-1:)-￡111(1'，1'-1:)。那么計算中1^關于(3的梯 度用J表達如下：
[0040]
[0041 ]其中η表示對應的聲速圖像中每一個離散化格點，同時qm，Pm對應qrn(r，t)和Pm(r， t)中坐標參數(shù)r為對應格點位置時的值。
[0042]本發(fā)明中，優(yōu)選地，步驟五的擬牛頓法兩種算法DFP法、BFGS法均可以在此應用，在 本問題中具體表述如下：
[0043] 擬牛頓法用近似的Hesse矩陣Bk代替真實的Hesse矩陣。那么Bk在每次迭代中的更 新公式如下，在DFP法中記=
[0044]
[0045] 其中，c(k)表示第k次迭代中的聲速圖，J(k)表示第k次迭代中的梯度圖。
[0046] 在BFGS法中，Bk的計算公式有所不同并區(qū)別如下：
[0047]
[0048] 那么迭代公式，也即最終得到的c(k+1)應為：
[0049] c(k+1) = c(k)+akPk
[0050] 其4
，ak為步長根據(jù)線性搜索法確定。具體而言線性搜索法 包括了進退法和黃金分割法。
[0051] 本發(fā)明中，優(yōu)選地，步驟六中利用迭代的過程不斷更新圖像的方法重建聲速圖像。
[0052] 本發(fā)明中，優(yōu)選地，步驟六中閾值的設定是相對的，在具體實驗分析中而定的，除 此以外，還可以設定固定的迭代次數(shù)。
[0053] 本發(fā)明中，針對全波反演法現(xiàn)有的低效問題，進行了一定的優(yōu)化，由于擬牛頓法充 分利用一階梯度而得到近似的二階梯度矩陣大大提高了效率，又在GHJ的應用下，計算時間 可以從一周左右縮短至幾個小時。
【附圖說明】
[0054]下面結合附圖和【具體實施方式】對本發(fā)明做更進一步的具體說明，本發(fā)明的有關敘 述和其他方面的優(yōu)點將會變得更加清楚。
[0055]圖1為本發(fā)明流程圖。
[0056]圖2為本發(fā)明系統(tǒng)截面示意圖。
【具體實施方式】
[0057] 如圖1所示，本發(fā)明公開了一種優(yōu)化全波反演法重建聲速圖像的方法，包括以下步 驟：
[0058] 步驟一，利用超聲數(shù)據(jù)采集設備，依次完成Μ個信號源單獨發(fā)射、所有的Μ個傳感器 接收的過程，直到所有信號源都發(fā)射過一次，得到傳感器所在位置的實測聲壓信號g m(r， t)，其中r為傳感器坐標向量參數(shù)，t為時間參數(shù)，m表示發(fā)射次序，m=l，2,. . .，M，M同時也表 示總共的發(fā)射次數(shù)；
[0059]步驟二，根據(jù)步驟一所用設備在計算機上完全模擬步驟一中的發(fā)射、接收過程作 為正向函數(shù)，使之輸入變量為聲速圖矩陣c，輸出變量為傳感器對應位置所記錄的聲壓信號 gm(r，t)與步驟一中的實測聲壓信號gm(r，t)的差的范數(shù)L;
[0060] 步驟三，開始迭代過程，首次迭代中給出聲速圖像CQ作為初始圖像，也即C = CQ，作 為步驟二中的正向函數(shù)的輸入量，得到輸出變量L;
[0061] 步驟四，計算出正向函數(shù)在聲速圖c位置的一階梯度圖J;
[0062] 步驟五，通過擬牛頓法解決最優(yōu)化問題，利用步驟四中的梯度圖J更新聲速圖c，得 到新的聲速圖；
[0063] 步驟六，將新的聲速圖作為初始圖像重復迭代過程，重復步驟三到六，直到得到滿 足條件退出迭代過程；
[0064] 步驟七，得到滿足步驟六的條件時對應的聲速圖像c即為優(yōu)化的全波反演法重建 的聲速圖像。
[0065] 如圖2所示，實驗中傳感器共Μ個在計算區(qū)域內環(huán)形排列。目標中心盡量置于計算 區(qū)域的中心位置。圖2所示意的結構既是實驗設備結構，也是重建過程中的正向函數(shù)的計算 區(qū)域及其相關信息的示意圖。
[0066] 本實例中，步驟一中選用了既可用于超聲波的發(fā)射，也可用做聲壓信息(也即指得 到聲壓信號gm)的記錄、保存?zhèn)鞲衅?。發(fā)射接收過程遵循以下原則:每次只有一個傳感器發(fā) 射，其余傳感器(包括其本身)接收，接收傳感器為Μ個，發(fā)射次數(shù)也為Μ次。
[0067]本實例中，步驟一中選用環(huán)形傳感器圍繞目標，其位置坐標可由r唯一表示。
[0068]本實例中，步驟二中的初始聲速圖選為一個矩陣，由普通的Ray-based方法得到。 該聲速圖矩陣每一點的值為對應位置的聲速值。
[0069]在本實例的首次迭代中，正向函數(shù)輸入變量為初始聲速圖矩陣co。
[0070] 本實例中，由于奈奎斯特定律限制，超聲信號s m ( r，t)的頻率不能超過
，其中co在表不聲速(m/s)，dx代表每一格點換算成實際的邊長(m)。故在本實例 中，發(fā)射超聲信號設定為頻率f = 〇. 8Mhz的正弦函數(shù)與高斯核函數(shù)的乘積，co = 151 Om/s并 滿足上述要求。格點長度大小設為波長的三分之一。具體實際的計算區(qū)域設定為N=256，Nz =5，N對應0.128m實際長度。那么dx = 0.5mm
[0071] 本實例中，范數(shù)表示為：
[0072]
[0073] 其中gm，gm*表表示步驟一得到的實測聲壓和在接下來的步驟三得到的模擬聲壓。
[0074] 本實例中，步驟四中的梯度的計算可以用伴隨狀態(tài)法(adjoint state method)計 算，首次迭代中c = co，否則根據(jù)迭代確定c的取值，L計算式表述如下：
[0075]
[0076] 其中ΩΜ(ΝΧΝΧΝζ)表示整個包括目標在內的計算區(qū)域，T表示時間總長度。T設定 為足夠超聲波傳播的實踐長度。dt的選取設定，
?據(jù)伴隨狀態(tài)法和L的表述 式，L的梯度計算式表述如下：
[0077]
[0078]其中qm(r，t)是下面波動方程的解。
[0079]
[0080] 其中1：111(1'，1:)可以表述為1：111(1'，1:) = 8111(1'，1'-1:)-￡111(1'，1'-1:)。那么計算中1^關于(3的近 似梯度用J表達如下：
[0081] rn=i lwJ? m=i ?=ζ
ι-?
[0082]其中η表示對應的聲速圖像中每一個離散化格點，同時qm，Pm對應qrn(r，t)和Pm(r， t)中坐標參數(shù)r為對應格點位置時的值。
[0083] 本實例中，步驟五的BFGS擬牛頓法在本問題中具體表述如下：
[0084] 在BFGS法中，Bk的計算公式如下：
[0085]
[0086] 那么迭代公式，也即最終得到的c(k+1)應為：
[0087] c(k+1) = c(k)+akPk
[0088] 其中，ak為步長根據(jù)線性搜索法確定。具體而言線性搜索法包括了進 退法和黃金分割法兩種方法共同應用，進退法求出區(qū)間范圍[a，b]，黃金分割法在區(qū)間[a， b]之間得出合適的步長ak。
[0089] 本實例中，步驟六中利用迭代的過程不斷更新圖像的方法重建聲速圖像。
[0090] 本實例中，步驟六中閾值的設定為1 X 10-9，每次迭代后保存結果。迭代過程中可以 隨時停止并查看結果，并且可以繼續(xù)上次的結果運行。
[0091] 本實例流程圖參照圖1。
[0092] 本發(fā)明提出了一種優(yōu)化全波反演法重建聲速圖像的方法，應當指出，步驟一種涉 及的實驗設備型號形式不對本專利構成限制；傳感器位置可以是環(huán)形也可以是其他形式， 具體使用何種分布形式不對本專利構成限制;模擬過程中的計算區(qū)域設定的大小等非關鍵 參數(shù)，不對本專利構成限制。應當指出，對于本技術領域的普通人員來說，在不脫離發(fā)明原 理的前提下還可以做出若干改進和潤飾，這些也應視為本發(fā)明的保護范圍。另外，本實例中 未明確的各組成部分均可用現(xiàn)有技術加以實現(xiàn)。
【主權項】
1. 一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，包括如下步驟： 步驟一，利用超聲數(shù)據(jù)采集設備，依次完成M個信號源單獨發(fā)射、所有的M個傳感器接收 的過程，直到所有信號源都發(fā)射過一次，得到傳感器所在位置的實測聲壓信號Sm(r，t)，其 中r為傳感器坐標向量參數(shù)，t為時間參數(shù)，m表示發(fā)射次序，m=l，2,. . .，M，M同時也表示總 共的發(fā)射次數(shù)； 步驟二，根據(jù)步驟一所用設備在計算機上完全模擬步驟一中的發(fā)射、接收過程作為正 向函數(shù)，使之輸入變量為聲速圖矩陣c，輸出變量為傳感器對應位置所記錄的聲壓信號gm (r，t)與步驟一中的實測聲壓信號￡m(r，t)的差的范數(shù)L; 步驟三，開始迭代過程，首次迭代中給出聲速圖像co作為初始圖像，也即C = CQ，作為步 驟二中的正向函數(shù)的輸入量，得到輸出變量L; 步驟四，計算出正向函數(shù)在聲速圖c位置的一階梯度圖J; 步驟五，通過擬牛頓法解決最優(yōu)化問題，利用步驟四中的梯度圖J更新聲速圖c，得到新 的聲速圖； 步驟六，將新的聲速圖作為初始圖像重復迭代過程，重復步驟三到六，直到得到滿足條 件退出迭代過程； 步驟七，得到滿足步驟六的條件時對應的聲速圖像c即為優(yōu)化的全波反演法重建的聲 速圖像。2. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 一中發(fā)射接收過程遵循以下原則:M個傳感器和M個信號源成對地位于同一位置，每次只有 一個信號源發(fā)射，所有傳感器接收，傳感器共M個，那么發(fā)射次數(shù)也應為M次。3. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 二中的正向函數(shù)模擬了計算區(qū)域、傳播介質、發(fā)射設備、接收設備以及聲波傳播過程。4. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 二中的正向函數(shù)，也即模擬過程中，接收到的聲壓信號與波源、傳播介質滿足方程：其中r為坐標向量參數(shù)，t為時間變量參數(shù)，Pm( r，t)為聲壓，m表示發(fā)射的順序次數(shù)，m = 1，2, . . .，M，M表示總共的發(fā)射次數(shù)，sm(r，t)為超聲信號，c(r)為聲速圖矩陣，c(r)中共有n =N X NX Nz個格點，時間變量參數(shù)t在計算中需要離散化表示;在傳感器點記錄聲壓信息的 采樣時間間隔以dt表示，T表示時間總長度。5. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 二中范數(shù)可以表示為：其中表表示步驟一得到的實測聲壓和在接下來的步驟三得到的模擬聲壓;與^ (r，t)區(qū)別是，gm表示在測量位置的聲壓信號，也就是r為傳感器位置時的pm(r，t)。6. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 三中得到L，求解聲速圖的問題就轉化為在接下來的步驟中解決如下表達式表述的問題：該式表達了一個非線性的最優(yōu)化問題；即為最后要求得的矩陣。7. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 三中除了得到輸出值L，函數(shù)產(chǎn)生的聲壓信息p m(r，t)也在此步驟中得到并保留，并被用于 接下來步驟中的運算。8. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 四中c指的是當前聲速圖，如果為首次迭代，那么c = c〇。9. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步驟 四中的梯度的計算可以用伴隨狀態(tài)法計算。10. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步 驟五的擬牛頓法兩種算法DFP法、BFGS法均可以在此應用，步長根據(jù)線性搜索法確定，具體 而言線性搜索法包括了進退法和黃金分割法。11. 根據(jù)權利要求1所述的一種優(yōu)化全波反演法重建聲速圖像的方法，其特征在于，步 驟六中閾值的設定是相對的，在具體實驗分析中而定的，除此以外，還可以設定固定的迭代 次數(shù)，利用迭代的過程不斷更新圖像的方法重建聲速圖像。
【文檔編號】A61B8/00GK106037795SQ201610333423
【公開日】2016年10月26日
【申請日】2016年5月17日
【發(fā)明人】袁杰, 朱昀浩, 孫輝, 尤琦, 余雙春, 王學鼎, 陶超, 劉曉峻
【申請人】南京大學