k-wave学习：声波传播的时域仿真（以kspaceFirstOrder2D为例）

概述

kspaceFirstOrder2D 模拟压缩波在二维同质或异质声学介质中的时域传播，给定四个输入结构：kgrid、medium、source 和sensor。计算基于一阶 k 空间模型，该模型考虑了幂律吸收以及异质声速和密度。如果指定了 medium.BonA，还将模拟累积非线性效应。在每个时间步长（由 kgrid.dt 和 kgrid.Nt 或 kgrid.t_array 定义），记录并存储由 sensor.mask 定义的位置处的声场参数。如果 kgrid.t_array 设置为 "auto"（自动），则该数组将使用 kWaveGrid 类的 makeTime 方法自动生成。为防止离开域的一侧的波从另一侧重新引入（使用 FFT 计算波方程空间导数的结果），实施了一个称为完全匹配层（PML）的各向异性吸收边界层。这样就可以使用较小的计算网格进行无限域模拟计算。

对于均质介质，该公式是精确的，时间步长仅受完全匹配层的有效性限制。对于异质介质，解决方案是一种跃迁式伪谱方法，其 k 空间校正提高了计算时间导数的精度。与传统的伪谱时域方法相比，这种方法允许在相同精度水平上采用更大的时间步长。计算网格在空间和时间上都是交错的。

通过为 source.p0 分配一个任意数值矩阵（与计算网格大小相同），可以指定初始压力分布。时变压力源同样可以通过为 source.p_mask 指定一个二进制矩阵（即由 1 和 0 组成的矩阵，其大小与计算网格相同）来指定，其中 1 代表构成压力源一部分的网格点。然后将时变输入信号分配到 source.p。这可以是单个时间序列（在这种情况下应用于所有源元素），也可以是使用 MATLAB 标准列式线性矩阵索引排序的源元素之后的时间序列矩阵。时变速度源可以类似方式指定，其中源位置由 source.u_mask 指定，时变输入速度分配给 source.ux 和 source.uy。

字段值以 sensor.mask 所定义的传感器位置的时间序列数组形式返回。这可以通过三种不同的方式来定义 (1) 二进制矩阵（即由 1 和 0 组成的矩阵，其尺寸与计算网格相同），代表计算网格内将收集数据的网格点。(2) 以 [x1; y1; x2; y2] 的形式表示矩形两个对角的网格坐标。这相当于使用二进制传感器掩码覆盖同一区域，但输出的索引方式不同，如下所述。(3) 网格内的一系列笛卡尔坐标，指定了每个时间步长存储压力值的位置。如果笛卡尔坐标与网格点坐标不完全一致，输出值将通过插值计算得出。笛卡尔坐标点必须以 2 乘 N 矩阵的形式给出，分别对应 x 和 y 位置，其中笛卡尔原点假定位于网格中心。如果不需要输出，传感器输入可以用空数组[]代替。

如果 sensor.mask 以笛卡尔坐标集的形式给出，计算出的 sensor_data 将以相同的顺序返回。如果以二进制矩阵形式给出 sensor.mask，则会使用 MATLAB 的标准列向线性矩阵索引顺序返回 sensor_data。在这两种情况下，记录数据的索引都是 sensor_data(sensor_point_index,time_index)。对于二进制传感器掩码，可以使用 unmaskSensorData 将特定时间的场值还原到计算网格中的传感器位置。如果以矩形对角列表的形式给出 sensor.mask，则记录数据的索引为 sensor_data(rect_index).p(x_index, y_index, time_index)，其中 x_index 和 y_index 对应于矩形内的网格索引，如果指定了多个矩形，则 rect_index 对应于矩形个数。

默认情况下，记录的声压场会直接传递到输出 sensor_data。不过，也可以通过将 sensor.record 设置为形式为 {'p', 'u', 'p_max', ...}的单元数组来记录其他声学参数。例如，通过设置 sensor.record = {'p', 'u'} 可以同时返回粒子速度和声压。如果给定了 sensor.record，输出的 sensor_data 将以结构体的形式返回，不同的输出将作为结构体字段被附加。例如，如果 sensor.record = {'p'、'p_final'、'p_max'、'u'}，输出将包含 sensor_data.p、sensor_data.p_final、sensor_data.p_max、sensor_data.ux 和 sensor_data.uy 字段。大多数输出参数都记录在给定的传感器位置上，如果使用的传感器掩码定义为相对的矩形角，则索引为 sensor_data.field(sensor_point_index,time_index)或 sensor_data(rect_index).field(x_index,y_index,time_index)。平均量（'p_max'、'p_rms'、'u_max'、'p_rms'、'I_avg'）、'all'量（'p_max_all'、'p_min_all'、'u_max_all'、'u_min_all'）和最终量（'p_final'、'u_final'）除外。平均量的索引为 sensor_data.p_max(sensor_point_index)，或 sensor_data(rect_index).p_max(x_index, y_index)（如果使用矩形角），而最终量和 "所有 "量则返回整个网格，并且始终以 sensor_data.p_final(nx, ny) 作为索引，与传感器屏蔽类型无关。

kspaceFirstOrder2D 也可用于时间反演图像重建，方法是将任意传感器表面上记录的时变压力分配给输入字段 sensor.time_reversal_boundary_data。然后，在 sensor.mask 给定的传感器表面上，该数据将作为时变 Dirichlet 边界条件按时间反转顺序执行。边界数据的索引必须为 sensor.time_reversal_boundary_data（sensor_point_index，time_index）。如果以直角坐标集的形式给出 sensor.mask，则必须以相同的顺序给出边界数据。然后，在每个时间步长内使用等效的二进制传感器掩码（使用近邻插值法计算）将压力值放入计算网格。如果 sensor.mask 以传感器点的二进制矩阵形式给出，则必须使用 MATLAB 的标准列向线性矩阵索引对边界数据进行排序。如果不需要额外输入，可以用空数组[]代替源输入。

在时间反转图像重建过程中，还可以通过分配吸收参数 medium.alpha_coeff 和 medium.alpha_power，并通过设置 medium.alpha_sign = [-1, 1]来反转吸收项的符号，从而实现声衰减补偿。这将迫使传播波根据吸收参数增长，而不是衰减。然后，应为 medium.alpha_filter 指定一个滤波器（可使用 getAlphaFilter 创建），对重建进行正则化处理。

注：要运行一个使用时间反转的简单光声图像重建示例（这犯了使用相同数值参数和模型进行数据模拟和图像重建的 "inverse crime"），可将 k 波模拟返回的 sensor_data 直接传递给 sensor.time_reversal_boundary_data，并删除输入字段 source.p0 和 source.p 或将其设置为零。

Input

注：运行初值问题（例如光声正演模拟）时必须分配的最小字段用 * 标记。

kgrid*	k-Wave grid object returned by kWaveGrid containing Cartesian and k-space grid fields
kgrid.t_array*	evenly spaced array of time values [s] (set to 'auto' by kWaveGrid)
medium.sound_speed*	sound speed distribution within the acoustic medium [m/s]
medium.sound_speed_ref	reference sound speed used within the k-space operator (phase correction term) [m/s]
medium.density*	density distribution within the acoustic medium [kg/m^3]
medium.BonA	parameter of nonlinearity
medium.alpha_power	power law absorption exponent
medium.alpha_coeff	power law absorption coefficient [dB/(MHz^y cm)]
medium.alpha_mode	optional input to force either the absorption or dispersion terms in the equation of state to be excluded; valid inputs are 'no_absorption' or 'no_dispersion'
medium.alpha_filter	frequency domain filter applied to the absorption and dispersion terms in the equation of state
medium.alpha_sign	two element array used to control the sign of absorption and dispersion terms in the equation of state
	以上为medium相关
source.p0*	initial pressure within the acoustic medium
source.p	time varying pressure at each of the source positions given by source.p_mask
source.p_mask	binary matrix specifying the positions of the time varying pressure source distribution
source.p_mode	optional input to control whether the input pressure is injected as a mass source or enforced as a dirichlet boundary condition; valid inputs are 'additive' (the default) or 'dirichlet'
source.ux	time varying particle velocity in the x-direction at each of the source positions given by source.u_mask
source.uy	time varying particle velocity in the y-direction at each of the source positions given by source.u_mask
source.u_mask	binary matrix specifying the positions of the time varying particle velocity distribution
source.u_mode	optional input to control whether the input velocity is applied as a force source or enforced as a dirichlet boundary condition; valid inputs are 'additive' (the default) or 'dirichlet'
	以上为source相关
sensor.mask*	binary matrix or a set of Cartesian points where the pressure is recorded at each time-step
sensor.record	cell array of t