k-Wave基础

新手，学习过程记录，主要参考k-Wave user manual，恐有错漏。

---

前言

k-Wave是一个开源的第三方Matlab工具箱，可以在一维、二维、三维的空间中对声波传播进行时域仿真。该工具箱功能全面，其核心是一个可以解释线性和非线性的波传播、非同质材料参数的任意分布和幂律声吸收的数学模型。

模拟函数的调用需要四个输入结构：网格点（kgrid）、介质（medium）、声源（source）和传感器（sensor）。这四个结构分别定义了计算网格点的属性、传播介质的材质属性、任意声源的位置和属性以及传感器点的位置和属性。

一、输入结构

1. kgrid

kgrid 定义了计算格点的属性。格点是为了便于模拟，将连续介质空间离散化的产物，kgrid 将连续介质均匀分为一个个格点。

创建kgrid需要使用函数makeGrid，以二维kgrid为例：

% create computational grid for a 2D simulation
kgrid = makeGrid(Nx, dx, Ny, dy);

其中Nx、Ny是对应维度上网格点的数量，dx、dy是网格点的空间尺寸，则 kgrid 在x方向上的实际尺寸应为kgrid.x_size = kgrid.Nx * kgrid.dx。

makeGrid函数会返回一个kWaveGrid类，该类包含的属性如下：

需要注意的是，生成的kgrid中默认包含一个PML（perfectly matched layer）层，该层位于网格点的内部边界，声波传播到PML层时会被全部吸收。二维的kgrid中PML层宽度为20个网格点，三维为10个网格点，模拟时注意不要将声源或是传感器放置在PML层内，否则会发生难以解释的现象。也可以在模拟时通过将 ‘PMLInside’ 选项设置为 false 来使PML层位于网格点外部，这样就可以不用考虑边界问题了。

在 kWaveGrid 对象被创建后，唯一可以调整的参数只剩下kgrid.t_array，它描述的是模拟的时间值数组，在不指定的情况下会自动设置。默认情况下，模拟的总时长为声波以最小声速穿过最长网格对角线所需的时间，而时间步长基于CFL数（默认为0.3）和介质的最大声速，公式为kgrid.dt = CFL*dx_min/c0_max。

利用 makeTime 函数可以对模拟的时间参数进行设置。

% create the time array using makeTime setting the CFL and end time
kgrid.t_array = makeTime(kgrid, medium.sound_speed, CFL, t_end);

% create the time array using makeTime setting the end time
kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);

2. medium

medium定义了每一个网格点处的介质属性。可以被定义的参数共有五个，分别是等熵声速（sound_speed）、环境质量密度（density）、非线性参数（BonA）、幂律吸收系数或前因数（alpha_coeff）和幂律吸收指数（alpha_power），具体见下表：

除了 alpha_power 参数必须是一个标量外，其余属性既可以是标量（介质为均匀同质时）也可以是矩阵（介质为非同质或不均匀时）。如果介质为非同质，则输入矩阵的大小必须和网格点的大小相同。

3. source

source 结构定义了介质中的声源，k-Wave中共有三种可以使用的声源。

第一种是初始声压分布，即在仿真开始前在网格中某个区域设置一个初始声压脉冲，使用预设函数可以很方便地画一些几何图形，比如二维网格中可以用 makeDisc 画实心圆，makeCircle 画圆圈，三维网格中则为 makeBall 和 makeSphere。以下代码在三维网格中定义了一个球形的初始声压。

% define an initial pressure distribution using makeBall
source.p0 = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, radius);

第二种是声压随时间变化的声源。该类型声源需要两个输入参数，一是声源的网格点位置，二是声源声压随时间变化的函数。声源的网格点位置使用二值的矩阵来表示（1表示声源点，0表示非声源点），该矩阵应与kgrid具有相同的尺寸，保存在 source.p_mask 中。输入信号可以为一维数组（长度为kgrid.Nt，记录每个时间点的声压强度，并被所有的声源点使用），也可以为二维矩阵（分别设置每一个声源的声压，每一列为一个声源的声压）。

k-Wave中声源点按从上到下、从左到右的顺序标号，对上图所设的source.p_mask，其标号顺序如下。

第三种声源按我个人的理解是在第二种声源的基础上，可以定义不同方向的声压强度。利用source.u_mask定义声源的网格点位置，然后用source.ux，source.uy 和 source.uz 定义不同方向的声压强度。

% define the source mask to be a line across the top of the grid
source.u_mask = zeros(Nx, Ny);
source.u_mask(1, :) = 1;

source 结构的详细属性如下：

4. sensor

sensor 定义了传感器相关的属性，仿真时，传感器会记录该处在每一个时间点的声压强度。传感器在网格点中的位置使用 sensor.mask 来定义（需要与 kgrid 尺寸相同），以下是一些示例。

% define a 2D binary sensor mask in the shape of a line
x_offset = 25; % [grid points]
width = 50; % [grid points]
sensor.mask = zeros(Nx, Ny);
sensor.mask(x_offset, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;

% define a 2D binary sensor mask in the shape of an arc using makeCircle
x_pos = Nx/2; % [grid points]
y_pos = Ny/2; % [grid points]
radius = 20; % [grid points]
arc_angle = pi/2; % [radians]
sensor.mask = makeCircle(Nx, Ny, x_pos, y_pos, radius, arc_angle);

也可以用坐标的方法来定义，默认与 kgrid 尺寸相同。

% define a 2D Cartesian sensor mask
x = -10:2:10; % [m]
y = -10:2:10; % [m]
sensor.mask = [x; y];

二、仿真

设置完输入参数之后，使用 kspaceFirstOrder2D 和 kspaceFirstOrder3D 函数分别对二维网格和三维网格进行仿真。

% run the simulation in 2D
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% run the simulation in 3D
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor);

仿真返回值是一个二维矩阵，大小为仿真时间点数*传感器个数，记录了每个传感器在每个时间点接收到的声压强度。

更多需求可以查看k-Wave官网的案例：
http://www.k-wave.org/documentation/k-wave.php