簡(jiǎn)介MATLAB 光聲仿真工具箱 K-Wave 1.2.1 完整資源包面向生物醫(yī)學(xué)成像、材料科學(xué)和水聲探測(cè)等領(lǐng)域的科研人員與工程師用于光聲效應(yīng)數(shù)值模擬、聲波傳播計(jì)算及光聲圖像重建。壓縮包共 619 個(gè)文件、6.51MB包含 201 個(gè) m 源程序與示例腳本、174 個(gè) html 幫助文檔、233 個(gè) png 演示圖像另有 gif、txt、xml、mat、css 等輔助文件目錄結(jié)構(gòu)清晰便于在 MATLAB 中直接加載并按文檔示例復(fù)現(xiàn)實(shí)驗(yàn)。工具箱基于有限差分法覆蓋光聲成像模擬、聲傳播模擬、數(shù)據(jù)采集與重建、可視化等完整流程支持點(diǎn)、線、面、體等多種光源及 PML 邊界條件可靈活設(shè)置聲速、密度、吸收系數(shù)等聲學(xué)參數(shù)。對(duì)初學(xué)者而言可從 html 文檔快速了解函數(shù)用法再對(duì)照 m 腳本修改光源與介質(zhì)參數(shù)觀察 png 輸出結(jié)果檢驗(yàn)聲場(chǎng)和重建效果從而縮短上手周期并支撐實(shí)驗(yàn)設(shè)計(jì)、參數(shù)測(cè)試與算法驗(yàn)證。已有 2385 人學(xué)習(xí)是深入理解光聲成像物理過程、開展相關(guān)仿真的實(shí)用工具。1. 光聲仿真到底在仿什么K-Wave-toolbox 1.2.1 能替你省下哪幾步做光聲成像的人第一年基本都耗在“聲”上。光學(xué)部分有現(xiàn)成的蒙特卡洛、擴(kuò)散方程求解器可光吸收完變成熱、熱變成壓力波、壓力波在組織里傳播到探頭這一段聲學(xué)傳播很少有人講清楚。K-Wave-toolbox 1.2.1 就是干這個(gè)的——它求解的是均勻或非均勻媒質(zhì)中的聲波方程輸入初始?jí)毫Ψ植紁0輸出傳感器位置上收到的時(shí)域壓力信號(hào)。換句話說你只需要算好“光在哪被吸收了多少”剩下的聲傳播、反射、折射、衰減它全包了。這套工具箱對(duì)兩類人價(jià)值最大一是做光聲成像重建算法的人需要大量仿真數(shù)據(jù)驗(yàn)證反演算法手工解析解根本算不了非均勻媒質(zhì)二是做系統(tǒng)設(shè)計(jì)的人要預(yù)估探頭布局、孔徑大小、中心頻率對(duì)圖像的影響。它用 k-space 偽譜法做時(shí)間推進(jìn)同樣的網(wǎng)格規(guī)模和精度要求下比有限元快一到兩個(gè)數(shù)量級(jí)。這篇筆記按“原理 → 安裝 → 最小算例 → 參數(shù)調(diào)優(yōu) → 踩坑 → 進(jìn)階驗(yàn)證”的順序講一遍所有命令我都按 1.2.1 版實(shí)測(cè)過老版本升級(jí)上來的人重點(diǎn)看第三章和第五章。2. k-space 偽譜法與 k-Wave 的三大數(shù)據(jù)結(jié)構(gòu)先搞清楚再動(dòng)手2.1 k-space 偽譜法為什么是光聲仿真的主流選擇光聲仿真本質(zhì)上是求解如下形式的聲波方程?2p/?t2 c2?2p 源項(xiàng)有限元法把空間離散成網(wǎng)格每個(gè)時(shí)間步要解一個(gè)大型稀疏線性方程組三維情況下自由度輕松上千萬內(nèi)存和時(shí)間都吃不消。有限差分法快一些但數(shù)值色散嚴(yán)重——波傳播幾十個(gè)網(wǎng)格后波形就畸變了。k-Wave 用偽譜法空間導(dǎo)數(shù)通過傅里葉變換在頻域里計(jì)算精度在單個(gè)網(wǎng)格內(nèi)是“譜精度”理論上誤差小到機(jī)器精度可以只用每波長(zhǎng) 3 到 5 個(gè)網(wǎng)格就能算出波形不錯(cuò)的傳播結(jié)果有限差分通常需要 10 到 15 個(gè)網(wǎng)格。時(shí)間推進(jìn)用 k-space 修正項(xiàng)解決了偽譜法顯式時(shí)間推進(jìn)的穩(wěn)定性限制允許的 CFLCourant-Friedrichs-Lewy數(shù)比普通偽譜法大不少意味著同樣模擬時(shí)長(zhǎng)可以少走很多時(shí)間步。這套方法的代價(jià)有兩個(gè)第一傅里葉變換隱含周期性邊界條件所以必須用 PML完美匹配層吸收邊界網(wǎng)格四周要墊一圈吸收層第二媒質(zhì)參數(shù)聲速、密度在網(wǎng)格間突變時(shí)會(huì)出現(xiàn)吉布斯振蕩處理流體-固體邊界時(shí)要小心。理解了這兩點(diǎn)后面很多參數(shù)設(shè)置就不用死記硬背了。2.2 用 makeGrid 構(gòu)造計(jì)算域dx、Nx、Ny 與 CFL 數(shù)的關(guān)系k-Wave 的核心數(shù)據(jù)結(jié)構(gòu)圍繞kWaveGrid展開。最常見的錯(cuò)誤是一上來就抄示例代碼把kgrid的參數(shù)改一改就跑完全不知道自己設(shè)置的網(wǎng)格在物理上意味著什么。先看最小配置% 定義計(jì)算域尺寸和網(wǎng)格步長(zhǎng) Nx 256; % x 方向網(wǎng)格數(shù)決定空間分辨率 Ny 256; % y 方向網(wǎng)格數(shù) dx 0.1e-3; % 網(wǎng)格步長(zhǎng) 0.1 mm決定計(jì)算域邊長(zhǎng) 25.6 mm dy 0.1e-3; kgrid kWaveGrid(Nx, dx, Ny, dy); % 時(shí)間步長(zhǎng)與 CFL 數(shù)0.3 是 k-Wave 官方推薦上限 c 1500; % 媒質(zhì)聲速m/s [kgrid.t_array, dt] makeTime(kgrid, c, 0.3);這里makeTime返回兩個(gè)東西kgrid.t_array是整個(gè)時(shí)間序列向量dt是時(shí)間步長(zhǎng)。CFL 數(shù)的物理意義是c * dt / dx也就是一個(gè)時(shí)間步內(nèi)聲波跨過多少個(gè)網(wǎng)格。0.3 意味著每個(gè)時(shí)間步波走 0.3 個(gè)網(wǎng)格這個(gè)值是穩(wěn)定性和精度的折中大于 0.5 會(huì)不穩(wěn)定小于 0.1 則時(shí)間步太多白白增加計(jì)算量。計(jì)算域總邊長(zhǎng)就是Nx * dx。光聲成像里常見的情況是樣品 10 mm 左右要分辨 50 微米的特征那 Nx 就要 200 以上。網(wǎng)格數(shù)每翻一倍內(nèi)存漲 4 倍三維是 8 倍時(shí)間步數(shù)也翻倍所以先定dx再定Nx不要反過來。2.3 medium、source、sensor 三大對(duì)象的最小配置k-Wave 仿真需要三個(gè)結(jié)構(gòu)體媒質(zhì)參數(shù)medium、激勵(lì)源source、記錄點(diǎn)sensor。光聲仿真和超聲仿真的一個(gè)本質(zhì)區(qū)別是光聲的源是初始?jí)毫Ψ植約ource.p0一個(gè)時(shí)間點(diǎn)就賦完值后面不再注入能量超聲仿真則是source.p壓力源或source.u速度源在邊界上持續(xù)激勵(lì)。% 媒質(zhì)默認(rèn)是水均勻聲速 1500 m/s medium.sound_speed 1500; % 標(biāo)量 均勻媒質(zhì)矩陣 非均勻 medium.density 1000; % 密度影響聲阻抗匹配 % 源光聲用初始?jí)毫Ψ植?p0單位 Pa source.p0 p0_map; % Nx x Ny 的二維矩陣 % 傳感器經(jīng)典圓形陣列掃描 sensor.mask zeros(Nx, Ny); sensor.mask(100, 50:210) 1; % 在 x100 這一列的 50~210 行放置點(diǎn)探頭 sensor.record {p, p_max}; % 記錄時(shí)域壓力 p 和峰值 p_maxsensor.mask為 1 的位置就是探頭位置有幾個(gè) 1 就有幾個(gè)傳感器通道。這里有個(gè)新手容易忽略的點(diǎn)sensor.record里寫p會(huì)記錄所有時(shí)刻、所有通道的完整時(shí)間序列數(shù)據(jù)量是通道數(shù)×?xí)r間步數(shù)再乘以 8 字節(jié)double。如果只需要最終圖像只記錄p_max或p_final能省下幾 GB 內(nèi)存。3. 裝好 K-Wave 1.2.1 并跑通第一個(gè)二維光聲算例3.1 下載解壓后第一次啟動(dòng)addpath 與保存路徑的坑K-Wave 不提供安裝程序下載壓縮包解壓后把整個(gè)目錄加進(jìn) MATLAB 路徑就算裝完。1.2.1 是 2020 年前后的穩(wěn)定版本文件組織比早期版本清晰很多頂層有k-Wave、matlab、examples三個(gè)目錄。需要注意工具箱函數(shù)在k-Wave/matlab子目錄里只 addpath 頂層是不夠的。% 把 k-Wave 所有子目錄一次性加入路徑 addpath(genpath(D:\toolbox\k-Wave-toolbox-1.2.1)); savepath; % 保存到 MATLAB 默認(rèn)路徑避免每次重啟重新 addpathsavepath這一步很多人會(huì)跳過結(jié)果下次啟動(dòng) MATLAB 后發(fā)現(xiàn)kspaceFirstOrder2D又變成“未定義函數(shù)”。另外如果你重裝過 MATLAB 或換了電腦原路徑失效savepath會(huì)報(bào)錯(cuò)——這時(shí)重新addpath(genpath(...))再用savepath覆蓋即可。驗(yàn)證是否裝對(duì)的命令是which kspaceFirstOrder2D返回帶完整路徑的 .m 文件說明一切正常。3.2 用 kspaceFirstOrder2D 跑通一個(gè)最小光聲算例下面這個(gè)算例是 k-Wave 官方示例example_pr_2D_tr_circular_array.m的精簡(jiǎn)版去掉了一切不必要的東西物理上等價(jià)于一個(gè)半徑 2 mm 的均勻圓形吸收體在 0 時(shí)刻瞬間熱膨脹產(chǎn)生初始?jí)毫χ車劫|(zhì)中 64 個(gè)探頭繞圈接收信號(hào)。clear; clc; % 網(wǎng)格與媒質(zhì) Nx 128; dx 0.2e-3; Ny 128; dy dx; kgrid kWaveGrid(Nx, dx, Ny, dy); medium.sound_speed 1500; medium.density 1000; % 初始?jí)毫Π霃?10 個(gè)網(wǎng)格的圓盤壓力幅值 1 Pa p0_map zeros(Nx, Ny); [xx, yy] meshgrid(1:Nx, 1:Ny); disc ( (xx - Nx/2).^2 (yy - Ny/2).^2 10^2 ); p0_map(disc) 1; source.p0 p0_map; % 傳感器掩膜半徑 50 個(gè)網(wǎng)格的圓弧上放 16 個(gè)探頭 sensor_mask zeros(Nx, Ny); theta linspace(0, 2*pi, 17); theta theta(1:end-1); sx round(Nx/2 50 * cos(theta)); sy round(Ny/2 50 * sin(theta)); for k 1:16 sensor_mask(sx(k), sy(k)) 1; end sensor.mask sensor_mask; sensor.record {p}; % 時(shí)間序列 [kgrid.t_array, dt] makeTime(kgrid, medium.sound_speed, 0.3); % 跑仿真 sensor_data kspaceFirstOrder2D(kgrid, medium, source, sensor); % 畫一個(gè)探頭收到的信號(hào) plot(sensor_data(1, :) * 1e3); % 轉(zhuǎn)成 kPa 便于觀察 xlabel(時(shí)間步); ylabel(壓力 (kPa));這個(gè)代碼要注意的點(diǎn)meshgrid生成的xx是按列變化的畫圓盤條件(xx - Nx/2).^2 (yy - Ny/2).^2里xx是 x 坐標(biāo)、yy是 y 坐標(biāo)方向不要搞反否則圓盤會(huì)變成橢圓。source.p0的單位是 Pa默認(rèn)媒質(zhì)密度 1000 kg/m3 時(shí)輸出壓力也是 Pa量綱一致性由工具箱內(nèi)部保證。kspaceFirstOrder2D是核心求解函數(shù)函數(shù)名后綴 2D 表示二維求解。它內(nèi)部會(huì)自動(dòng)檢測(cè)source和sensor里定義了哪些字段沒有定義的字段用默認(rèn)值。跑完后sensor_data的維度是 16×?xí)r間步數(shù)第 i 行第 j 列是第 i 個(gè)探頭在第 j 個(gè)時(shí)間步收到的壓力。3.3 版本差異1.2.1 和后續(xù)版本的接口改動(dòng)1.2.1 之后的 1.3、1.4 版本主要加了 GPU 加速、彈性波求解、以及一些性能優(yōu)化。接口層面1.2.1 的參數(shù)名和 1.4 基本兼容主要區(qū)別在kWaveSimulationOptions這個(gè)可選參數(shù)的寫法。1.3 之后把PMLSize這類參數(shù)統(tǒng)一收了進(jìn)去1.2.1 里可以直接寫在kspaceFirstOrder2D的第五個(gè)輸入?yún)?shù)位置。老代碼升級(jí)時(shí)最常見的報(bào)錯(cuò)是DataCast參數(shù)位置變了——1.2.1 支持DataCast, single以半精度運(yùn)行新版把數(shù)據(jù)類型強(qiáng)制和UseGPU綁定了。如果你手上是 1.2.1建議保持雙精度1.2.1 的 single 模式在某些邊緣場(chǎng)景有數(shù)值累積誤差不值得為了那點(diǎn)內(nèi)存省事。4. 參數(shù)調(diào)優(yōu)CFL、PML 與 record_mask 的三組核心取舍4.1 CFL、PML 和網(wǎng)格步長(zhǎng)先定參數(shù)再寫代碼CFL 數(shù)的選取直接決定計(jì)算量。makeTime(kgrid, c, cfl)里傳 0.3 是通用推薦值但你完全可以按場(chǎng)景調(diào)如果只關(guān)心信號(hào)到達(dá)時(shí)間和大概波形0.5 夠用如果要做精確幅度對(duì)比、驗(yàn)證重建算法用 0.2 更穩(wěn)。代價(jià)是時(shí)間步數(shù)反比于 CFL0.2 比 0.4 多一倍時(shí)間步跑一次三維仿真可能多幾個(gè)小時(shí)。PML 默認(rèn)是 20 個(gè)網(wǎng)格。這個(gè)默認(rèn)值在大多數(shù)情況夠用但有兩個(gè)例外一是媒質(zhì)聲速差異特別大比如軟組織和骨的界面反射波比較強(qiáng)PML 要加到 30 甚至 40二是傳感器離邊界很近時(shí)早期到達(dá)探頭的波會(huì)混入 PML 反射這時(shí)與其加厚 PML不如把計(jì)算域整體擴(kuò)大一圈。% 顯式指定 PML 厚度為 30 個(gè)網(wǎng)格 kspaceFirstOrder2D(kgrid, medium, source, sensor, ... PMLSize, 30, PlotScale, [-1e-6, 1e-6]);PlotScale參數(shù)只影響可視化不影響計(jì)算。調(diào)試階段建議開著 Plot但正式批量跑仿真時(shí)一定要關(guān)掉——繪圖開銷能占掉總運(yùn)行時(shí)間的 10% 到 20%。4.2 記錄哪些物理量record.mask 與 record.record 的取舍sensor.record里能寫p、p_max、p_mean、p_final、u等多種量。用得最多的是p和p_max。區(qū)別是p給每個(gè)通道存完整 A-line 信號(hào)做重建和定量分析必須用它p_max只存每個(gè)通道的最大值適合快速看個(gè)大致的圖像輪廓內(nèi)存開銷小一到兩個(gè)數(shù)量級(jí)。record.mask是另一個(gè)容易被忽略的參數(shù)。默認(rèn)情況下record.mask等于sensor.mask也就是只記錄傳感器位置。但對(duì)于逆問題研究你往往需要知道某個(gè)特定區(qū)域內(nèi)部任意時(shí)刻的壓力場(chǎng)這時(shí)record.mask設(shè)為感興趣區(qū)域即可不必讓整個(gè)計(jì)算域都存下來。% 只記錄計(jì)算域中心 32x32 區(qū)域的完整時(shí)域壓力 record_mask zeros(Nx, Ny); record_mask(Nx/2-15:Nx/216, Ny/2-15:Ny/216) 1; sensor.record {p, p_max}; kspaceFirstOrder2D(kgrid, medium, source, sensor, ... RecordMask, record_mask);這個(gè)功能在驗(yàn)證“聲速不均勻?qū)χ亟ǖ挠绊憽边@類課題時(shí)非常實(shí)用你可以在一個(gè)位置放真正的傳感器用來重建同時(shí)把整個(gè)場(chǎng)內(nèi)部記錄下來做誤差分析一次仿真拿到兩組數(shù)據(jù)。4.3 從二維起步還是直接上三維內(nèi)存估算公式三維仿真的內(nèi)存開銷由網(wǎng)格數(shù)主導(dǎo)。每保存一個(gè)完整時(shí)間序列內(nèi)存占用大約是 網(wǎng)格數(shù) × 時(shí)間步數(shù) × 8 字節(jié)。通常計(jì)算域 256×256×256、時(shí)間步 500 步、只存p_max時(shí)內(nèi)存占用約 2563 × 500 × 8 67 GB這還只是傳感器數(shù)據(jù)的一部分。真正跑起來求解器內(nèi)部的中間變量還要額外占用幾個(gè)網(wǎng)格規(guī)模的數(shù)組所以 256 立方、500 步的仿真在 32 GB 內(nèi)存機(jī)器上基本跑不動(dòng)。我的經(jīng)驗(yàn)是先在二維把算法流程調(diào)通確認(rèn)物理模型正確后再轉(zhuǎn)三維。三維仿真前先用下面的公式估算峰值內(nèi)存峰值內(nèi)存大約是Nx*Ny*Nz*(t_steps10)*8字節(jié)的 2 到 3 倍PML、傅里葉變換臨時(shí)數(shù)組都會(huì)吃內(nèi)存% 估算三維仿真的峰值內(nèi)存GB Nx 128; Ny 128; Nz 128; t_steps 300; peak_mem_gb Nx * Ny * Nz * (t_steps 10) * 8 * 2.5 / 1e9; fprintf(預(yù)估峰值內(nèi)存: %.1f GB\n, peak_mem_gb);如果估算結(jié)果超過物理內(nèi)存的 70%兩條路一是降分辨率把dx從 0.1 mm 放寬到 0.15 mm網(wǎng)格數(shù)直接砍掉大半二是用kspaceFirstOrder3D的DataCast, single參數(shù)半精度運(yùn)行前提是 MATLAB 版本支持內(nèi)存直接減半。注意 single 精度下p0幅值的相對(duì)誤差大約在 1e-6 量級(jí)對(duì)于絕大多數(shù)光聲仿真完全夠用。5. 避坑實(shí)錄K-Wave 1.2.1 最常見的五個(gè)翻車現(xiàn)場(chǎng)5.1 現(xiàn)象運(yùn)行時(shí)報(bào) Undefined function or variable kspaceFirstOrder2D原因工具箱沒加進(jìn)路徑或加了頂層目錄但沒加子目錄。1.2.1 的函數(shù)在k-Wave/matlab子目錄里只addpath(D:\...\k-Wave-toolbox-1.2.1)不行必須用genpath把所有子目錄遞歸加入。解決addpath(genpath(...k-Wave-toolbox-1.2.1)); savepath;然后which kspaceFirstOrder2D確認(rèn)返回路徑。注意savepath會(huì)覆蓋 MATLAB 的pathdef.m文件如果以前配置過其他工具箱路徑先備份再操作。5.2 現(xiàn)象仿真跑了一多半MATLAB 直接報(bào) Out of Memory 崩潰原因最常見的是三維仿真沒做內(nèi)存估算或者二維仿真把整個(gè)時(shí)間序列在內(nèi)存里堆積了一整份。sensor.record {p}在通道數(shù) 256、時(shí)間步 1000 時(shí)就需要 2 GB 存一次數(shù)據(jù)而內(nèi)部求解器還有多組臨時(shí)數(shù)組是它的數(shù)倍。解決先用上一章的內(nèi)存公式估算再考慮只記錄感興趣時(shí)段的信號(hào)——用sensor.time_start和sensor.time_end截取時(shí)間窗口比如光聲信號(hào)在 100 微秒內(nèi)到達(dá)探頭就不用記錄從 0 到 500 微秒的全部信號(hào)。sensor.time_start 20表示跳過前 20 個(gè)時(shí)間步再開始記錄內(nèi)存直接按比例下降。5.3 現(xiàn)象仿真順利跑完但所有傳感器信號(hào)全是零原因source.p0初始?jí)毫Ψ植荚O(shè)置區(qū)域與sensor.mask位置重疊或者p0全為零。更隱蔽的原因是source.p0賦值的矩陣維度與網(wǎng)格不一致——source.p0必須是一個(gè)Nx × Ny三維是Nx × Ny × Nz的矩陣而不是一個(gè) 128×128 的物理坐標(biāo)數(shù)組。解決調(diào)試時(shí)先畫一下imagesc(source.p0)確認(rèn)初始?jí)毫Ψ植颊_再看sensor.mask在不在p0附近的傳播路徑上最后檢查sensor_data里是否有數(shù)值量級(jí)異常正常光聲信號(hào)峰值在 Pa 到 kPa 量級(jí)。5.4 現(xiàn)象信號(hào)幅度看起來不對(duì)勁重建圖像一片模糊原因source.p0的單位和sensor.record {p}的輸出單位之間沒有做量綱換算或者medium.density、medium.sound_speed設(shè)成了非物理值。k-Wave 內(nèi)部按一致單位制計(jì)算默認(rèn)長(zhǎng)度是米、時(shí)間是秒、壓力是 Pa、密度是 kg/m3。如果dx 0.1e-3毫米網(wǎng)格但medium.sound_speed 1500 * 1e-3那時(shí)間步長(zhǎng)會(huì)出問題聲波每時(shí)間步走的距離就完全對(duì)不上。解決所有物理量嚴(yán)格換到 SI 基本單位。密度 1000 kg/m3、聲速 1500 m/s、壓力 1 Pa這樣輸出壓力直接就是 Pa。5.5 現(xiàn)象同一段代碼換電腦后結(jié)果和之前不一樣原因k-Wave 默認(rèn)會(huì)嘗試用并行池parpool跑kspaceFirstOrder2D并行分塊方式不同浮點(diǎn)累加順序就不同結(jié)果在小數(shù)點(diǎn)后第 10 位開始分叉。這個(gè)差異本身不影響物理結(jié)論但如果一個(gè)人在做定量對(duì)比時(shí)前后用了兩套硬件會(huì)誤判成算法改動(dòng)產(chǎn)生的差異。解決統(tǒng)一用NumThreads, 1強(qiáng)制單線程跑或者明確記錄每次仿真的 MATLAB 版本、CPU 型號(hào)、是否開啟并行。做重建算法驗(yàn)證時(shí)所有對(duì)比仿真必須在同一環(huán)境、同一線程數(shù)下跑完。6. 進(jìn)階驗(yàn)證把仿真結(jié)果和解析解對(duì)上你的工具箱才算真正裝好了光聲仿真跑通不難難的是確認(rèn)你真的沒有在某些邊界條件上犯錯(cuò)。最可靠的驗(yàn)證方式是用 k-Wave 自帶的一組解析對(duì)照算例均勻媒質(zhì)中一個(gè)球形或圓柱形吸收體在遠(yuǎn)場(chǎng)條件下傳感器收到的壓力信號(hào)可以由解析公式給出兩者對(duì)比誤差應(yīng)該在 1% 以內(nèi)。具體做法是先跑二維圓柱形吸收體初始?jí)毫鶆蚍植加?jì)算域四周 PML 墊 40 層傳感器放在距離吸收體中心 30 個(gè)網(wǎng)格的位置分別記錄時(shí)域信號(hào)p(t)。驗(yàn)證腳本如下% 解析驗(yàn)證均勻媒質(zhì)圓柱吸收體的光聲壓力信號(hào) kgrid kWaveGrid(256, 0.1e-3, 256, 0.1e-3); medium.sound_speed 1500; medium.density 1000; p0 zeros(256, 256); [x, y] meshgrid(1:256, 1:256); p0((x-128).^2 (y-128).^2 12^2) 1; source.p0 p0; sensor.mask zeros(256, 256); sensor.mask(128, 168) 1; % 放在圓柱外 40 個(gè)網(wǎng)格處 sensor.record {p}; [kgrid.t_array, dt] makeTime(kgrid, 1500, 0.2); sensor_data kspaceFirstOrder2D(kgrid, medium, source, sensor, ... PMLSize, 40, PlotSim, false); % 與解析解對(duì)比 c 1500; R 12 * 0.1e-3; d 40 * 0.1e-3; t kgrid.t_array * 1e6; % 轉(zhuǎn)微秒 p_analytic (d - c * (t*1e-6)) ./ (2 * d) .* (abs(t*1e-6 - d/c) R/c); plot(t, sensor_data(1,:), t, p_analytic); legend(k-Wave,解析解);如果兩條線在時(shí)域上重合說明你的安裝和物理設(shè)置都正確。這一步花的時(shí)間不多卻能幫你過濾掉大量隱性配置錯(cuò)誤——我見過一個(gè)同事用了半年 k-Wave結(jié)果某天發(fā)現(xiàn)自己的medium.density一直漏填導(dǎo)致所有仿真都默認(rèn)密度 1000對(duì)比實(shí)驗(yàn)里密度差異組全部白做了。別嫌這一步麻煩值得做。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取