k-wave三维模拟及AS模拟

三维模拟

三维模拟与一维和二维模拟类似。介质离散化同样由 kWaveGrid 执行，Z 维有两个额外的输入。同样，介质属性矩阵也多了一个维度。矩阵索引约定（x, y, z）是二维矩阵索引约定（x, y）的延伸，矩阵行和列数据分别对应 x 轴和 y 轴。

同样，传感器掩膜可以是直角坐标列表、二进制掩膜或立方体两个对角的网格坐标。k-Wave 工具箱中包含多个用于创建三维图形的函数。例如，makeSphere 返回使用中点圆算法创建的单像素球壳的二进制地图，makeCartSphere 返回球面上任意数量点的均匀分布阵列的笛卡尔位置。在本例中，明确创建了一个线性传感器阵列。

通过调用 kspaceFirstOrder3D 并输入上文定义的输入值，即可进行计算。默认情况下，会显示水平面、中线面和额面的波场可视化效果以及状态栏。

函数运行时，状态更新和计算参数会打印到命令行。通过将可选输入 "DataCast "设置为 “single”，可以减少计算时间（更多详情请参见优化 k-Wave 性能示例）。

时间循环结束后，函数将返回在 sensor.mask 所定义的每个传感器位置上记录的时间序列。排序同样取决于使用的是直角坐标传感器掩码还是二进制传感器掩码。可视化效果如下。

如果使用的是笛卡尔传感器掩码，则在每个时间步长内，传感器点上的声场值都是通过插值获得的。插值模式可通过将可选输入参数 "CartInterp "设置为 “线性”（默认）或 “最近”（本例中使用）来更改。在模拟过程中，使用笛卡尔传感器掩模和二进制传感器掩模之间的性能差异很小。但是，如果使用线性插值的笛卡尔传感器掩膜，三角测量点的计算会大大延长预计算时 间，尤其是在三维模拟中（请参阅下面将 "CartInterp "设置为 "线性 "后运行的同一示例的输出）。可以使用函数 cart2grid 和 grid2cart 将笛卡尔传感器掩码转换为二进制传感器掩码（反之亦然）。

在二维和三维模拟中，如果可选输入 "CartInterp "设置为 “最近”，则笛卡尔点将在内部映射到二进制网格上。如果多个笛卡尔传感器点映射到同一位置，重复的点将被丢弃。因此，如果笛卡尔传感器掩码具有密集的点阵列，返回的时间序列数可能与原始的笛卡尔传感器点数不一致。被丢弃的重复点的详细信息会打印到命令行中。这种情况可以通过提高计算网格的分辨率或减少笛卡尔传感器点的堆积来避免。

轴对称坐标系模拟

创建 k 空间网格并定义介质属性

轴对称模拟的执行方式与二维模拟类似。不过，在轴对称坐标系中，x 维对应轴向，y 维对应径向，如下图所示。坐标系关于 x 轴是旋转对称的，因此 y 轴上的一个点对应三维空间中的一个连续圆。

定义好网格参数后，kWaveGrid 将以与二维模拟相同的方式进行介质离散化。相比之下，对于 kspaceFirstOrder2D，笛卡尔点 y = 0 位于计算网格的中间。对于异质声波传播介质，介质属性以与计算网格大小相同的二维矩阵形式给出。在本例中，介质属性被划分为两个半空间。

定义初始压力分布和传感器掩膜

与二维矩阵一样，初始压力分布也是通过矩阵设置的，该矩阵包含计算域内每个网格点的初始压力值。不过，对于轴对称代码来说，初始压力是围绕矩阵的第一列旋转对称的。例如，如果将初始压力设置为第一列中的一个网格点（即在 x 轴上，y = 0），这相当于在三维模拟中设置一个点源。反之，如果将初始压力设置为任何其他列中的单个网格点，则相当于在三维空间中设置了一个环。在本例中，初始压力被设置为半个圆盘的形状，相当于三维中的一个球。

同样，传感器掩码定义了在每个时间步长记录压力场的位置，可以是笛卡尔坐标列表、二进制掩码或矩形两个对角的网格坐标。在本例中，定义了一个圆形的笛卡尔传感器掩码。初始压力和传感器掩膜的可视化效果如下所示。请注意，对于轴对称代码，在径向（y）方向上，完全匹配层（PML）仅应用于域的外边缘。默认 PML 大小为 20 个网格点。

运行模拟

向 kspaceFirstOrderAS 传递四个输入结构（kgrid、介质、源和传感器）后，计算就开始了。默认情况下，将显示传播波场的可视化效果和状态栏。函数运行时，状态更新和计算参数会打印到命令行。

代码

%对称 Axisymmetric Coordinate
% create the computational grid
Nx = 128;           % number of grid points in the axial (x) direction
Ny = 64;            % number of grid points in the radial (y) direction
dx = 0.1e-3;        % grid point spacing in the axial (x) direction [m]
dy = 0.1e-3;        % grid point spacing in the radial (y) direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500 * ones(Nx, Ny);       % [m/s]
medium.sound_speed(Nx/2:end, :) = 1800;         % [m/s]
medium.density = 1000 * ones(Nx, Ny);           % [kg/m^3]
medium.density(Nx/2:end, :) = 1200;             % [kg/m^3]

% create initial pressure distribution in the shape of a disc - this is
% generated on a 2D grid that is doubled in size in the radial (y)
% direction, and then trimmed so that only half the disc is retained
source.p0 = 10 * makeDisc(Nx, 2 * Ny, Nx/4 + 8, Ny + 1, 5);
source.p0 = source.p0(:, Ny + 1:end);

% define a Cartesian sensor mask with points in the shape of a circle
sensor.mask = makeCartCircle(40 * dx, 50);

% remove points from sensor mask where y < 0
sensor.mask(:, sensor.mask(2, :) < 0) = [];

% run the simulation
input_args = {'PlotLayout', true,'RecordMovie', true, 'MovieName', 'AS'};
sensor_data = kspaceFirstOrderAS(kgrid, medium, source, sensor, input_args{:});

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;