利用U-Net網(wǎng)絡(luò)增強(qiáng)骨密度全波形反演

謝輝武 楊 艷

(武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院 武漢 430072)

0 引言

全波形反演(Full waveform inversion,FWI)是基于聲波方程約束的優(yōu)化問題。全波形反演最初由Tarantola 提出，主要是時(shí)域反演。當(dāng)然作為波動(dòng)方程約束的問題，也被推廣到了頻域。全波形反演最初是用在地質(zhì)勘探中，使用地震波來進(jìn)行反演，之后也用在了工程和醫(yī)學(xué)等領(lǐng)域。為了解決反演中的病態(tài)問題，目前最常用的是Tikhonnv正則化方法。使用更凸的目標(biāo)函數(shù)可能減小反演對(duì)準(zhǔn)確初始模型的需求，如Warner等[1]提出的自適應(yīng)反演，Yang[2]提出了使用新的目標(biāo)函數(shù)可以減小反演的病態(tài)問題。Bernard 等[3]嘗試將全波形反演應(yīng)用在骨密度測(cè)量中。但是己有的研究無法解決骨模型中數(shù)值差異大、更容易陷入局部極值的問題，反演結(jié)果中皮質(zhì)骨的聲速值只有2800 m/s，距離皮質(zhì)骨的實(shí)際聲速3500～4200 m/s還有很大的差距。

1 方法

本節(jié)將介紹如何將神經(jīng)網(wǎng)絡(luò)用于增強(qiáng)全波形反演。

1.1 逆時(shí)偏移成像

逆時(shí)偏移成像是一種基于波動(dòng)理論的成像方法。成像的步驟為：首先通過模擬零時(shí)刻到最大時(shí)刻的發(fā)射點(diǎn)波場(chǎng)傳播得到模型的發(fā)射點(diǎn)波場(chǎng)。然后將接收點(diǎn)接收到的波場(chǎng)從最大時(shí)刻逆時(shí)傳播到零時(shí)刻，得到接收點(diǎn)波場(chǎng)。最后按照一定的成像條件使用兩種波場(chǎng)進(jìn)行成像。該方法最早是由Whitemore[4]在1983年的SEG 會(huì)議上提出的。本文使用的成像條件是簡(jiǎn)單且常用的非時(shí)空移相關(guān)成像法，將震源波場(chǎng)與接收波場(chǎng)進(jìn)行零延時(shí)互相關(guān)就可得到偏移成像結(jié)果。

正向傳播使用有限差分法數(shù)值模擬波動(dòng)方程：

其中，ρ是密度，K=ρv2是彈性模量，p是波場(chǎng)值，s是發(fā)射點(diǎn)波場(chǎng)。

正向傳播可以得到接收數(shù)據(jù)。將接收數(shù)據(jù)逆時(shí)延拓，在合適的初始模型上傳播可以得到接收點(diǎn)波場(chǎng)R，將其與發(fā)射點(diǎn)波場(chǎng)執(zhí)行互相關(guān)運(yùn)算可以得到偏移成像。在地質(zhì)成像的過程中因?yàn)槟Ｐ洼^大，常常需要使用一些時(shí)間換空間的算法，如使用邊界存儲(chǔ)等方法減小存儲(chǔ)量[5]。與地質(zhì)成像不同，骨骼的模型大小較小，儲(chǔ)存波場(chǎng)所需的空間也不是很大。在偏移成像時(shí)可以直接存儲(chǔ)正向和反向傳播的波場(chǎng)，所以不需要使用邊界儲(chǔ)存等用時(shí)間換空間的算法。

因?yàn)殚L骨模型較為簡(jiǎn)單，在進(jìn)行逆時(shí)偏移時(shí)選取的初始背景模型可以是勻質(zhì)模型，速度假定為軟組織平均速度1500 m/s，密度為1000 kg/m3。

1.2 神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)

Ronneberger 等[6]提出了一種網(wǎng)絡(luò)結(jié)構(gòu)UNet。這種網(wǎng)絡(luò)結(jié)構(gòu)由下采樣部分和上采樣部分構(gòu)成，下采樣獲取信息，上采樣重構(gòu)圖像。為了在重構(gòu)時(shí)可以獲取原始的位置信息，采用了跳躍連接層。U-Net 的輸出可以用來進(jìn)行圖像分割，之前主要應(yīng)用于醫(yī)學(xué)圖像分割，近年來也有將其應(yīng)用到反演中。如Yang 等[7]提出的通過地震波形來生成地下速度模型，Zhu 等[8]使用U-Net 來提取走時(shí)。本文使用U-Net 提取逆時(shí)偏移成像的信息，重建出骨密度測(cè)量中皮質(zhì)骨，骨髓和軟組織的分布。

1.2.1 數(shù)據(jù)集

人體皮質(zhì)骨聲速為3500～4200 m/s，密度為1800 kg/m3左右。皮質(zhì)骨的厚度為4 mm 左右?？梢愿鶕?jù)這些數(shù)據(jù)生成大量有著不同組織厚度的模型，并給相同的組織在合理范圍內(nèi)提供不同的速度與密度值，使用有限差分法來正向模擬聲波的傳播，將接收點(diǎn)相應(yīng)位置的數(shù)據(jù)收集起來進(jìn)行逆時(shí)互相關(guān)偏移成像。將不同模型生成的偏移成像作為神經(jīng)網(wǎng)絡(luò)的輸入。本文使用1500 m/s、1000 kg/m3的背景值，使用225 種不同厚度和位置的簡(jiǎn)單平值骨骼模型，將其速度和密度賦值為不同的3 組，共得到675 種簡(jiǎn)單模型。再使用隨機(jī)方法生成背景上的675 種不同形狀和物性的模型，最后使用的樣本一共有1350 種。因?yàn)槟M使用是波動(dòng)方程的正向推理，當(dāng)差分精度合適時(shí)可以得到足夠精度的接收數(shù)據(jù)。超聲檢測(cè)裝置使用線性陣探頭，可以在手臂長骨的上下面各放置一個(gè)。假設(shè)探頭的間距是2 mm，在61 mm×51 mm 的模型中可以有62 個(gè)模擬源。表1介紹了一些正向模擬時(shí)的參數(shù)。

表1 超聲正向模擬時(shí)參數(shù)設(shè)置Table 1 Parameter setting during ultrasonic forward simulation

使用的超聲源為雷克子波，它可以看作是高斯函數(shù)的二次導(dǎo)數(shù)。它的表達(dá)式為公式(3)。波形和頻譜圖如圖1所示。

圖1 雷克子波波形與頻譜圖Fig.1 Ricker wavlet wavefomr and spectrum

1.2.2 網(wǎng)絡(luò)結(jié)構(gòu)

本文采用的U-Net 網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。輸入是1通道的偏移成像結(jié)果，下采樣過程每個(gè)單元由2個(gè)卷積層和其后的ReLU 激活層構(gòu)成，如果需要也可以在每個(gè)單元的輸入前加上一個(gè)參數(shù)標(biāo)準(zhǔn)化層(BN 層)和在卷積層之間添加Dropout 層。本文因?yàn)榫W(wǎng)絡(luò)任務(wù)簡(jiǎn)單，訓(xùn)練過程也沒有出現(xiàn)過擬合的現(xiàn)象，所以并沒有使用Dropout 層。每個(gè)單元之后經(jīng)過2×2 的最大池化層，張量尺寸減小。上采樣使用反卷積來完成。在反卷積的過程中，由于輸入張量的大小會(huì)出現(xiàn)奇數(shù)，反卷積后的張量大小與下采樣過程中的大小不匹配。而進(jìn)行跳躍連接需要兩端大小相同，這時(shí)可以使用簡(jiǎn)單的裁剪將上采樣過程中的張量大小裁剪為與下采樣過程中對(duì)應(yīng)深度的張量大小相同。經(jīng)過網(wǎng)絡(luò)后輸出為3 通道的分類結(jié)果，最后采用交叉熵作為損失函數(shù)。將之前用模擬方法生成的數(shù)據(jù)集送入網(wǎng)絡(luò)進(jìn)行訓(xùn)練。該網(wǎng)絡(luò)參數(shù)量為31, 054, 275。其中有31, 042, 499 個(gè)參數(shù)是需要訓(xùn)練的。該網(wǎng)絡(luò)的訓(xùn)練量并不大，甚至可以在一般的CPU上完成訓(xùn)練。

圖2 U-Net 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 Neural network structure diagram

該網(wǎng)絡(luò)訓(xùn)練的正確率和損失變化如圖3所示?？梢钥吹诫S著訓(xùn)練次數(shù)增加，正確率有了很大的提高，可以達(dá)到99%以上。

圖3 神經(jīng)網(wǎng)絡(luò)訓(xùn)練效果Fig.3 Neural network training results

表2展示了一些訓(xùn)練時(shí)的網(wǎng)絡(luò)超參數(shù)，U-Net的卷積核濾波器基礎(chǔ)數(shù)目為64，每一層的數(shù)目為64×2depth，depth表示網(wǎng)絡(luò)層在U-Net中的深度。訓(xùn)練時(shí)一共用了100個(gè)Epoch，學(xué)習(xí)率設(shè)為1.0×10?5。

表2 U-Net 超參數(shù)設(shè)置Table 2 U-Net hyperparameters

1.3 全波形反演的目標(biāo)函數(shù)

全波形反演的優(yōu)化目標(biāo)是讓檢測(cè)波場(chǎng)和通過預(yù)測(cè)模型模擬得到的接收波場(chǎng)均方差最小，因?yàn)樵搩?yōu)化問題在數(shù)學(xué)上是病態(tài)的，為了結(jié)果的穩(wěn)定需要添加上正則化，常用的有Tikhnonv 正則化[9]。本文提出可以將通過神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的模型m0加入目標(biāo)函數(shù)作為約束，看作是另外一種正則化。

式(4)中，pobs和pcal分別代表觀測(cè)的波場(chǎng)和模擬計(jì)算的波場(chǎng)，m是反演模型。而不同的α代表相應(yīng)的正則系數(shù)，用來控制該項(xiàng)對(duì)結(jié)果的影響力。因?yàn)閙0的數(shù)值并不一定準(zhǔn)確，所以需要將其系數(shù)調(diào)整為合適的值。為了增強(qiáng)分布約束，可以在反演開始的幾個(gè)迭代上使用較大的αcons，然后再將其調(diào)小，使波形匹配占據(jù)主要地位。

全波形反演可以使用各種優(yōu)化方法，如梯度下降法、高斯-牛頓法等。因?yàn)榉囱輰?duì)速度和密度的敏感性不同，實(shí)際反演中還需要一定的反演策略，比如說對(duì)密度使用更小的更新步長。本文綜合考慮反演準(zhǔn)確度和計(jì)算速度后決定使用高斯-牛頓法，其中觀測(cè)波場(chǎng)對(duì)模型的導(dǎo)數(shù)可以使用伴隨狀態(tài)法(格林函數(shù)法)求解[10]。

2 模擬實(shí)驗(yàn)

本文主要驗(yàn)證了以下幾種情況的對(duì)比。

2.1 正常數(shù)值的模型

如圖4所示，使用的模型長度為6 cm，寬度為5 cm，網(wǎng)絡(luò)間距為1 mm。皮質(zhì)骨聲速為4000 m/s，密度為1800 kg/m3。經(jīng)過偏移成像和神經(jīng)網(wǎng)絡(luò)后的預(yù)測(cè)結(jié)果如圖5所示，可以看到模型被準(zhǔn)確地分為軟組織、皮質(zhì)骨、骨骼3類?？芍獙?duì)于這種簡(jiǎn)單模型，神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)結(jié)果非常準(zhǔn)確。當(dāng)不加約束，只使用Tikhnonv 正則化時(shí)，結(jié)果如圖6(a)和圖6(b)所示；而添加約束后的結(jié)果如圖6(c)和圖6(d)所示。根據(jù)對(duì)比可以明顯看出當(dāng)約束與模型相同時(shí)對(duì)反演結(jié)果有很大的提升，中間的骨髓區(qū)域結(jié)果完整。接下來討論對(duì)于骨質(zhì)疏松，聲速和密度都有所下降的結(jié)果。

圖4 正常的骨骼模型Fig.4 Normal bono model

圖5 神經(jīng)網(wǎng)絡(luò)對(duì)圖3模型的分類預(yù)測(cè)Fig.5 The neural network predict result for Fig.3

圖6 正常模型兩種反演方法的結(jié)果對(duì)比Fig.6 The normal model results comparison between two methods

2.2 低密度模型

如圖7所示，使用的模型大小與之前相同，皮質(zhì)骨聲速為3500 m/s，密度為1700 kg/m3。經(jīng)過偏移成像和神經(jīng)網(wǎng)絡(luò)后的分類結(jié)果與圖5相同。然后執(zhí)行全波形反演，最后的結(jié)果如圖8所示。因?yàn)榇藭r(shí)模型中數(shù)值的變化沒有之前的大，沒有加約束的結(jié)果也還可以，不過在中間的骨髓區(qū)域加入了約束的結(jié)果形狀和數(shù)值都優(yōu)于未加約束的結(jié)果。對(duì)比圖8(a)與圖8(c)可以看出速度圖像的形狀更好，在骨髓區(qū)域的數(shù)值也優(yōu)于未加約束的結(jié)果。圖8(b)和圖8(c)的對(duì)比更加明顯，密度的結(jié)果圖形狀更好。圖9展示了速度和密度的皮質(zhì)骨橫向剖面和反演區(qū)域中間的縱向剖面結(jié)果對(duì)比。經(jīng)過圖片的對(duì)比可以看出約束后骨髓區(qū)域的數(shù)值更加優(yōu)秀。密度區(qū)域與假設(shè)的理論模型更加接近。而且進(jìn)行了約束的反演對(duì)于正則系數(shù)的健壯性有很大的提高。

圖7 低密度模型Fig.7 Low density bone model

圖8 低密度模型兩種反演方法的結(jié)果對(duì)比Fig.8 Results comparison between two methods for low density model

圖9 低密度模型兩種反演方法剖面對(duì)比Fig.9 Comparison of profile results of two inversion methods for low-density model

3 結(jié)論

與骨密度檢測(cè)的“金標(biāo)準(zhǔn)”雙能X 射線比較來說，聲學(xué)法儀器便宜，無輻射，有可能小型化后方便家庭進(jìn)行提前檢查。而定量聲學(xué)檢測(cè)法的準(zhǔn)確度還需要實(shí)踐的檢驗(yàn)。如果使用本文的方法可以對(duì)全波形骨密度測(cè)量的準(zhǔn)確度有一定的提高，還可以減小反演對(duì)初值的依賴。不過現(xiàn)有的骨密度檢測(cè)的標(biāo)準(zhǔn)是骨礦物質(zhì)密度，使用的是骨頭對(duì)光子的吸收能力相對(duì)值。使用聲學(xué)檢測(cè)需要在未來進(jìn)一步研究聲學(xué)的相對(duì)標(biāo)準(zhǔn)或它與光學(xué)法的結(jié)果對(duì)應(yīng)關(guān)系。