CN104036140A - 一种用于声学不均匀媒介的磁声耦合成像声压求解方法 - Google Patents

Abstract

一种用于声学不均匀媒介的磁声耦合成像声压求解方法，在求解过程中采用离散化方式，考虑了被测样本和声耦合媒介的声速和密度的不均匀性带来的声波传播过程中的反射，透射和折射问题，包括：确定声学不均匀介质中的声压求解基本方程；根据时域有限差分的稳定性条件和数值色散条件确定网格划分方式和抽样方式；对声压求解方程进行离散差分化；根据差分化基本方程编写程序，进行声压求解。本发明能够有效考虑由于媒介的声学不均匀性导致的声波在传播过程中存在复杂的反射，折射问题，进而求解出更接近真实状态的声压。各点的声压分量仅仅与和它相邻的质点振动速度有关，随着时间推进直接模拟声波在样本内部的传播及其和周围物体的相互作用。

Description

一种用于声学不均匀媒介的磁声耦合成像声压求解方法

技术领域

本发明涉及一种磁声耦合成像中的声压求解方法。特别是涉及一种用于声学不均匀媒介的磁声耦合成像声压求解方法。

背景技术

磁声耦合成像是一种结合电阻抗成像和超声成像的技术的新型生物组织电特性成像技术。其原理是：将被测生物组织置于稳恒磁场中，对生物组织给予稳恒磁场方向平行的脉冲磁场，进而在其内部产生感应电流(或者通过电极直接往生物组织内部注入脉冲电流)；被测组织内部的注入电流在稳恒磁场的作用下受洛伦兹力影响产生与激励信号同频率的位移，进而产生高频振动，向外界辐射超声波；包含物体内部的声信号在生物组织内部通过一系列复杂的反射、折射等变化传递到体表，被换能器所接收并用于重建内部电导率信息。

由磁声耦合成像的原理可知，一个电流分布为(变量上端的波浪线指代变量包含时间量)的媒介置于一个静态磁场B₀中，可以得到声压的波动方程如下：

${\▿}^2\tilde{p}(r,t)-\frac{1}{c_s^2}\·\frac{{\∂}^2\tilde{p}(r,t)}{{\∂}^{2}t}=\▿\·(\tilde{J}\×B_0)---\left(1\right)$

上式中c_s指代媒质内的声速。由该式可知，在已知稳恒磁场前提下，具有一定声学特性的媒介产生及传播的声压只与电流密度随时间变化量有关。故利用声换能器检测到的声压与能够反映沿此方向内部电流密度变化，进而反映媒介内部电特性分布。

在磁声耦合成像正问题解析过程中，生物组织内部拥有复杂的声学特性，其声学不均匀性存在于组织的各个部分；复杂的声学构造导致振动声波在组织内部的传播过程存在反射，折射，透射等现象。然而，目前的研究(如方程1)均建立在声学均匀前提下，忽视了声学参数对声传播的影响，在理想状态下求解生物组织内部的电特性，影响内部电导率特征的提取。

发明内容

本发明所要解决的技术问题是，提供一种考虑了由于声学不均匀特性带来的声传播过程中的反射，折射，透射等问题的用于声学不均匀媒介的磁声耦合成像声压求解方法。

本发明所采用的技术方案是：一种用于声学不均匀媒介的磁声耦合成像声压求解方法，在求解过程中采用离散化方式，考虑了被测样本和声耦合媒介的声速和密度的不均匀性带来的声波传播过程中的反射，透射和折射问题，具体包括如下阶段：

1)确定声学不均匀介质中的声压求解基本方程，包括如下步骤：

(1)通过磁声耦合成像声源的产生机制确定非均匀媒介中的声学波动方程；

(2)根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

2)根据时域有限差分的稳定性条件和数值色散条件确定网格划分方式和抽样方式，包括如下步骤：

(1)根据磁声耦合成像原理及激励信号特征和声学参数，在数值色散条件下确定空间网格尺寸；

(2)根据空间网格尺寸，空间维度和非均匀分布的声学参数，设定在三维空间各方向的网格尺寸，依据时域有限差分解的稳定性条件确定时间网格尺寸；

(3)根据步骤(1)和步骤(2)确定的空间网格和时间网格，设定时间和空间交替抽样方式，将声压和质点振速采用时间和空间各相隔半个步长在四维结构中交替抽样；

3)对声压求解方程进行离散差分化

根据阶段2)所述网格划分方式，对阶段1)中方程的各个变量进行离散化，并依据声传播过程中的边界条件，确定差分格式，建立差分化基本方程；

4)根据阶段3)的差分化基本方程编写程序，进行声压求解。

阶段1)所述的非均匀媒介中的声学波动方程，是基于声学物态方程、连续性方程和运动方程变换后，结合磁声耦合成像声源产生机制，得到在不均匀媒介中的声压和振速表达式，公式如下：

$\frac{\∂p}{\∂t}=-{\mathrm{\ρc}}^2(\frac{{\∂v}_x}{\∂x}+\frac{{\∂v}_y}{\∂y}+\frac{{\∂v}_z}{\∂z})+c^2\·\mathrm{\η}\left(t\right)[\▿\·(J\×B)]$

$\left\{\begin{array}{c}\frac{{\∂v}_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂x}\\ \frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂y}\\ \frac{\∂v_z}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p}{\∂z}\end{array}\right.$

其中，v_x,v_y,v_z表示质点在x,y,z方向声的振动速度，p表示声压，J表示电流密度，B表示磁通密度，ρ表示物质密度，c表示声速，η(t)表示激励信号。

阶段2)所述的磁声耦合成像采用的激励为单脉冲，所述的确定空间网格时采用的频率为2倍于原始激励的频率，表达式如下：

f_max≥2f_in

其中，f_in为激励信号的频率，即在磁声耦合成像过程中所关注的频率，f_max为空间网格划分过程中采用最大理论频率。

阶段3)中所述的差分化基本方程如下：

$\begin{array}{c}{p}^{n+1}(i,j,k)=p^n(i,j,k)-\mathrm{\ρ}(i,j,k)c^2(i,j,k)\·\mathrm{\Δt}\·\frac{1}{\mathrm{\Δs}}\·\left\{\right[v_x^{n+1/2}(i+\frac{1}{2},j,k)-v_x^{n+1/2}(i-\frac{1}{2},j,k)]\\ +[v_y^{n+1/2}(i,j+\frac{1}{2},k)-v_y^{n+1/2}(i,j-\frac{1}{2},k)]+[v_z^{n+1/2}(i,j+\frac{1}{2},k)-v_z^{n+1/2}(i,j-\frac{1}{2},k)]\}\\ +c^2(i,j,k)\·{\mathrm{\η}}^n[\▿\·(J(i,j,k)\×B(i,j,k))]\·\mathrm{\Δt}\end{array}$

$v_x^{n+1/2}(i-\frac{1}{2},j,k)=v_x^{n-1/2}(i-\frac{1}{2},j,k)-\frac{\mathrm{\Δt}}{\mathrm{\ρ}(i,j,k)\·\mathrm{\Δs}}[p^n(i,j,k)-p^n(i-1,j,k)]$

$v_y^{n+1/2}(i,j-\frac{1}{2},k)=v_y^{n-1/2}(i,j-\frac{1}{2},k)-\frac{\mathrm{\Δt}}{\mathrm{\ρ}(i,j,k)\·\mathrm{\Δs}}[p^n(i,j,k)-p^n(i,j-1,k)].$

$v_z^{n+1/2}(i,j,k-\frac{1}{2})=v_y^{n-1/2}(i,j,k-\frac{1}{2})-\frac{\mathrm{\Δt}}{\mathrm{\ρ}(i,j,k)\·\mathrm{\Δs}}[p^n(i,j,k)-p^n(i,j,k-1)]$

阶段4)中所述的程序包括如下步骤：

(1)确定磁声耦合成像模型的初始条件；

在磁声耦合成像模型中，根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

(2)确定自由空间大小及时间长度T

依据实验需要确定算法模拟的空间大小，自由空间大小要满足磁声信号传递及检测需要，并依据自由空间大小声组织及检测环境声学特性确定时间长度T；

(3)确定时间和空间网格尺寸

在确定组织及环境声学特性前提下，依据阶段2)中给出的空间网格和时间网格确定方式，确定编程过程中采用的空间和时间网格尺寸；

(4)确定初始时刻的声压

将步骤(1)中确定的初始条件代入阶段3)中所述的差分化基本方程求解得到初始时刻t1＝0的声压p；

(5)求解振速

依据t1时刻的声压，求解t2＝t1+0.5时刻的振速v_x,v_y,v_z；

(6)求解声压

依据t2时刻的振速，求解t1＝t2+0.5时刻的声压p，存储数据；

(7)循环判定

如果t1<T，转到步骤(5)，否则结束。

步骤(1)所述的初始条件包括：电流密度J，磁通密度B，物质密度ρ，声速媒介声速分布c和激励信号η(t)。

本发明的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，利用空间离散化的方式能够有效考虑由于媒介的声学不均匀性导致的声波在传播过程中存在复杂的反射，折射问题，进而求解出更接近真实状态的声压。同时，计算过程中把包含时间变量和空间变量的磁声成像基本方程转化为yee网格中的差分方程，各点的声压分量仅仅与和它相邻的质点振动速度有关，随着时间推进直接模拟声波在样本内部的传播及其和周围物体的相互作用。

附图说明

图1是本发明方法的流程图；

图2是本发明方法中空间抽样方式示意图；

图3是本发明的算法编程实现流程图。

具体实施方式

下面结合实施例和附图对本发明的一种用于声学不均匀媒介的磁声耦合成像声压求解方法做出详细说明。

本发明的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，在求解过程中采用离散化方式，考虑了被测样本和声耦合媒介的声速和密度的不均匀性带来的声波传播过程中的反射，透射和折射问题，具体包括如下阶段：

1)确定声学不均匀介质中的声压求解基本方程，包括如下步骤：

(1)通过磁声耦合成像声源的产生机制确定非均匀媒介中的声学波动方程；

所述的非均匀媒介中的声学波动方程，是根据磁声耦合成像原理及基于声学物态方程、连续性方程和运动方程变换后，结合磁声耦合成像声源产生机制，得到在不均匀媒介中的声压和振速表达式，公式如下：

$\frac{\&PartialD;p}{\&PartialD;t}=-{\mathrm{\&rho;c}}^2(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_y}{\&PartialD;y}+\frac{{\&PartialD;v}_z}{\&PartialD;z})+c^2\&CenterDot;\mathrm{\&eta;}\left(t\right)[\&dtri;\&CenterDot;(J\&times;B)]---\left(2\right)$

$\left\{\begin{array}{c}\frac{{\&PartialD;v}_x}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;x}\\ \frac{\&PartialD;v_y}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;y}\\ \frac{\&PartialD;v_z}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;z}\end{array}\right.---\left(3\right)$

其中，v_x,v_y,v_z表示质点在x,y,z方向声的振动速度，p表示声压，J表示电流密度，B表示磁通密度，ρ表示物质密度，c表示声速，η(t)表示激励信号。上述方程组中，确定初始条件电流密度J，磁通密度B，物质密度ρ，声速媒介声速分布c和激励信号η(t)，求解声压分布p。

(2)根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

2)根据时域有限差分的稳定性条件和数值色散条件确定网格划分方式和抽样方式，包括如下步骤：

(1)根据磁声耦合成像原理及激励信号特征和声学参数，在数值色散条件下确定空间网格尺寸，

所述的磁声耦合成像采用的激励为单脉冲，所述的确定空间网格时采用的频率为2倍于原始激励的频率，空间离散网格尺寸的确定，除了受到资源限制外，还需要考虑的主要因素是色散误差，假定磁声耦合成像实验中激励的频率为f_in，考虑磁声耦合成像原理，激励信号为单脉冲，故系统中最大频率f_max应当满足下式：

f_max≥2f_in (4)

其中，f_in为激励信号的频率，即在磁声耦合成像过程中所关注的频率，f_max为空间网格划分过程中采用最大理论频率，

故得到对应的最小波长λ_min。设定在笛卡尔坐标系中3个方向的网格尺寸相等，均为Δs；故为满足色散误差要求，空间网格符合如下：

$\mathrm{\&Delta;s}\&le;\frac{{\mathrm{\&lambda;}}_{\min }}{N};(N\&GreaterEqual;10)---\left(5\right)$

非均匀媒介构成的磁声成像空间中，选用courant稳定条件，得到时域有限差分的时间空间关系如下： $\mathrm{\&Delta;t}\&le;\frac{1}{c_{\max }\sqrt{3*{\left(\frac{1}{\mathrm{\&Delta;s}}\right)}^2}}---\left(6\right)$

其中，c_max为磁声耦合成像实验中的最大声速。具体抽样方式如图2所示，采用声压和振速在时间和空间上分布相差半个步长进行交替抽样。

(2)根据空间网格尺寸，空间维度和非均匀分布的声学参数，设定在三维空间各方向的网格尺寸，依据时域有限差分解的稳定性条件确定时间网格尺寸；

(3)根据步骤(1)和步骤(2)确定的空间网格和时间网格，设定时间和空间交替抽样方式，将声压和质点振速采用时间和空间各相隔半个步长在四维结构中交替抽样；

图2是本发明的空间抽样方式。使用笛卡尔网格划分法对声传播过程中质点振动分量和声压分量在空间和时间上采取交替抽样的离散方式。依据yee元胞理论，我们将质点振动分量和声压在空间轴上按如图2所示放置。并把质点振速ν和声压P在时间上相差半步进行计算，进而实现在空间加时间的四维模型中实现交替抽样。

3)对声压求解方程进行离散差分化

根据阶段2)所述网格划分方式，对阶段1)中方程的各个变量进行离散化，并依据声传播过程中的边界条件，确定差分格式，建立差分化基本方程，

所述的差分化基本方程如下：

$\begin{array}{c}{p}^{n+1}(i,j,k)=p^n(i,j,k)-\mathrm{\&rho;}(i,j,k)c^2(i,j,k)\&CenterDot;\mathrm{\&Delta;t}\&CenterDot;\frac{1}{\mathrm{\&Delta;s}}\&CenterDot;\left\{\right[v_x^{n+1/2}(i+\frac{1}{2},j,k)-v_x^{n+1/2}(i-\frac{1}{2},j,k)]\\ +[v_y^{n+1/2}(i,j+\frac{1}{2},k)-v_y^{n+1/2}(i,j-\frac{1}{2},k)]+[v_z^{n+1/2}(i,j+\frac{1}{2},k)-v_z^{n+1/2}(i,j-\frac{1}{2},k)]\}\\ +c^2(i,j,k)\&CenterDot;{\mathrm{\&eta;}}^n[\&dtri;\&CenterDot;(J(i,j,k)\&times;B(i,j,k))]\&CenterDot;\mathrm{\&Delta;t}\end{array}---\left(7\right)$

$v_x^{n+1/2}(i-\frac{1}{2},j,k)=v_x^{n-1/2}(i-\frac{1}{2},j,k)-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i-1,j,k)]$

$v_y^{n+1/2}(i,j-\frac{1}{2},k)=v_y^{n-1/2}(i,j-\frac{1}{2},k)-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i,j-1,k)]---\left(8\right).$

$v_z^{n+1/2}(i,j,k-\frac{1}{2})=v_y^{n-1/2}(i,j,k-\frac{1}{2})-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i,j,k-1)]$

4)根据阶段3)的差分化基本方程编写程序，进行声压求解，所述的程序包括如下步骤：

(1)确定磁声耦合成像模型的初始条件，所述的初始条件包括：电流密度J，磁通密度B，物质密度ρ，声速媒介声速分布c和激励信号η(t)；

在磁声耦合成像模型中，根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

(2)确定自由空间大小及时间长度T

依据实验需要确定算法模拟的空间大小，自由空间大小要满足磁声信号传递及检测需要，并依据自由空间大小声组织及检测环境声学特性确定时间长度T；

(3)确定时间和空间网格尺寸

在确定组织及环境声学特性前提下，依据阶段2)中给出的空间网格和时间网格确定方式，确定编程过程中采用的空间和时间网格尺寸；

(4)确定初始时刻的声压

将步骤(1)中确定的初始条件代入阶段3)中所述的差分化基本方程求解得到初始时刻t1＝0的声压p；

(5)求解振速

依据t1时刻的声压，求解t2＝t1+0.5时刻的振速v_x,v_y,v_z；

(6)求解声压

依据t2时刻的振速，求解t1＝t2+0.5时刻的声压p，存储数据；

(7)循环判定

如果t1<T，转到步骤(5)，否则结束。

Claims (6)

1.一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，在求解过程中采用离散化方式，考虑了被测样本和声耦合媒介的声速和密度的不均匀性带来的声波传播过程中的反射，透射和折射问题，具体包括如下阶段：

1)确定声学不均匀介质中的声压求解基本方程，包括如下步骤：

(1)通过磁声耦合成像声源的产生机制确定非均匀媒介中的声学波动方程；

(2)根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

2)根据时域有限差分的稳定性条件和数值色散条件确定网格划分方式和抽样方式，包括如下步骤：

(1)根据磁声耦合成像原理及激励信号特征和声学参数，在数值色散条件下确定空间网格尺寸；

(2)根据空间网格尺寸，空间维度和非均匀分布的声学参数，设定在三维空间各方向的网格尺寸，依据时域有限差分解的稳定性条件确定时间网格尺寸；

(3)根据步骤(1)和步骤(2)确定的空间网格和时间网格，设定时间和空间交替抽样方式，将声压和质点振速采用时间和空间各相隔半个步长在四维结构中交替抽样；

3)对声压求解方程进行离散差分化

根据阶段2)所述网格划分方式，对阶段1)中方程的各个变量进行离散化，并依据声传播过程中的边界条件，确定差分格式，建立差分化基本方程；

4)根据阶段3)的差分化基本方程编写程序，进行声压求解。

2.根据权利要求1所述的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，阶段1)所述的非均匀媒介中的声学波动方程，是基于声学物态方程、连续性方程和运动方程变换后，结合磁声耦合成像声源产生机制，得到在不均匀媒介中的声压和振速表达式，公式如下：

$\frac{\&PartialD;p}{\&PartialD;t}=-{\mathrm{\&rho;c}}^2(\frac{{\&PartialD;v}_x}{\&PartialD;x}+\frac{{\&PartialD;v}_y}{\&PartialD;y}+\frac{{\&PartialD;v}_z}{\&PartialD;z})+c^2\&CenterDot;\mathrm{\&eta;}\left(t\right)[\&dtri;\&CenterDot;(J\&times;B)]$

$\left\{\begin{array}{c}\frac{{\&PartialD;v}_x}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;x}\\ \frac{\&PartialD;v_y}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;y}\\ \frac{\&PartialD;v_z}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;p}{\&PartialD;z}\end{array}\right.$

其中，v_x,v_y,v_z表示质点在x,y,z方向声的振动速度，p表示声压，J表示电流密度，B表示磁通密度，ρ表示物质密度，c表示声速，η(t)表示激励信号。

3.根据权利要求1所述的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，阶段2)所述的磁声耦合成像采用的激励为单脉冲，所述的确定空间网格时采用的频率为2倍于原始激励的频率，表达式如下：

f_max≥2f_in

其中，f_in为激励信号的频率，即在磁声耦合成像过程中所关注的频率，f_max为空间网格划分过程中采用最大理论频率。

4.根据权利要求1所述的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，阶段3)中所述的差分化基本方程如下：

$\begin{array}{c}{p}^{n+1}(i,j,k)=p^n(i,j,k)-\mathrm{\&rho;}(i,j,k)c^2(i,j,k)\&CenterDot;\mathrm{\&Delta;t}\&CenterDot;\frac{1}{\mathrm{\&Delta;s}}\&CenterDot;\left\{\right[v_x^{n+1/2}(i+\frac{1}{2},j,k)-v_x^{n+1/2}(i-\frac{1}{2},j,k)]\\ +[v_y^{n+1/2}(i,j+\frac{1}{2},k)-v_y^{n+1/2}(i,j-\frac{1}{2},k)]+[v_z^{n+1/2}(i,j+\frac{1}{2},k)-v_z^{n+1/2}(i,j-\frac{1}{2},k)]\}\\ +c^2(i,j,k)\&CenterDot;{\mathrm{\&eta;}}^n[\&dtri;\&CenterDot;(J(i,j,k)\&times;B(i,j,k))]\&CenterDot;\mathrm{\&Delta;t}\end{array}$

$v_x^{n+1/2}(i-\frac{1}{2},j,k)=v_x^{n-1/2}(i-\frac{1}{2},j,k)-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i-1,j,k)]$

$v_y^{n+1/2}(i,j-\frac{1}{2},k)=v_y^{n-1/2}(i,j-\frac{1}{2},k)-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i,j-1,k)].$

$v_z^{n+1/2}(i,j,k-\frac{1}{2})=v_z^{n-1/2}(i,j,k-\frac{1}{2})-\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}(i,j,k)\&CenterDot;\mathrm{\&Delta;s}}[p^n(i,j,k)-p^n(i,j,k-1)]$

5.根据权利要求1所述的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，阶段4)中所述的程序包括如下步骤：

(1)确定磁声耦合成像模型的初始条件；

在磁声耦合成像模型中，根据成像目标和实验环境确定声学参数，依据实验条件确定稳恒磁场磁通密度和激励信号，利用电磁场仿真软件求解得到电流密度分布；

(2)确定自由空间大小及时间长度T

依据实验需要确定算法模拟的空间大小，自由空间大小要满足磁声信号传递及检测需要，并依据自由空间大小声组织及检测环境声学特性确定时间长度T；

(3)确定时间和空间网格尺寸

在确定组织及环境声学特性前提下，依据阶段2)中给出的空间网格和时间网格确定方式，确定编程过程中采用的空间和时间网格尺寸；

(4)确定初始时刻的声压

将步骤(1)中确定的初始条件代入阶段3)中所述的差分化基本方程求解得到初始时刻t1＝0的声压p；

(5)求解振速

依据t1时刻的声压，求解t2＝t1+0.5时刻的振速v_x,v_y,v_z；

(6)求解声压

依据t2时刻的振速，求解t1＝t2+0.5时刻的声压p，存储数据；

(7)循环判定

如果t1<T，转到步骤(5)，否则结束。

6.根据权利要求5所述的一种用于声学不均匀媒介的磁声耦合成像声压求解方法，其特征在于，步骤(1)所述的初始条件包括：电流密度J，磁通密度B，物质密度ρ，声速媒介声速分布c和激励信号η(t)。