kspaceSecondOrder 模拟函数

概述

kspaceSecondOrder 模拟线性压缩波在一维、二维或三维均质声学介质中的时域传播，给定四个输入结构：kgrid、介质、声源和传感器。计算基于具有幂律吸收的介质的精确二阶 k 空间模型。在每个时间步长（由 kgrid.dt 和 kgrid.Nt 或 kgrid.t_array 定义），记录并存储由 sensor.mask 定义的位置处的压力。如果 kgrid.t_array 设置为 “auto”（自动），则该数组将通过 kWaveGrid 类的 makeTime 方法自动生成。为防止波浪缠绕，可以通过将可选输入项 "ExpandGrid "设置为 true，自动将计算域扩大 2 倍。

可以通过为 source.p0 指定一个任意数值矩阵（与计算网格大小相同）来指定初始压力分布。同样，也可以使用 source.dp0dt 指定初始压力梯度。压力将以 sensor.mask 所定义的传感器位置的时间序列数组形式返回。该数组被指定为二进制矩阵（即与计算网格大小相同的 "1 "和 "0 "矩阵），代表计算网格内将收集数据的网格点。传感器数据使用 MATLAB 的标准列式线性矩阵索引排序，记录数据的索引为 sensor_data(sensor_position,time)。也可以通过将 sensor.record 设置为 {‘p’, ‘p_final’} 来获取整个计算网格的最终压力场。在这种情况下，输出的 sensor_data 将以结构体形式返回，其输出将作为结构字段 sensor_data.p 和 sensor_data.p_final 附加到结构体中。

与一阶模拟函数 kspaceFirstOrder1D、kspaceFirstOrder2D 和 kspaceFirstOrder3D 相比，kspaceSecondOrder 只适用于同质介质，功能较少。 不过，它对于均质吸收介质是精确的，计算效率更高，并允许指定初始压力梯度。

Inputs & Outputs

详见http://www.k-wave.org/documentation/kspaceSecondOrder.php

模拟函数比较

本示例简要比较了 kspaceFirstOrder2D 和 kspaceSecondOrder 模拟函数。它以同质传播介质和使用二进制传感器掩码示例为基础。

关于 k-Wave 中使用的二阶模型的更详细讨论，请参见 Treeby, B. E. 和 Cox, B. T., “A k-space Green’s function solution for acoustic initial value problems in homogeneous media with power law absorption,” J. Acoust. 美国声学学会》，第 129 卷，第 6 期，第 3652-3660 页，2011 年。

关于建模函数

在前面的示例中，声学模拟都是使用 kspaceFirstOrder2D 进行的。该函数基于使用三个耦合一阶偏微分方程（质量守恒、动量守恒和压力密度关系）对粒子速度、声学密度和声学压力进行的顺序计算。对于均质介质，这些方程也可以合并为一个二阶声波方程，并使用格林函数方法求解。kspaceSecondOrder 函数是专门针对初值问题的格林函数求解的高效数值实现。对于均质介质，两种方法的结果相同。

与 kspaceFirstOrder2D 相比，kspaceSecondOrder 的计算效率更高，时间步长也可以任意大（因为求解是精确的）。它还允许定义初始压力和初始压力梯度（参见设置初始压力梯度示例）。 幂律吸收的编码也更加精确（请参阅 “幂律吸收建模示例”）。不过，kspaceSecondOrder 的功能略少于其一阶对应功能（kspaceFirstOrder1D、kspaceFirstOrder2D 和 kspaceFirstOrder3D）。例如，它只支持同质介质和二进制传感器掩码，不支持时间变化源或全部可视化选项。

二阶代码没有实现吸收边界层（有关 kspaceFirstOrder2D 中使用的完全匹配层的更多详情，请参阅控制吸收边界层示例）。相反，为了延迟波包裹的出现，它扩大了计算网格的尺寸。虽然这会降低计算效率，但却能保持解法的精确性。网格扩展可通过设置可选输入参数 "ExpandGrid "为 true 来控制。请注意，如果模拟时间长于波在网格扩展后的传播时间，则仍会出现包裹波。

运行模拟

kspaceSecondOrder 函数（适用于所有维度）采用与一阶模拟函数相同的输入结构。例如，在运行均质传播介质示例时，首先使用 cart2grid 将笛卡尔传感器掩码转换为二进制传感器掩码，然后调用 kspaceSecondOrder，并将 "ExpandGrid "设置为 true。

% convert the cartesian sensor mask to a binary sensor mask
sensor.mask = cart2grid(kgrid, sensor.mask);

% run the simulation using the second order code
sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'ExpandGrid', true);

数值比较

当介质不吸水且未使用波浪缠绕补偿时，两个模拟函数得出的结果与机器精度相同。在第一个传感器位置记录的信号和两者之间的差异图如下（在示例 m 文件中设置 example_number = 1）。

当介质为吸收性介质且使用完全匹配层减轻波包裹时，使用 kspaceFirstOrder2D 计算出的数值结果存在微小误差。这一误差主要是由于幂律吸收项的计算造成的（见 “幂律吸收建模示例”）。下面给出了在第一个传感器位置记录的信号和它们之间的差值图（在示例 m 文件中设置 example_number = 4）。

设置初始压力梯度

kspaceSecondOrder 函数允许定义初始压力及其时间梯度的初始条件。初始条件定义为与计算网格大小相同的数值矩阵。初始压力分配给 source.p0，初始压力梯度分配给 source.dp0dt。

在本例中，makeDisc 用于创建由两个直径不同的小圆盘组成的初始源分布。首先分配给源结构的 p0 字段，然后分配给 dp0dt 字段，最后分配给这两个字段。可选输入 "MeshPlot（网格绘图）"设置为 “true”，以实现压力场的三维可视化。如果 "PlotFrames（绘图帧）"设置为 true，则会在新图中生成单个帧。下面是两个不同 t 值下的压力场图。

补充：如何导出传感器获得的数据

要将传感器记录的数据导出，可以按照以下步骤进行操作：

1. 确保已经加载了包含传感器数据的MATLAB变量。

2. 使用MATLAB的save函数将数据保存为MAT文件：

save('sensor_data.mat', 'sensor_data');

这将在当前工作目录中创建一个名为"sensor_data.mat"的MAT文件，并将传感器数据保存在其中。

3. 如果想将数据导出为其他格式（如CSV、TXT等），可以使用MATLAB的dlmwrite函数：

dlmwrite('sensor_data.csv', sensor_data, 'delimiter', ',');

这将在当前工作目录中创建一个名为"sensor_data.csv"的CSV文件，并将传感器数据保存在其中。可以根据需要调整分隔符和其他参数。

4. 如果想将数据导出为其他格式（如Excel文件），可以使用MATLAB的writematrix函数（需要MATLAB R2019a或更高版本）：

writematrix(sensor_data, 'sensor_data.xlsx');

这将在当前工作目录中创建一个名为"sensor_data.xlsx"的Excel文件，并将传感器数据保存在其中。

根据需要，可以选择适合的导出格式，并在导出时调整文件名、分隔符等参数。

简单作图

要绘制第一个传感器记录的时变压力图，可以在MATLAB命令窗口中，加载示例数据集（或者自己的数据集）并设置相关参数，例如：

load('example_data.mat'); % 加载示例数据集
sensor_index = 1; % 第一个传感器的索引

1. 提取第一个传感器的压力数据，并获取时间步长：

sensor_data = sensor_data(:, sensor_index); % 提取第一个传感器的压力数据
dt = kgrid.dt; % 获取时间步长

2. 绘制第一个传感器的时变压力图：

figure;
t = (0:length(sensor_data)-1) * dt; % 时间向量
plot(t, sensor_data);
xlabel('时间（s）');
ylabel('压力');
title('第一个传感器的时变压力图');

这样，就可以绘制出第一个传感器记录的时变压力图。根据需要，可以调整图形的标题、横轴和纵轴标签等。