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 the acoustic parameters to record in the form sensor.record = {'p', 'u', ...}; valid inputs are:   'p' (acoustic pressure)   'p_max' (maximum pressure)   'p_min' (minimum pressure)   'p_rms' (RMS pressure)   'p_final' (final pressure field at all grid points)   'p_max_all' (maximum pressure at all grid points)   'p_min_all' (minimum pressure at all grid points)   'u' (particle velocity)   'u_max' (maximum particle velocity)   'u_min' (minimum particle velocity)   'u_rms' (RMS particle velocity)   'u_final' (final particle velocity field at all grid points)   'u_max_all' (maximum particle velocity at all grid points)   'u_min_all' (minimum particle velocity at all grid points)   'u_non_staggered' (particle velocity on non-staggered grid)   'I' (time varying acoustic intensity)   'I_avg' (average acoustic intensity)
sensor.record_start_index	time index at which the sensor should start recording the data specified by sensor.record (default = 1)
sensor.time_reversal_boundary_data	time varying pressure enforced as a Dirichlet boundary condition over sensor.mask
sensor.frequency_response	two element array specifying the center frequency and percentage bandwidth of a frequency domain Gaussian filter applied to the sensor_data
sensor.directivity_angle	matrix of directivity angles (direction of maximum response) for each sensor element defined in sensor.mask. The angles are in radians where 0 = max sensitivity in x direction (up/down) and pi/2 or -pi/2 = max sensitivity in y direction (left/right)
sensor.directivity_size	equivalent element size (the larger the element size the more directional the response)

注意：对于异质介质，medium.sound_speed（介质声速）和 medium.density （介质密度）必须以矩阵形式给出，尺寸与 kgrid 相同。对于均质介质，可以标量形式给出。如果介质是均质的，且不需要速度输入或输出，则不必指定 medium.density。

 Optional Inputs

可选的 "字符串"，可用于修改默认计算设置的值对。 

Input	Valid Settings	Default	Description
'CartInterp'	'linear' 'nearest'	'linear'	Interpolation mode used to extract the pressure when a Cartesian sensor mask is given. If set to 'nearest' and more than one Cartesian point maps to the same grid point, duplicated data points are discarded and sensor_data will be returned with less points than that specified by sensor.mask.
'CreateLog'	(Boolean scalar)	false	Boolean controlling whether the command line output is saved using the diary function with a date and time stamped filename.
'DataCast'	(string of data type)	'off'	String input of the data type that variables are cast to before computation. For example, setting to 'single' will speed up the computation time (due to the improved efficiency of ifftn for this data type) at the expense of a loss in precision. This variable is also useful for utilising GPU parallelisation through libraries such as the Parallel Computing Toolbox by setting 'DataCast' to 'gpuArray-single'.
'DataRecast'	(Boolean scalar)	false	Boolean controlling whether the output data is cast back to double precision. If set to false, sensor_data will be returned in the data format set using the 'DataCast' option.
'DisplayMask'	(binary matrix) or 'off'	sensor.mask	Binary matrix overlaid onto the animated simulation display. Elements set to 1 within the display mask are set to black within the display.
'HDFCompressionLevel'	(Integer scalar)	0	Compression level used for writing the input HDF5 file when using 'SaveToDisk' or kspaceFirstOrder2DC. Can be set to an integer between 0 (no compression) and 9 (maximum compression). The compression is lossless. Increasing the compression level will reduce the file size if there are portions of the medium that are homogeneous, but will also increase the time to create the HDF5 file.
'LogScale'	(Boolean scalar) or (numeric scalar)	false	Boolean controlling whether the pressure field is log compressed before display. The data is compressed by scaling both the positive and negative values between 0 and 1 (truncating the data to the given plot scale), adding a scalar value (compression factor) and then using the corresponding portion of a log10 plot for the compression (the negative parts are remapped to be negative thus the default color scale will appear unchanged). The amount of compression can be controlled by adjusting the compression factor which can be given in place of the Boolean input. The closer the compression factor is to zero, the steeper the corresponding part of the log10 plot used, and the greater the compression (the default compression factor is 0.02).
'MeshPlot'	(Boolean scalar)	false	Boolean controlling whether imagesc to plot the pressure field. When 'MeshPlot' is set to true, the default display mask is set to 'off'.
'MovieArgs'	(string cell array)	{}	Settings for VideoWriter. Parameters must be given as {'param', value, ...} pairs within a cell array (default = {}), where 'param' corresponds to a writable property of a VideoWriter object.
'MovieName'	(string)	'date-time-kspaceFirstOrder2D'	Name of the movie produced when 'RecordMovie' is set to true.
'MovieProfile'	(string)	'Uncompressed AVI'	Profile input passed to VideoWriter.
'PlotFreq'	(integer numeric scalar)	10	The number of iterations which must pass before the simulation plot is updated.
'PlotLayout'	(Boolean scalar)	false	Boolean controlling whether a four panel plot of the initial simulation layout is produced (initial pressure, sensor mask, sound speed, density).
'PlotPML'	(Boolean scalar)	true	Boolean controlling whether the perfectly matched layer is shown in the simulation plots. If set to false, the PML is not displayed.
'PlotScale'	(numeric two element vector) or 'auto'	[-1, 1]	[min, max] values used to control the scaling for imagesc (visualisation). If set to 'auto', a symmetric plot scale is chosen automatically for each plot frame.
'PlotSim'	(Boolean scalar)	true	Boolean controlling whether the simulation iterations are progressively plotted.
'PMLAlpha'	(numeric scalar or three element vector)	2	Absorption within the perfectly matched layer in Nepers per grid point.
'PMLInside'	(Boolean scalar)	true	Boolean controlling whether the perfectly matched layer is inside or outside the grid. If set to false, the input grids are enlarged by PMLSize before running the simulation.
'PMLSize'	(integer numeric scalar or three element vector)	10	Size of the perfectly matched layer in grid points. By default, the PML is added evenly to all sides of the grid, however, both PMLSize and PMLAlpha can be given as two element arrays to specify the x and y properties, respectively. To remove the PML, set the appropriate PMLAlpha to zero rather than forcing the PML to be of zero size.
'RecordMovie'	(Boolean scalar)	false	Boolean controlling whether the displayed image frames are captured and stored as a movie using VideoWriter.
'Smooth'	(Boolean scalar or three element vector)	[true, false, false]	Boolean controlling whether source.p0, medium.sound_speed, and medium.density are smoothed using smooth before computation. 'Smooth' can either be given as a single Boolean value or as a 3 element array to control the smoothing of source.p0, medium.sound_speed, and medium.density, independently.

修改方法

手动修改input_args，并在kspaceFirstOrder2D中输入input_args{:}参数。

应用示例 

 通过修改PlotLayout参数，可控制是否绘制初始模拟布局的四幅图（初始压力、传感器掩膜、声速、密度）。PlotLayout参数默认值为false，通过上方代码修改为true后可得到：

 Output

If sensor.record is not defined by the user:

sensor_data

time varying pressure recorded at the sensor positions given by sensor.mask

If sensor.record is defined by the user:

sensor_data.p	time varying pressure recorded at the sensor positions given by sensor.mask (returned if 'p' is set)
sensor_data.p_max	maximum pressure recorded at the sensor positions given by sensor.mask (returned if 'p_max' is set)
sensor_data.p_min	minimum pressure recorded at the sensor positions given by sensor.mask (returned if 'p_min' is set)
sensor_data.p_rms	rms of the time varying pressure recorded at the sensor positions given by sensor.mask (returned if 'p_rms' is set)
sensor_data.p_final	final pressure field at all grid points within the domain (returned if 'p_final' is set)
sensor_data.p_max_all	maximum pressure recorded at all grid points within the domain (returned if 'p_max_all' is set)
sensor_data.p_min_all	minimum pressure recorded at all grid points within the domain (returned if 'p_min_all' is set)
sensor_data.ux	time varying particle velocity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'u' is set)
sensor_data.uy	time varying particle velocity in the y-direction recorded at the sensor positions given by sensor.mask (returned if 'u' is set)
sensor_data.ux_max	maximum particle velocity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'u_max' is set)
sensor_data.uy_max	maximum particle velocity in the y-direction recorded at the sensor positions given by sensor.mask (returned if 'u_max' is set)
sensor_data.ux_min	minimum particle velocity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'u_min' is set)
sensor_data.uy_min	minimum particle velocity in the y-direction recorded at the sensor positions given by sensor.mask (returned if 'u_min' is set)
sensor_data.ux_rms	rms of the time varying particle velocity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'u_rms' is set)
sensor_data.uy_rms	rms of the time varying particle velocity in the y-direction recorded at the sensor positions given by sensor.mask (returned if 'u_rms' is set)
sensor_data.ux_final	final particle velocity field in the x-direction at all grid points within the domain (returned if 'u_final' is set)
sensor_data.uy_final	final particle velocity field in the y-direction at all grid points within the domain (returned if 'u_final' is set)
sensor_data.ux_max_all	maximum particle velocity in the x-direction recorded at all grid points within the domain (returned if 'u_max_all' is set)
sensor_data.uy_max_all	maximum particle velocity in the y-direction recorded at all grid points within the domain (returned if 'u_max_all' is set)
sensor_data.ux_min_all	minimum particle velocity in the x-direction recorded at all grid points within the domain (returned if 'u_min_all' is set)
sensor_data.uy_min_all	minimum particle velocity in the y-direction recorded at all grid points within the domain (returned if 'u_min_all' is set)
sensor_data.ux_non_staggered	time varying particle velocity in the x-direction recorded at the sensor positions given by sensor.mask after shifting to the non-staggered grid (returned if 'u_non_staggered' is set)
sensor_data.uy_non_staggered	time varying particle velocity in the y-direction recorded at the sensor positions given by sensor.mask after shifting to the non-staggered grid (returned if 'u_non_staggered' is set)
sensor_data.Ix	time varying acoustic intensity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'I' is set)
sensor_data.Iy	time varying acoustic intensity in the y-direction recorded at the sensor positions given by sensor.mask (returned if 'I' is set)
sensor_data.Ix_avg	average acoustic intensity in the x-direction recorded at the sensor positions given by sensor.mask (returned if 'I_avg' is set)
sensor_data.Iy_avg	average acoustic intensity in the z-direction recorded at the sensor positions given by sensor.mask (returned if 'I_avg' is set)

 应用示例1：Recording The Particle Velocity

1.Setting the acoustic variables that are recorded

默认情况下，k-Wave 的一阶模拟函数会返回每个时间步在 sensor.mask 指定位置记录的压力场。这些数据将直接分配给输出参数 sensor_data（除非使用对角定义的传感器掩码）。还可以通过设置 sensor.record 的值来控制传感器掩码记录的声学变量。所需的字段参数在单元格数组中以字符串形式列出。例如，要同时记录声压和粒子速度，sensor.record 应设置为 {'p', 'u'}。如果设置了 sensor.record 的值，模拟返回的输出 sensor_data 将被定义为一个结构体，其中记录的声学变量将作为结构体字段附加。例如，如果 sensor.record = {'p'、'p_max'、'u'}，则单个输出变量的访问方式为 sensor_data.p、sensor_data.p_max、sensor_data.ux、sensor_data.uy（其他维度类似）。其他 sensor.record 选项的完整列表请参见 kspaceFirstOrder2D。

% set the acoustic variables that are recorded
sensor.record = {'p', 'u'};

2.Defining the time array

在前面的示例中，kWaveGrid 用于定义空间离散的属性。创建 kWaveGrid 对象后，用户可以修改的唯一参数是模拟所用的时间步长和步数。这些参数由 kgrid.Nt 和 kgrid.dt 定义。如果将这两个参数设置为 "auto"（默认值），时间阵列将通过调用 kWaveGrid 类的 makeTime 方法在 kspaceFirstOrder2D 中自动计算。这将总时间设置为声波以最小声速穿过最长网格对角线所需的时间。时间步长基于 0.3 的 Courant-Friedrichs-Lewy (CFL) 数字和介质中的最大声速，即 kgrid.dt = CFL * dx_min / c0_max（有关 CFL 数字与精度和稳定性之间关系的详细讨论，请参阅 k-Wave 手册）。

除了将 kgrid.Nt 和 kgrid.dt 设置为 "自动 "外，还可以手动设置时间阵列。kWaveGrid 类有两个方法可以使用： (1) makeTime（根据上述声速和 CFL 设置时间参数）和 (2) setTime（直接设置 Nt 和 dt）。这里的时间阵列是使用 makeTime 定义的，使用的是默认的 CFL 和 6 秒的模拟时间。

% create time array
t_end = 6e-6;       % [s]
kgrid.makeTime(kgrid, medium.sound_speed, [], t_end);

创建后，可以使用 kgrid.Nt 和 kgrid.dt 查询时间点的数量和时间步长，也可以使用 kgrid.t_array 返回完整的时间数组。请注意，时间数组总是均匀分布且单调递增的。 

3.Running the simulation

模拟运行方式与之前的示例相同。在本例中，以一个小圆形声源的位置为中心定义了四个传感器点，并将传感器设置为同时记录声压和粒子速度。初始压力分布图和传感器掩膜图如下。

4.实验结果

直接使用kgrid.makeTime(kgrid, medium.sound_speed, [], t_end)会报错：

修改为

kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);

可正确运行。 

直接使用plot()绘制sensor_data.p&sensor_data.ux，得到结果明显错误：

按照官方示例，应该可以直接得到在四个传感器位置记录的压力和速度信号图如下，

但是我没有找到生成的方法，直接遍历每个传感器得到其压力随时间变化的图像，曲线的趋势和理论得到的相符：

应用示例2：保存影片文件

1.Setting the 'RecordMovie' flag

在前面的示例中，可选输入 "PlotLayout（绘图布局）"和 "PlotPML "用于更改 kspaceFirstOrder2D 的默认行为。这里使用了其他几个可选输入来将模拟动画保存到影片文件中。将 "RecordMovie（录制动画）"设置为 true 后，显示的图像帧将被保存并导出为带有日期和时间戳的影片文件名。用户也可以通过将 "MovieName "输入设置为字符串来定义文件名。影片帧通过 getframe 捕捉，如果模拟窗口前有其他窗口移动，这些帧也会被捕获。
 

% set the input arguments
input_args = {'RecordMovie', true, 'MovieName', 'example_movie'};

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

2.Controlling the movie settings

由于影片帧是模拟仿真动画的帧，因此可以通过更改动画设置来控制影片属性。例如，可以使用 "PlotScale "设置图像缩放比例，也可以通过 "PlotFreq "控制模拟图更新前的迭代次数。需要注意的是，动画使用的默认色图的零点设置为白色（参见 getColorMap），因此使用设置为 [-a, a] 的绘图缩放参数将获得最佳视觉效果。

% 设置输入参数
input_args = {..., 'PlotScale', [-2, 2], 'PlotFreq', 5, ...}；

VideoWriter 的配置文件和属性也可以通过可选输入 "MovieProfile "和 "MovieArgs "进行修改。下面是一个更改默认帧频并将影片保存为 MP4 的示例。

% 设置输入参数
input_args = {..., 'MovieProfile', 'MPEG-4', 'MovieArgs', {'FrameRate', 10}, ...}；

影片保存在与示例 m 文件相同的目录中。

请注意，默认情况下，k-Wave 会将 "MovieProfile "设置为 "Uncompressed AVI"，这意味着生成的电影文件可能会非常大。如果使用的是最新版本的 MATLAB，建议使用可选输入 "MovieProfile"，即 "MPEG-4"。另外，为了在创建电影后减小文件大小，可以使用开源软件包 Hand Brake 压缩和转换文件。