候英颯，孫正，2*，孫美晨

1.華北電力大學(xué)電子與通信工程系，保定 071003；2.華北電力大學(xué)河北省電力物聯(lián)網(wǎng)技術(shù)重點實驗室，保定 071003

0 引言

光聲層析成像（photoacoustic tomography，PAT）是一種新型的生物醫(yī)學(xué)成像技術(shù)。用短脈沖激光作為激勵源，生物組織吸收光能量并產(chǎn)生熱膨脹，進(jìn)而激發(fā)出寬帶超聲波（即光聲信號）。用超聲探測器收集光聲信號，通過求解聲學(xué)逆問題，獲得顯示成像目標(biāo)形態(tài)結(jié)構(gòu)的光吸收能量分布圖或者初始聲壓分布圖（Park等，2022）。在此基礎(chǔ)上求解光學(xué)逆問題，可以實現(xiàn)定量成像。

常規(guī)的PAT 圖像重建方法有反投影（back projection，BP）法（Xu 和Wang，2005）、時間反演（timereversal，TR）法（Sun等，2016）、基于傅里葉變換的重建（Fourier transform-based reconstruction，F(xiàn)TR）方法、延遲求和（delay and sum，DAS）法以及迭代重建技術(shù)（iterative reconstruction technology，IRT）（Rosenthal 等，2013；Lutzweiler 和Razansky，2013）等。為了簡化問題，這些方法通常對成像場景進(jìn)行理想化的假設(shè)。例如：超聲波在介質(zhì)中以恒定的速度進(jìn)行無損傳播，超聲探測器被簡化成理想的點探測器，能夠以足夠高的時間和空間采樣率在全視角內(nèi)采集到完備光聲信號等。然而，在實際應(yīng)用中，上述條件往往難以滿足。例如：超聲波在不同成分的介質(zhì)中的傳播速度通常存在較大差異，在重建圖像的過程中，恒定聲速（speed of sound，SoS）的假設(shè)會導(dǎo)致圖像質(zhì)量明顯下降。解決該問題的方法有兩種，一是在圖像重建中引入有關(guān)聲速分布的先驗知識，補償由于聲速變化引起的測量光聲場的像差，改善圖像的聚焦效果。此類方法需要已知聲速分布的先驗知識，然而多數(shù)情況下，在進(jìn)行光聲成像前無法獲得該先驗知識（Huang 等，2013；Haltmeier 和Nguyen，2017）。另一種方法是根據(jù)光聲信號，同時重建聲速分布和初始聲壓或光吸收能量密度（absorbed optical energy density，AOED）分布，無需有關(guān)聲速的先驗知識（Zhang等，2008；Matthews 等，2018）。目前的主流方法是采用基于模型的迭代框架，即選擇適當(dāng)?shù)膬?yōu)化策略，迭代求解前向成像模型的逆問題（Huang 等，2016；Sun 和Sun，2021）。此類方法的重建精度很大程度上依賴于有關(guān)成像目標(biāo)的先驗假設(shè)模型。為了降低問題的病態(tài)程度，需要合理選擇正則化方法及其參數(shù)。此外，求解過程非常耗時。

近年來，深度學(xué)習(xí)已成為醫(yī)學(xué)成像領(lǐng)域的研究熱點（唐朝生 等，2021；李書林 等，2022）。為了減輕逆問題求解中的病態(tài)，使用有監(jiān)督學(xué)習(xí)進(jìn)行圖像重建也成為該領(lǐng)域的最新技術(shù)，主要包括非迭代學(xué)習(xí)法和學(xué)習(xí)迭代法兩類（Farnia 等，2020；Guo 等，2021）。按照深度神經(jīng)網(wǎng)絡(luò)輸入和輸出的不同，非迭代學(xué)習(xí)法分為信號增強、圖像增強和優(yōu)化重建3 類（Lan 等，2019b，2021）。信號增強是采用“信號”到“信號”的轉(zhuǎn)換思路，利用深度學(xué)習(xí)技術(shù)提高原始光聲信號的信噪比，然后采用常規(guī)方法重建圖像（Awasthi 等，2020）。圖像增強采用“圖像”到“圖像”的轉(zhuǎn)換思路，首先采用常規(guī)方法重建初始圖像，然后用神經(jīng)網(wǎng)絡(luò)對初始圖像進(jìn)行優(yōu)化（Vu 等，2020；Sun和Yan，2020）。這類方法的成像質(zhì)量受到初始重建精度的限制，而采用常規(guī)方法重建圖像時，可能無法恢復(fù)原始信號中隱含的細(xì)節(jié)信息。特別是當(dāng)探測器在有限角度內(nèi)采集到不完備的光聲信號時，采用該方法很難完整地恢復(fù)出掃描盲區(qū)中的結(jié)構(gòu)。優(yōu)化重建方法是利用神經(jīng)網(wǎng)絡(luò)直接學(xué)習(xí)“信號”到“圖像”的映射（Guan 等，2020）。網(wǎng)絡(luò)的輸入是原始信號或預(yù)處理后的信號，輸出是高質(zhì)量圖像（Lan 等，2019a；Kim 等，2020）。此類方法需要利用全連接或大卷積核的卷積神經(jīng)網(wǎng)絡(luò)解決輸入信號矩陣的非對稱性問題，難以獲得復(fù)雜目標(biāo)的準(zhǔn)確重建（Guo 等，2021；Lan等，2020）。

學(xué)習(xí)迭代法是將迭代重建技術(shù)與深度學(xué)習(xí)相結(jié)合，神經(jīng)網(wǎng)絡(luò)從訓(xùn)練數(shù)據(jù)中學(xué)習(xí)得到先驗約束，用于優(yōu)化迭 代過程（Adler 和?ktem，2017；Shan 等，2019）。例如：利用神經(jīng)網(wǎng)絡(luò)訓(xùn)練正則化工具（Antholzer等，2019；Li等，2020），或者通過深度梯度下降（deep gradient descent，DGD）進(jìn)行迭代更新（Boink 等，2020；Sun 等，2021），確保每次迭代都收斂于最優(yōu)值。相較于傳統(tǒng)的迭代重建技術(shù)，學(xué)習(xí)迭代法可以通過更少的迭代獲得更高質(zhì)量的圖像。但是由于神經(jīng)網(wǎng)絡(luò)需要從訓(xùn)練數(shù)據(jù)中學(xué)習(xí)先驗約束，因而此類方法的時間成本較高，訓(xùn)練數(shù)據(jù)的質(zhì)量也會對成像結(jié)果產(chǎn)生影響。

本文搭建一個基于學(xué)習(xí)迭代策略的聯(lián)合迭代重建網(wǎng)絡(luò)（joint iterative reconstruction network，JIRNet），根據(jù)有限角度稀疏光聲信號，同時重建AOED分布圖和聲速分布圖，解決由于不完備測量數(shù)據(jù)和不均勻聲速導(dǎo)致的圖像質(zhì)量下降問題。采用仿真、仿體和在體數(shù)據(jù)訓(xùn)練和測試網(wǎng)絡(luò)。通過與傳統(tǒng)非學(xué)習(xí)方法、非迭代學(xué)習(xí)方法和基于DGD 的學(xué)習(xí)迭代方法比較，證明在稀疏測量和組織聲學(xué)特性不均勻情況下，該方法重建光聲圖像的優(yōu)勢。

1 基于迭代梯度網(wǎng)絡(luò)的聯(lián)合重建方法

基于學(xué)習(xí)迭代策略搭建卷積神經(jīng)網(wǎng)絡(luò)模型，將AOED 和聲速的梯度下降信息融入到網(wǎng)絡(luò)的訓(xùn)練中，通過反向傳播梯度下降求解非線性最小二乘問題，實現(xiàn)從不完備光聲信號到高質(zhì)量AOED 分布圖和聲速分布圖的映射。

1.1 成像模型

短脈沖激光照射組織、成像平面內(nèi)的AOED 與入射光輻射通量和組織光吸收系數(shù)之間的關(guān)系為

式中，A(r)和μa(r)分別是成像平面中位置r處的AOED 和光吸收系數(shù)，Φ(r)是光輻射通量。組織中的光吸收體基于光聲效應(yīng)產(chǎn)生的初始聲壓為

式中，p0(r)是初始聲壓，Γ是熱膨脹系數(shù)。光聲信號在組織中的傳播用如下波動方程描述，具體為

式中，p(r，t)是位置r在時刻t的聲壓（pressure），I(t)是入射激光脈沖能量函數(shù)，c(r)是聲速（speed of sound，SoS），βe和Cp分別是組織的等壓膨脹系數(shù)和一定壓力下的比熱容。

從AOED到產(chǎn)生聲壓的過程用算子表示為

式中，p是到達(dá)探測器表面的聲壓矩陣，A是AOED矩陣，H(c)是與聲速c有關(guān)的算子。

求解如下非線性最小二乘問題，根據(jù)聲壓信號同時重建AOED和聲速分布，具體為

1.2 JIR-Net的結(jié)構(gòu)

如圖1 所示，JIR-Net 的輸入是在成像平面內(nèi)采集的聲壓信號矩陣和預(yù)先設(shè)置的常數(shù)聲速，輸出是AOED分布圖和聲速分布圖。

圖1 本文方法流程圖Fig.1 Flow chart of the proposed method

如圖2 所示，網(wǎng)絡(luò)由4 個結(jié)構(gòu)單元組成，每個單元的輸入是前一個單元輸出的AOED 分布圖、歸一化聲速分布圖以及二者的梯度。一個單元由特征提取、特征融合和重建3個模塊組成，如圖3所示。

圖2 JIR-Net結(jié)構(gòu)示意圖Fig.2 Architecture of JIR-Net

圖3 結(jié)構(gòu)單元示意圖Fig.3 Schematic of a structural unit

特征提取模塊分別對4 個輸入進(jìn)行兩次濾波核尺寸是5 × 5 ×S、步長為1的same卷積，S是當(dāng)前層輸入的特征圖數(shù)量，四路輸入的特征通道數(shù)均為16。

特征融合模塊的作用是：將AOED 分布圖及其梯度進(jìn)行像素值求和，將歸一化聲速分布圖及其梯度進(jìn)行像素值求和，然后將得到的兩個特征圖沿通道維度拼接起來，并進(jìn)行3次卷積操作。

重建模塊通過兩次卷積操作，得到AOED 和聲速的更新值，并通過遠(yuǎn)跳連接，將該單元的輸入與當(dāng)前層輸出的更新值進(jìn)行像素值求和，最后的ReLU激活函數(shù)保證更新后的AOED 和速度分布是非負(fù)分布。

第i個結(jié)構(gòu)單元輸入的AOED梯度為

式中，pk是訓(xùn)練集中第k個樣本輸入的聲壓信號測量值，Ai-1，k和ci-1，k分別是經(jīng)過i-1次迭代后第k個樣本輸出的AOED 估計值和反歸一化后的聲速估計值。為方便起見，文中將?d1，i[pk，H(ci-1，k)Ai-1，k]簡寫為?d1，i(?)。

聲速梯度利用伴隨波動方程求得。伴隨波動方程的表達(dá)式（Huang等，2016）為

式中，q(r，t)是伴隨波，s(r，t)是伴隨源，具體為

式中，rl是探測器的第l個測量位置，L是探測器測量位置的總數(shù)，(rl，t)是探測器在位置rl處時刻t測得的聲壓，p(rl，t)是聲壓理論值，由式（3）求得，δ(r-rl)是狄拉克函數(shù)。式（7）的初始條件為

式中，T是采集信號的時間長度。根據(jù)式（3）和式（7）得到聲速的梯度，具體為

根據(jù)式（10）計算第i個結(jié)構(gòu)單元中聲速的梯度?d2，i[pk，Η(ci-1，k)Ai-1，k]，對其歸 一化后 表示為?d2，i，N[pk，Η(ci-1，k)Ai-1，k]，簡寫為?d2，i，N(?)。對第i個單元輸入的聲速分布ci-1歸一化，記為ci-1，N。

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

1.3.1 計算機仿真數(shù)據(jù)集

建立包含不同組織成分的計算機仿真模型，其中4 個模型的幾何結(jié)構(gòu)如圖4 所示，并設(shè)置其組織特性參數(shù)，如表1 所示。采用蒙特卡羅方法模擬入射光在組織中的傳輸，得到成像平面中的光輻射通量（Wang 等，1995）。進(jìn)而根據(jù)式（1），由光輻射通量和光吸收系數(shù)得到AOED（absorbed optical energy density）。采用MATLAB k-Wave 工具包仿真得到探測器接收的光聲信號（Treeby 和Cox，2010），其中加入30 dB 的高斯隨機噪聲。探測器采樣頻率和中心頻率分別設(shè)置為250 MHz 和5 MHz。將不同稀疏度的仿真光聲信號和常數(shù)聲速共同作為樣本輸入，將仿真AOED 分布圖和預(yù)設(shè)的聲速分布圖作為期望輸出，構(gòu)建數(shù)據(jù)集。部分樣本示例如圖4 所示。

表1 仿真模型的組織特性參數(shù)Table 1 Characteristic parameters of simulated phantoms

圖4 仿真數(shù)據(jù)集中的樣本圖像示例（Amp·表示幅值A(chǔ)mplitude，a.u.表示任意單位arbitray units）Fig.4 Examples of samples in simulation datasets

構(gòu)建了3 個由不同稀疏度仿真信號和不同常數(shù)聲速作為樣本輸入的仿真數(shù)據(jù)集，如表2 所示。用“（視角，采樣率）”表示數(shù)據(jù)稀疏度，如（180°，0.7）表示掃描角度是180°、采樣率為0.7。將每個數(shù)據(jù)集中的樣本隨機打亂，按照6∶3∶1 的比例劃分為訓(xùn)練集、驗證集和測試集，利用驗證集選擇最佳網(wǎng)絡(luò)模型。

表2 仿真數(shù)據(jù)集的構(gòu)成Table 2 Description of simulation datasets

1.3.2 仿體數(shù)據(jù)集

制作了14 個包含不同嵌入物的柱形仿體，不同材料的聲速和光吸收系數(shù)如表3 所列，制作細(xì)節(jié)如表4所列，部分仿體實物照片如圖5所示。

表3 制作仿體的不同材料的聲速和光吸收系數(shù)Table 3 SoS and absorption coefficients of different materials used to fabricate phantoms

表4 仿體的構(gòu)成Table 4 Composition of phantoms

圖5 仿體實物照片F(xiàn)ig.5 Photos of phantoms

每個仿體的總長均為80 mm，橫截面直徑為20 mm。用于采集仿體光聲信號的成像系統(tǒng)采用Nd：YAG 激光器輸出的全波波長（680～900 nm）可調(diào)諧激光作為光源，脈沖重復(fù)率是10 Hz，脈沖寬度是7 ns，最大脈沖能量是120 mJ。超聲探測裝置由256個中心頻率為5 MHz、帶寬為60%的聚焦超聲換能器陣元，排列成覆蓋角度為270°的環(huán)形陣列。對每個仿體的掃描長度約為10 mm，對4個仿體（即表3中編號為I-1、I-3、II-8 和III-11 的仿體）得到的28 個掃描切片構(gòu)成測試集，對其余10個仿體得到的69個掃描切片構(gòu)成訓(xùn)練集。樣本輸入的常數(shù)聲速均為1 500 m·s-1。為了避免產(chǎn)生過擬合，采用數(shù)據(jù)增強技術(shù)（包括隨機翻轉(zhuǎn)、旋轉(zhuǎn)和平移等）對訓(xùn)練集進(jìn)行擴(kuò)增，最終的訓(xùn)練集中包含4 140個樣本。

1.3.3 在體數(shù)據(jù)集

在體數(shù)據(jù)集由對活體小鼠進(jìn)行全身掃描采集的光聲信號構(gòu)成。成像系統(tǒng)采用Nd：YAG 激光器輸出的715 nm、730 nm、760 nm、800 nm、850 nm 可調(diào)諧激光作為光源，在730 nm 處的最大入射脈沖能量約為70 mJ，脈沖重復(fù)率是10 Hz，脈沖寬度是8 ns。超聲探測裝置是由128個中心頻率為5 MHz、帶寬為60%的聚焦超聲換能器陣元構(gòu)成的環(huán)形陣列。利用前述每種波長的激光對俯臥位小鼠進(jìn)行全身掃描，共獲得244 個切片，其中200 個切片數(shù)據(jù)用做構(gòu)建訓(xùn)練集，44個切片構(gòu)建測試集，樣本輸入的常數(shù)聲速均設(shè)置為1 500 m·s-1。增強后的訓(xùn)練集中包含4 392個樣本。

1.4 訓(xùn)練網(wǎng)絡(luò)

采用自適應(yīng)矩估計（adaptive moment estimation，Adam）梯度下降法（Kingma 和Ba，2017）沿各層參數(shù)梯度方向逐個單元進(jìn)行參數(shù)優(yōu)化。第i個單元的損失函數(shù)為

式中，L1，i和L2，i分別是第i個單元重建AOED分布圖（即任務(wù)1）和聲速分布圖（即任務(wù)2）的損失函數(shù)，具體為

式中，M是樣本總數(shù)，Atrue，k和ctrue，k，N分別是第k個樣本期望輸出的AOED 分布圖和歸一化聲速分布圖，Ai，k和ci，k，N分別是 第i個單元 輸出的 第k個樣本 的AOED估計值和歸一化聲速估計值，具體為

設(shè)置學(xué)習(xí)率為0.001，批處理大小為batch size=64，epoch=200。損失函數(shù)曲線如圖6 所示，可見隨著epoch 的增加，損失函數(shù)值逐漸減小，當(dāng)epoch >200時，變化趨于平緩。

圖6 采用稀疏度為（180°，0.9）的仿真數(shù)據(jù)集訓(xùn)練JIR-Net時，4個結(jié)構(gòu)單元的損失函數(shù)曲線Fig.6 Loss of the four structural units when JIR-Net is trained using the simulation dataset with a sparsity of（180°，0.9）（（a）Unit 1；（b）Unit 2；（c）Unit 3；（d）Unit 4）

2 方法的性能評價

分別采用計算機仿真、仿體和在體成像數(shù)據(jù)定量評價JIR-Net 根據(jù)不完備光聲信號重建圖像的性能。仿真和仿體實驗的目的是在已知目標(biāo)特性參數(shù)真實值的情況下，定量評估重建精度，并討論算法中的相關(guān)參數(shù)對其性能的影響。在體實驗的目的是驗證JIR-Net在真實場景中應(yīng)用的可行性。

將JIR-Net 與經(jīng)典TR 法（Sun 等，2016）、交替優(yōu)化（alternative optimization，AO）法（Sun 和Sun，2021）、U-Net 后處理法（Sun 和Yan，2020）（簡稱U-Net）和DGD 網(wǎng)絡(luò)（Sun 等，2021）進(jìn)行比較。AO 法是通過交替優(yōu)化迭代求解信號測量值和理論值之間誤差最小的非線性最小二乘問題。U-Net 后處理法是將采用TR 法重建的低質(zhì)量圖像作為U-Net 的輸入，最終輸出優(yōu)化后的圖像。DGD 網(wǎng)絡(luò)共包含5 個結(jié)構(gòu)單元。采用相同的數(shù)據(jù)集訓(xùn)練JIR-Net、DGD 和U-Net，JIR-Net 和DGD 的訓(xùn)練參數(shù)為batch size=64和epoch=200。U-Net 的訓(xùn)練參數(shù)為batch size=1和epoch=200。U-Net、DGD 和JIR-Net 的訓(xùn)練時間分別約為2.4 h、6.5 h和17 h。

在仿真實驗中，以仿真AOED 分布圖和預(yù)設(shè)的聲速分布圖作為標(biāo)準(zhǔn)圖像，采用峰值信噪比（peak signal noise ratio，PSNR）和結(jié)構(gòu)相似度（structural similarity，SSIM）作為定量指標(biāo)，評價重建圖像質(zhì)量。重建圖像和標(biāo)準(zhǔn)圖像之間的SSIM 越大，PSNR 越高，表明重建質(zhì)量越佳。仿體和在體實驗中，由于無法獲得AOED 的真實值，因此采用對比度（contrast ratio，CR）和對比 噪聲比（contrast-to-noise ratio，CNR）（單位：dB）作為量化指標(biāo)評價AOED 圖像質(zhì)量。CR衡量圖像灰度反差的程度，CNR是圖像中有用信號區(qū)域和背景區(qū)域之間對比度的度量。CR 和CNR 越高，表示重建圖像質(zhì)量越好。此外，仿體實驗中，不同材料的聲速真實值是已知的（如表4 所列），因此采用PSNR和SSIM定量評價聲速分布的重建精度。在體實驗中，僅通過光聲掃描無法獲得聲速的真實值，因此采用視覺評估定性評價聲速分布的重建質(zhì)量。

實驗中，仿真、TR 和AO 法實現(xiàn)的編程環(huán)境是MATLAB（R2016a，the Math Works，Inc.，Natick，Massachusetts），計算機配置為12th Gen Intel（R）Core（TM）i5-12500H CPU、16 GB RAM和Windows 11 64 位操作系統(tǒng)。神經(jīng)網(wǎng)絡(luò)的搭建、訓(xùn)練和測試的運行環(huán)境為Ubuntu 16.04，所用GPU是NVIDIA公司生產(chǎn)的GeForce 3090Ti，顯存為16 GB。操作系統(tǒng)是Ubuntu 16.04，編程語言是Python 3.7，深度學(xué)習(xí)框架是TensorFlow 2.6.0和Keras 2.0。

3 實驗結(jié)果與分析

分別對計算機仿真、仿體和在體3 種成像數(shù)據(jù)對網(wǎng)絡(luò)性能進(jìn)行測試，并對實驗結(jié)果進(jìn)行分析。

3.1 仿真實驗結(jié)果

對于不同稀疏度的仿真光聲信號，分別采用不同方法重建圖像的結(jié)果如圖7所示。由圖7可知，在全角度掃描時，信號的采樣率越低，AO和TR重建圖像的整體亮度也越低，這是由于在稀疏采樣時將未探測到的聲壓設(shè)為零造成的。采用深度學(xué)習(xí)方法重建圖像的視覺效果和評價指標(biāo)明顯優(yōu)于AO 和TR法。當(dāng)采樣率較高時，3個網(wǎng)絡(luò)都可以重建出較高質(zhì)量的圖像。但是當(dāng)采樣率較低時（如0.3 和0.1），JIR-Net 重建圖像的質(zhì)量則明顯優(yōu)于DGD 和U-Net。而且JIR-Net重建圖像的質(zhì)量并沒有隨著信號采樣率減小而明顯降低，這表明經(jīng)過訓(xùn)練的JIR-Net對測量信號的稀疏度具有一定的魯棒性。當(dāng)進(jìn)行有限角度掃描時，TR重建質(zhì)量隨著掃描角度的減小而降低，角度越小，盲區(qū)越大，圖像中目標(biāo)結(jié)構(gòu)缺失的部分就越多。與TR相比，雖然AO和U-Net可以在一定程度上改善圖像質(zhì)量，但是隨著掃描角度的減小，其能夠補全的盲區(qū)也越少。DGD雖然可以補全盲區(qū)，但由于假設(shè)聲速恒定，與真實聲速分布并不相符，因此存在目標(biāo)嚴(yán)重變形的問題。而采用JIR-Net重建的圖像，無論視覺效果還是評價指標(biāo)均優(yōu)于其他方法，且掃描角度越小，優(yōu)勢越明顯。180°掃描時，JIR-Net重建圖像的SSIM 比DGD 高約8.4%，比U-Net 高 約18.3%，PSNR比DGD高約17.4%，比U-Net高約74.7%。

圖7 根據(jù)不同稀疏度的仿真信號重建圖像的結(jié)果Fig.7 Results of image reconstruction from simulated signals with different sparsity（（a）AOED images；（b）speed of sound images；（c）AOED metrics；（d）SoS metrics）

3.2 仿體實驗結(jié)果

對仿體的圖像重建結(jié)果如圖8 所示，與實物照片對比，采用5 種方法重建的圖像都能夠可靠地顯示嵌入目標(biāo)的位置。但是采用AO 和TR 法重建圖像的整體亮度明顯低于3 種深度學(xué)習(xí)方法重建的圖像。在DGD、U-Net、AO和TR重建圖像中，可以清晰地觀察到嵌入物周圍的偽影。而JIR-Net 重建的圖像具有最佳的視覺效果和最高的對比度。

圖8 仿體的圖像重建結(jié)果Fig.8 Results of image reconstruction of phantoms（（a）AOED images；（b）SoS images；（c）AOED metrics；（d）SoS metrics）

3.3 在體實驗結(jié)果

圖9 是包含復(fù)雜結(jié)構(gòu)的活體小鼠4 個頭部和胸腹部切片的成像結(jié)果。由圖9 可知，與仿真和仿體實驗結(jié)果類似，采用AO 和TR 法重建圖像的整體亮度以及橫截面邊緣的對比度都明顯低于深度學(xué)習(xí)方法重建的圖像。在JIR-Net 重建圖像中，箭頭標(biāo)示細(xì)節(jié)結(jié)構(gòu)的對比度明顯高于采用其他方法重建圖像的對比度。這證明了在重建復(fù)雜目標(biāo)方面，JIR-Net 優(yōu)于非學(xué)習(xí)方法以及DGD 和U-Net方法。

圖9 在體小鼠實驗中的圖像重建結(jié)果Fig.9 Results of image reconstruction for in vivo mice（（a）AOED images；（b）SoS images；（c）AOED metrics）

4 討論

結(jié)合仿真實驗結(jié)果，分析討論初始聲速以及網(wǎng)絡(luò)結(jié)構(gòu)單元個數(shù)對重建圖像質(zhì)量的影響。并對所提方法相對于其他成像方案的優(yōu)勢和不足進(jìn)行討論。

4.1 初始聲速對重建圖像質(zhì)量的影響

為了分析初始聲速對JIR-Net 性能的影響，在其他條件不變的情況下，分別設(shè)置初始聲速為1 404 m·s-1、1 449 m·s-1、1 500 m·s-1和1 716 m·s-1。對全角度密集采樣仿真光聲信號進(jìn)行圖像重建，結(jié)果如圖10所示。

由圖10 可知，當(dāng)采用不同的初始聲速時，JIRNet重建圖像的視覺效果和評價指標(biāo)變化并不顯著。由表1 可知，本文建立仿真模型時，將不同組織成分的平均聲速設(shè)定在1 540～1 650 m·s-1范圍內(nèi)。當(dāng)初始聲速與實際聲速比較接近時（如1 500 m·s-1），JIR-Net、DGD、U-Net 和AO 重建圖像的質(zhì)量差別不大，且明顯優(yōu)于TR法。隨著初始聲速與實際聲速的差別增大，JIR-Net 相較于其他方法表現(xiàn)出明顯的優(yōu)勢。例如，當(dāng)初始聲速為1 716 m·s-1時，遠(yuǎn)大于光聲信號傳播過程中的實際聲速。此時，TR 重建中聲波的反演時間變短，導(dǎo)致重建出的目標(biāo)區(qū)域不完全，且存在嚴(yán)重畸變。將TR 重建圖像作為U-Net 的輸入時，輸出的圖像質(zhì)量也沒有得到明顯改善。對于DGD 網(wǎng)絡(luò)，由于所設(shè)初始聲速與真實值差別較大，每次求解梯度時均會產(chǎn)生誤差，逐次迭代后誤差累積，導(dǎo)致圖像中的畸變嚴(yán)重。此外，由圖10（b）（d）可知，不同初始聲速條件下，JIR-Net 重建聲速的精度優(yōu)于AO法，SSIM和PSNR可分別提高約34.4%和22.6%。

4.2 結(jié)構(gòu)單元數(shù)對網(wǎng)絡(luò)性能的影響

為了分析JIR-Net 中的結(jié)構(gòu)單元數(shù)對網(wǎng)絡(luò)性能的影響，在其他條件和參數(shù)不變的情況下，分別將結(jié)構(gòu)單元數(shù)設(shè)定為1、2、3、4 和5。對稀疏度為（180°，0.9）的仿真光聲信號進(jìn)行圖像重建，評價指標(biāo)如圖11 所示。由圖11 可知，結(jié)構(gòu)單元數(shù)越少，重建圖像質(zhì)量越低。增加結(jié)構(gòu)單元數(shù)雖然會改善圖像質(zhì)量，但也會導(dǎo)致網(wǎng)絡(luò)訓(xùn)練時間明顯延長，如表5 所示。當(dāng)單元數(shù)大于4 時，重建圖像的評價指標(biāo)變化趨于平緩，表明此時再增加結(jié)構(gòu)單元也不會顯著改善圖像質(zhì)量，但是訓(xùn)練時間會大幅延長。因此，本文最終采用4 個結(jié)構(gòu)單元構(gòu)成JIR-Net。

表5 JIR-Net由不同數(shù)量結(jié)構(gòu)單元組成時的訓(xùn)練時長Table 5 Time cost in training JIR-Net consisting of different number of structural units

圖11 當(dāng)JIR-Net包含不同數(shù)量的結(jié)構(gòu)單元時，根據(jù)稀疏度為（180°，0.9）的仿真信號重建圖像的評價指標(biāo)Fig.11 Evaluation metrics of images reconstructed from the simulated signal with the sparsity of（180°，0.9）when JIR-Net consists of different numbers of structural units

4.3 所提方法相對于其他成像方案的優(yōu)勢和不足

JIR-Net 采用學(xué)習(xí)迭代框架實現(xiàn)有限角度稀疏測量條件下的高質(zhì)量圖像重建。與非迭代學(xué)習(xí)方法相比，這種方案以犧牲時間成本為代價獲得了高質(zhì)量的成像結(jié)果。訓(xùn)練JIR-Net 單個結(jié)構(gòu)單元的時間約為DGD 的2.6 倍，U-Net 的7 倍。采用訓(xùn)練后的JIR-Net、DGD 和U-Net 重建一幀256 × 256 像素的圖像，所需的時間分別為18.4 s、14 s和9 s。時間成本的主要差別在于JIR-Net 和DGD 需要計算AOED 和聲速的梯度，而U-Net不用。

此外，通過同時重建AOED 分布圖和聲速分布圖，本文方法解決了由于不均勻聲速所致圖像質(zhì)量下降的問題。該方法不需要任何有關(guān)成像目標(biāo)的先驗知識。與基于模型的非學(xué)習(xí)方法相比，成像質(zhì)量和精度不依賴于前向模型，無需對成像目標(biāo)和場景進(jìn)行理想化的假設(shè)，可直接應(yīng)用于不同的掃描幾何和成像場景。

5 結(jié)論

本文提出一種聯(lián)合迭代重建網(wǎng)絡(luò)，通過將AOED 分布和聲速分布的梯度下降信息融合到網(wǎng)絡(luò)訓(xùn)練中，實現(xiàn)利用不完備光聲信號重建高質(zhì)量圖像的目的。仿真、仿體和在體實驗結(jié)果表明，采用該網(wǎng)絡(luò)重建的AOED 分布圖，其視覺效果和評價指標(biāo)均優(yōu)于經(jīng)典TR 方法、AO 方法、基于U-Net 的非迭代學(xué)習(xí)方法和基于DGD 的學(xué)習(xí)迭代方法。而且光聲信號測量數(shù)據(jù)的稀疏度越低，優(yōu)勢越明顯。討論結(jié)果表明，所提方法的重建精度對初始聲速不敏感，具有較好的穩(wěn)定性。網(wǎng)絡(luò)結(jié)構(gòu)單元個數(shù)對重建精度有一定的影響，雖然增加結(jié)構(gòu)單元可以提高重建質(zhì)量，但是也會進(jìn)一步增加時間成本，需要在圖像質(zhì)量和時間成本之間做出權(quán)衡。

在未來的工作中，將從以下3 個方面改善網(wǎng)絡(luò)性能：1）構(gòu)建規(guī)模更大、質(zhì)量更高的數(shù)據(jù)集，盡可能多地包含可能觀察到的圖像特征，并使仿真訓(xùn)練集適應(yīng)使用真實光傳播模型生成的圖像。2）進(jìn)一步優(yōu)化網(wǎng)絡(luò)結(jié)構(gòu)并改進(jìn)訓(xùn)練策略，以更快的速度重建更高質(zhì)量的圖像。3）改變網(wǎng)絡(luò)輸入信號的類型，對原始光聲信號測量數(shù)據(jù)進(jìn)行預(yù)處理后再送入網(wǎng)絡(luò)，提高網(wǎng)絡(luò)的學(xué)習(xí)效率。

致謝：本項工作中活體小鼠光聲成像數(shù)據(jù)的采集得到了中國醫(yī)學(xué)科學(xué)院藥物研究所的支持，仿體成像數(shù)據(jù)的采集得到了北京化工大學(xué)的支持，在此表示感謝。