CN104688224A - 一种应用于声学非均匀媒介磁声耦合成像重建方法 - 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)在扫描过程中采用平面线阵换能器阵列采集样本内部发出的磁声信号；

2)求解换能器之间的声穿越时间，包括：

(1)设定换能器检测的声压满足的条件；(2)进行声压补偿；(3)求解对称位置换能器所检测声压的相关系数；(4)求解两个轴线相互重叠的换能器聚焦区轴线部位的声穿越时间；(5)求解得到整个扫描过程的声穿越时间；

3)求解被测样本的声速分布图，包括：

(1)建立对称位置换能器连线的声速平均矩阵；(2)修正平均声速矩阵；(3)划分各单元网格；(4)求解声穿越时间；(5)对结果再次修正之后得到每个网格实际声速；

4)求解原始声源分布，包括：

(1)声速重排及声速空间扩展；(2)反转声源求解及声压转置；(3)声压数据重采样；(4)求解原始声源；(5)提取样本所在区域声源分布。

阶段1)在扫描过程中每一个换能器都要满足如下条件：

(1)换能器扫描过中完成一次扫描后进行旋转，旋转次数要满足2M-1次；

(2)每一个角度的聚焦区域移动次数为N，移动间距Δd保持不变，且满足N×Δd≥R，

其中R为能容纳样本的最小正方形区域的边长，M和N是大于等于1的整数；

(3)扫描过程中，换能器旋转半径r要满足：

$r\≥\sqrt{2}R;$

(4)处于对称位置的换能器，在聚焦区轴线重叠的区域内每一个点满足如下条件：

L₁+L₂≈2r

其中L₁为所述点到左边换能器的距离，L₂为所述点到右边换能器的距离，r为换能器扫描过程中的旋转半径。

阶段2)所述的求解换能器之间的声穿越时间，包括：

(1)设定换能器检测的声压满足如下条件：

$\tilde{p}(r_0,t)=-\frac{1}{4\mathrm{\π}}{\∫\∫}_{\mathrm{\Ω}}\frac{\▿\·(J\left(r^{`}\right)\×B\left(r^{`}\right))}{|r_0-r^{`}|}\mathrm{d\Ω}{|}_{T(r_0,r^{`})=t}$

其中r₀为换能器聚焦区轴线的延长线与换能器的交点，t为换能器接收的时间轴，为换能器所接收的声压，J为感应电流密度，B为稳恒磁场强度，r`为样本上某一点的位置，T(r₀,-r₀)为重叠的聚焦区轴线间的声传播时间。

(2)进行声压补偿

定义换能器接收信号的修正值S(r₀,t)满足：

$S(r_0,t)=\tilde{p}(r_0,t)\·t$

(3)求解对称位置换能器所检测的声压的相关系数

相对于位置r₀的换能器接收的声压修正后为S(r₀,t)，设定与r₀位置对称的信号采集点位置为-r₀，故所述信号采集点-r₀接收的声压修正后为S(-r₀,t)，考虑声的传播特性，位置r₀和位置-r₀的声压信号修正值在时间和幅度上均具有高度相关，定义满足：

$R_{r_0}\left(t\right)={\∫}_{-\∞}^{\∞}S(r_0,\mathrm{\τ})S(-r_0,-(\mathrm{\τ}-t))\mathrm{d\τ}$

为聚焦区域轴线相互重叠的一组换能器所接收信号的相关系数；

(4)求解两个换能器聚焦区轴线相互重叠部位的声穿越时间

由相关性原理及卷积定理可以得到，当t＝T(r₀,-r₀)时，式T(r₀,-r₀)可以得到最大值，故由下式得到两个换能器聚焦区轴线相互重叠部位的声穿越时间

$T(r_0,-r_0)=\arg \max \left(R_{r_0}\left(t\right)\right);$

(5)求解得到整个扫描过程的声穿越时间

对阶段1)扫描所得到的M×N对重叠的聚焦区轴线位置，依次重复步骤(1)～步骤(4)，得到M×N组穿越时间。

阶段3)所述的求解被测样本的声速分布图，具体包括：

(1)建立对称位置换能器连线的传播时间矩阵

对于扫描过程中换能器动态聚焦区的位置分布，定义传播时间矩阵T_NM的任何一个元素要满足如下公式：

T_i＝T(r_i,-r_i)；1≤i≤N×M

其中r_i,-r_i代表第i对聚焦区轴线延长线与换能器交点所处的位置；

(2)修正传播时间矩阵T_NM，得到声穿越时间T'_MN

设定耦合剂的所处区域作为背景，背景声速为c₀，对于扫描过程中换能器的动态聚焦区位置分布，得到对称位置换能器之间的距离等同于旋转直径，故排除背景声速c₀影响后的声穿越时间T'_MN的元素满足：

$T_{\mathrm{MN}}^{\′}=T_{\mathrm{MN}}-\frac{2r}{c_0}\·A_{\mathrm{MN}}$

A_MN为元素全为1的M×N矩阵；

由上式得到声穿越时间T'_MN，T'_MN是被测样本被换能器聚焦区轴线穿过部位的声穿越时间再减去相同距离耦合剂的声穿越时间；

(3)划分各单元网格

将一个感兴趣的包含被测样本的扫描空间的区域分成n×n的网格，网格大小为Δs，每个网格的声穿越时间为t_j；1≤j≤n²，网格大小Δs应当满足如下：

$\mathrm{\Δs}\≤\frac{c_{\min }}{N{f}_{\max }};N\≥10$

上式中f_max为磁声耦合成像自由空间中的最大声波频率，满足：

f_max≥2f_in

f_in为磁声耦合成像的激励信号频率，C_min为磁声耦合系统中的最小声速；

(4)求解声穿越时间

根据步骤(3)的网格，按如下公式进行迭代重建，求解每一个网格的声穿越时间值t_j：

$t_j^{(k+1)}=t_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}\frac{T_i^{\′}-\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}{t}_m}{\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}^2};1\≤i\≤N\×M,1\≤j\≤n^2;$

k为迭代次数，λ为松弛因子，a为权重因子，当网格m对声穿越时间t_j有贡献时取值为1，否则为0；

(5)对结果再次修正之后得到每个网格实际声速

步骤(4)求解得到的声穿越时间为去除背景声速之后的声穿越时间，对结果再次修正之后得到每个网格实际声速c_j；1≤j≤n²满足：

$c_j=\frac{\mathrm{\Δs}}{t_j}+c_0;1\≤j\≤n^2.$

步骤(4)所述的迭代重建步骤包括：

a)给第j个网格的声穿越时间t_j赋初值：

$t_j^{\left(k\right)}=t_j^{\left(0\right)},1\≤j\≤n^2$

式中，初值选取依据人体组织的平均声速和耦合剂的背景声速，满足：

$t_j^{\left(0\right)}=\frac{\mathrm{\Δs}}{{\stackrel{\‾}{c}}_p}-\frac{\mathrm{\Δs}}{c_0}$

其中为人体组织的平均声速；

b)计算第i对换能器间的声穿越时间修正后的估计值：

$\widetilde{T_i^{\′}}=\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^{\left(k\right)}$

c)计算误差：

$\mathrm{\Δi}=T_i^{\′}-\widetilde{T_i^{\′}}$

d)计算第j个未知量的修正值：

$C_{\mathrm{ij}}^V=\mathrm{\Δi}\frac{a_{\mathrm{ij}}}{\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}$

e)对t_j的值进行修正：

$T_j^{(k+1)}=T_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}{C}_{\mathrm{ij}}^V$

f)k＝k+1；重复步骤b)到f)，直到完成所有声穿越时间方程，即完成一轮迭代；

g)以上一轮迭代作为初值，重复b)到f)进行新一轮迭代，直到得到符合收敛要求的结果。

阶段4)所述的求解原始声源分布，具体包括：

(1)声速重排及声速空间扩展

利用阶段3)所求得到的每个网格实际声速c_j；1≤j≤n²，进行声速重排，得到n×n的声速分布矩阵将所述声速分布矩阵向外扩充为一个m×m；m>n的自由空间，声速分布为C_m×m；扩充部分的网格的声速为c₀；

(2)反转声源求解及声压转置

对扫描过程中换能器接收到的声压信号为p(r,t)，其中r表示空间位置；依据时间反演理论，将采集的声压信号反转作为声源，定义用于反演的声源项为q(r,t)满足：

q(r,t)＝p(r,T_d-t)；1≤i≤N×M

T_d表示检测过程中，空间每一个点的检测时间，故得基于时间反演的波动方程如下：

${\▿}^2{p}_0(r,t)-\frac{1}{c{\left(r\right)}^2}\frac{{\∂}^2{p}_0(r,t)}{\∂t^2}=q(r,t)$

c(r)为声速分布C_m×m圆柱坐标形式，p₀(r,T_d)为需要求解的原始声源分布；

(3)声压数据重采样

定义时间间隔Δt： $\mathrm{\Δt}\≤\frac{1}{c_{\max }\sqrt{3*{\left(\frac{1}{\mathrm{\Δs}}\right)}^2}}$

其中，c_max为磁声耦合成像实验中的最大声速，

依据时间间隔Δt对反演的声源项q(r,t)进行重采样，采样结果用反演声源Q(t)表示：

$Q\left(t\right)=q^l=q\left(\mathrm{l\Δt}\right);1\≤l\≤\frac{T_d}{\mathrm{\Δt}};$

(4)时间和空间交替迭代求解原始声源

依据声波波动方程，结合磁声耦合成像声源产生机制，得到不均匀媒介中的声压和振动速度的表达式，如下：

$\frac{\∂p_0}{\∂t}=-\mathrm{\ρ}{C_{m\×m}}^2(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y})+q$

$\left\{\begin{array}{c}\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂x}\\ \frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂y}\end{array}\right.$

其中，v_x,v_y表示质点二维平面x,y方向的振动速度，p₀表示声压，ρ表示物质密度，c表示声速，t表示波动持续时间。采用时域有限差分的方法，在时间和空间上交替离散上述表达式，利用反演声源Q求解下一个时刻的质点振动速度，并利用质点振动速度求解声压，依此交替；

(5)提取样本所在区域声源分布

依据自由空间大小，将原始声源分布p₀(r,T_d)，以自由空间中心为基准点，半径为r之外的所有区域声压置零，所得到的声压值P₀即为初始声源分布。

步骤(4)所述的不均匀媒介中的声压和振动速度的表达式具体求解过程如下：

a)确定所述表达式求解的初始条件，所述的初始条件包括：磁通密度B，物质密度ρ，

声速媒介声速分布C_m×m和采集的声压信号Q(lΔt)；

b)离散化不均匀媒介中的声压和振动速度的表达式

$\begin{array}{c}{p}_0^{l+1}(i,j)=p_0^l(i,j)-\mathrm{\ρ}(i,j)c^2(i,j)\·\mathrm{\Δt}\·\frac{1}{\mathrm{\Δs}}\·\left\{\right[v_x^{l+1/2}(i+\frac{1}{2},j)-v_x^{l+1/2}(i-\frac{1}{2},j)]\\ +[v_y^{l+1/2}(i,j+\frac{1}{2})-v_y^{l+1/2}(i,j-\frac{1}{2})]\}+c^2(i,j)\·q^l\·\mathrm{\Δt}\end{array}$

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

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

c)确定初始时刻的声压

将过程a)中确定的初始条件代入过程b)中所给出的声压方程中，求解得到初始时刻t_a＝0的声压p₀；

d)将过程c)所求t_a时刻的声压p₀代入过程b)中所给出的振动速度方程中，求解时刻的振速v_x,v_y；

e)依据过程d)所求t_b时刻的振速v_x,v_y，将换能器检测的反演声源Q代入过程b)中所给出的声压方程中，求解时刻的声压p，存储数据；

f)如果t_a＝T_d，转到过程d)，否则结束；

当时t_a＝T_d时，依据上述步骤所求解得到的声压分布p(r,T_d)，即包含原始声源分布。

本发明的一种应用于声学非均匀媒介磁声耦合成像重建方法，可以有效消除样本的声学特性不均匀带来的重建误差，同时也能够通过信号的相关性得到样本的声速分布，为声源的精确重建提供了必要的条件。

附图说明

图1是本发明应用于声学非均匀媒介磁声耦合成像重建方法的流程图；

图中：1、脉冲磁场方向，2、稳恒磁场方向，3、旋转的圆周轨道，4、被测样本，5、换能器，6、聚焦区移动方向，7、动态聚焦区

图2是本发明方法中采用的扫描方式示意图；

图3是本发明方法中声穿越时间求解流程图；

图4是本发明方法中声速分布图求解流程图；

图5是本发明方法中各网格声穿越时间求解流程图；

图6是本发明方法中初始声源求解流程图；

图7是本发明方法中基于时域有限差分求解初始声源流程图；

图8是本发明方法中采用线性换能器阵列采集声压信号的示意图

图中：1、线性换能器阵列，2、被测样本，3、换能器轴线及其延长线，4线性换能器阵列的声场示意图，5、接收换能器示意图，6、聚焦区轴线重叠时对称换能器检测到的声压信号。

具体实施方式

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

本发明的一种应用于声学非均匀媒介磁声耦合成像重建方法，是在信号采集过程中使用线阵换能器阵列采集相互对称的声压数据，利用对称位置换能器采集数据之间的相关性求解声穿越时间，结合代数迭代重建算法重建样本声速分布，最后利用所求声速分布求解样本初始声源分布，本发明的扫描方式采如图2所示。

本发明的一种应用于声学非均匀媒介磁声耦合成像重建方法具体如图1所示，包括如下阶段：

1)在扫描过程中采用平面线阵换能器阵列采集样本内部发出的磁声信号；

在扫描过程中每一个换能器都要满足如下条件：

(1)换能器扫描过中完成一次扫描后进行旋转，旋转次数要满足2M-1次；

(2)每一个角度的聚焦区域移动次数为N，移动间距Δd保持不变，且满足N×Δd≥R，

其中R为能容纳样本的最小正方形区域的边长，M和N是大于等于1的整数；

(3)扫描过程中，换能器旋转半径r要满足：

$r\≥\sqrt{2}R;$

(4)处于对称位置的换能器，在聚焦区轴线重叠的区域内每一个点满足如下条件：

L₁+L₂≈2r

其中L₁为所述点到左边换能器的距离，L₂为所述点到右边换能器的距离，r为换能器扫描过程中的旋转半径。

2)求解换能器之间的声穿越时间，如图3所示，包括：

(1)依据声传播理论，对于一个点源波随传播距离而衰减，在磁声偶成像系统中，设定换能器检测的声压满足如下条件：

$\tilde{p}(r_0,t)=-\frac{1}{4\mathrm{\π}}{\∫\∫}_{\mathrm{\Ω}}\frac{\▿\·(J\left(r^{`}\right)\×B\left(r^{`}\right))}{|r_0-r^{`}|}\mathrm{d\Ω}{|}_{T(r_0,r^{`})=t}$

其中r₀为换能器聚焦区轴线的延长线与换能器的交点，t为换能器接收的时间轴，为换能器所接收的声压，J为感应电流密度，B为稳恒磁场强度，r`为样本上某一点的位置，T(r₀,-r₀)为重叠的聚焦区轴线间的声传播时间。

(2)进行声压补偿

声传播过程中，声压随传播距离的增大而衰减，对于放置在对称位置的换能器，由于聚焦区域内的点到两个换能器的距离不一样，故相同的声源发出的信号传递到换能器之后表现出不同大小的声压，为确保换能器所检测声压之间的相关性，采用如下方式对声压进行补偿，定义换能器接收信号的修正值S(r₀,t)满足：

$S(r_0,t)=\tilde{p}(r_0,t)\·t$

(3)求解对称位置换能器所检测的声压的相关系数

相对于位置r₀的换能器接收的声压修正后为S(r₀,t)，设定与r₀位置对称的信号采集点位置为-r₀，故所述信号采集点-r₀接收的声压修正后为S(-r₀,t)，考虑声的传播特性，位置r₀和位置-r₀的声压信号修正值在时间和幅度上均具有高度相关，如图8所示，定义满足：

$R_{r_0}\left(t\right)={\∫}_{-\∞}^{\∞}S(r_0,\mathrm{\τ})S(-r_0,-(\mathrm{\τ}-t))\mathrm{d\τ}$

为聚焦区域轴线相互重叠的一组换能器所接收信号的相关系数。

(4)求解两个换能器聚焦区轴线相互重叠部位的声穿越时间

由相关性原理及卷积定理可以得到，当t＝T(r₀,-r₀)时，式T(r₀,-r₀)可以得到最大值，故由下式得到两个换能器聚焦区轴线相互重叠部位的声穿越时间

$T(r_0,-r_0)=\arg \max \left(R_{r_0}\left(t\right)\right);$

(5)求解得到整个扫描过程的声穿越时间

对阶段1)扫描所得到的M×N对重叠的聚焦区轴线位置，依次重复步骤(1)～步骤(4)，得到M×N组穿越时间。

3)求解被测样本的声速分布图，如图4所示，包括：

(1)建立对称位置换能器连线的传播时间矩阵

对于扫描过程中换能器动态聚焦区的位置分布，定义传播时间矩阵T_NM的任何一个元素要满足如下公式：

T_i＝T(r_i,-r_i)；1≤i≤N×M

其中r_i,-r_i代表第i对聚焦区轴线延长线与换能器交点所处的位置；

(2)修正传播时间矩阵T_NM，得到声穿越时间T'_MN

设定耦合剂的所处区域作为背景，背景声速为c₀，对于扫描过程中换能器的动态聚焦区位置分布，得到对称位置换能器之间的距离等同于旋转直径，故排除背景声速c₀影响后的声穿越时间T'_MN的元素满足：

$T_{\mathrm{MN}}^{\′}=T_{\mathrm{MN}}-\frac{2r}{c_0}\·A_{\mathrm{MN}}$

A_MN为元素全为1的M×N矩阵；

由上式得到声穿越时间T'_MN，T'_MN是被测样本被换能器聚焦区轴线穿过部位的声穿越时间再减去相同距离耦合剂的声穿越时间；

(3)划分各单元网格

将一个感兴趣的包含被测样本的扫描空间的区域分成n×n的网格，网格大小为Δs，每个网格的声穿越时间为t_j；1≤j≤n²，网格大小Δs应当满足如下：

$\mathrm{\Δs}\≤\frac{c_{\min }}{N{f}_{\max }};N\≥10$

上式中f_max为磁声耦合成像自由空间中的最大声波频率，满足：

f_max≥2f_in

f_in为磁声耦合成像的激励信号频率，C_min为磁声耦合系统中的最小声速；

(4)求解声穿越时间，如图5所示：

根据步骤(3)的网格，按如下公式进行迭代重建，求解每一个网格的声穿越时间值t_j：

$t_j^{(k+1)}=t_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}\frac{T_i^{\′}-\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}{t}_m}{\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}^2};1\≤i\≤N\×M,1\≤j\≤n^2;$

k为迭代次数，λ为松弛因子，a为权重因子，当网格m对声穿越时间t_j有贡献时取值为1，否则为0。

所述的迭代重建步骤包括：

a)给第j个网格的声穿越时间t_j赋初值：

$t_j^{\left(k\right)}=t_j^{\left(0\right)},1\≤j\≤n^2$

式中，初值选取依据人体组织的平均声速和耦合剂的背景声速，满足：

$t_j^{\left(0\right)}=\frac{\mathrm{\Δs}}{{\stackrel{\‾}{c}}_p}-\frac{\mathrm{\Δs}}{c_0}$

其中为人体组织的平均声速；

b)计算第i对换能器间的声穿越时间修正后的估计值：

$\widetilde{T_i^{\′}}=\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^{\left(k\right)}$

c)计算误差：

$\mathrm{\Δi}=T_i^{\′}-\widetilde{T_i^{\′}}$

d)计算第j个未知量的修正值：

$C_{\mathrm{ij}}^V=\mathrm{\Δi}\frac{a_{\mathrm{ij}}}{\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}$

e)对t_j的值进行修正：

$T_j^{(k+1)}=T_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}{C}_{\mathrm{ij}}^V$

f)k＝k+1；重复步骤b)到f)，直到完成所有声穿越时间方程，即完成一轮迭代；

g)以上一轮迭代作为初值，重复b)到f)进行新一轮迭代，直到得到符合收敛要求的结果。

(5)对结果再次修正之后得到每个网格实际声速

步骤(4)求解得到的声穿越时间为去除背景声速之后的声穿越时间，对结果再次修正之后得到每个网格实际声速c_j；1≤j≤n²满足：

$c_j=\frac{\mathrm{\Δs}}{t_j}+c_0;1\≤j\≤n^2.$

4)求解原始声源分布，求解流程如图6所示，包括：

(1)声速重排及声速空间扩展

利用阶段3)所求得到的每个网格实际声速c_j；1≤j≤n²，进行声速重排，得到n×n的声速分布矩阵将所述声速分布矩阵向外扩充为一个m×m；m>n的自由空间，声速分布为C_m×m；扩充部分的网格的声速为c₀；

(2)反转声源求解及声压转置

对扫描过程中换能器接收到的声压信号为p(r,t)，其中r表示空间位置；依据时间反演理论，将采集的声压信号反转作为声源，定义用于反演的声源项为q(r,t)满足：

q(r,t)＝p(r,T_d-t)；1≤i≤N×M

T_d表示检测过程中，空间每一个点的检测时间，故得基于时间反演的波动方程如下：

${\▿}^2{p}_0(r,t)-\frac{1}{c{\left(r\right)}^2}\frac{{\∂}^2{p}_0(r,t)}{\∂t^2}=q(r,t)$

c(r)为声速分布C_m×m圆柱坐标形式，p₀(r,T_d)为需要求解的原始声源分布；

(3)声压数据重采样

定义时间间隔Δt

$\mathrm{\Δt}\≤\frac{1}{c_{\max }\sqrt{3*{\left(\frac{1}{\mathrm{\Δs}}\right)}^2}}$

其中，c_max为磁声耦合成像实验中的最大声速，

依据时间间隔Δt对反演的声源项q(r,t)进行重采样，采样结果用反演声源Q(t)表示：

$Q\left(t\right)=q^l=q\left(\mathrm{l\Δt}\right);1\≤l\≤\frac{T_d}{\mathrm{\Δt}};$

(4)时间和空间交替迭代求解原始声源，如图7所示

依据声波波动方程，结合磁声耦合成像声源产生机制，得到不均匀媒介中的声压和振动速度的表达式，如下：

$\frac{\∂p_0}{\∂t}=-\mathrm{\ρ}{C_{m\×m}}^2(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y})+q$

$\left\{\begin{array}{c}\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂x}\\ \frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂y}\end{array}\right.$

其中，v_x,v_y表示质点二维平面x,y方向的振动速度，p₀表示声压，ρ表示物质密度，c表示声速，t表示波动持续时间。采用时域有限差分的方法，在时间和空间上交替离散上述表达式，利用反演声源Q求解下一个时刻的质点振动速度，并利用质点振动速度求解声压，依此交替；

所述的不均匀媒介中的声压和振动速度的表达式具体求解过程如下：

a)确定所述表达式求解的初始条件，所述的初始条件包括：磁通密度B，物质密度ρ，

声速媒介声速分布C_m×m和采集的声压信号Q(lΔt)；

b)离散化不均匀媒介中的声压和振动速度的表达式

$\begin{array}{c}{p}_0^{l+1}(i,j)=p_0^l(i,j)-\mathrm{\ρ}(i,j)c^2(i,j)\·\mathrm{\Δt}\·\frac{1}{\mathrm{\Δs}}\·\left\{\right[v_x^{l+1/2}(i+\frac{1}{2},j)-v_x^{l+1/2}(i-\frac{1}{2},j)]\\ +[v_y^{l+1/2}(i,j+\frac{1}{2})-v_y^{l+1/2}(i,j-\frac{1}{2})]\}+c^2(i,j)\·q^l\·\mathrm{\Δt}\end{array}$

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

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

c)确定初始时刻的声压

将过程a)中确定的初始条件代入过程b)中所给出的声压方程中，求解得到初始时刻t_a＝0的声压p₀；

d)将过程c)所求t_a时刻的声压p₀代入过程b)中所给出的振动速度方程中，求解时刻的振速v_x,v_y；

e)依据过程d)所求t_b时刻的振速v_x,v_y，将换能器检测的反演声源Q代入过程b)中所给出的声压方程中，求解时刻的声压p，存储数据；

f)如果t_a＝T_d，转到过程d)，否则结束；

当时t_a＝T_d时，依据上述步骤所求解得到的声压分布p(r,T_d)，即包含原始声源分布。

(5)提取样本所在区域声源分布

依据自由空间大小，将原始声源分布p₀(r,T_d)，以自由空间中心为基准点，半径为r之外的所有区域声压置零，所得到的声压值P₀即为初始声源分布。

Claims (7)

1.一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，是在信号采集过程中使用线阵换能器阵列采集相互对称的声压数据，利用对称位置换能器采集数据之间的相关性求解声穿越时间，结合代数迭代重建算法重建样本声速分布，最后利用所求声速分布求解样本初始声源分布，具体包括如下阶段：

1)在扫描过程中采用平面线阵换能器阵列采集样本内部发出的磁声信号；

2)求解换能器之间的声穿越时间，包括：

(1)设定换能器检测的声压满足的条件；(2)进行声压补偿；(3)求解对称位置换能器所检测声压的相关系数；(4)求解两个轴线相互重叠的换能器聚焦区轴线部位的声穿越时间；(5)求解得到整个扫描过程的声穿越时间；

3)求解被测样本的声速分布图，包括：

(1)建立对称位置换能器连线的声速平均矩阵；(2)修正平均声速矩阵；(3)划分各单元网格；(4)求解声穿越时间；(5)对结果再次修正之后得到每个网格实际声速；

4)求解原始声源分布，包括：

(1)声速重排及声速空间扩展；(2)反转声源求解及声压转置；(3)声压数据重采样；(4)求解原始声源；(5)提取样本所在区域声源分布。

2.根据权利要求1所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，阶段1)在扫描过程中每一个换能器都要满足如下条件：

(1)换能器扫描过中完成一次扫描后进行旋转，旋转次数要满足2M-1次；

(2)每一个角度的聚焦区域移动次数为N，移动间距Δd保持不变，且满足N×Δd≥R，其中R为能容纳样本的最小正方形区域的边长，M和N是大于等于1的整数；

(3)扫描过程中，换能器旋转半径r要满足：

$r\≥\sqrt{2}R;$

(4)处于对称位置的换能器，在聚焦区轴线重叠的区域内每一个点满足如下条件：

L₁+L₂≈2r

其中L₁为所述点到左边换能器的距离，L₂为所述点到右边换能器的距离，r为换能器扫描过程中的旋转半径。

3.根据权利要求1所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，阶段2)所述的求解换能器之间的声穿越时间，包括：

(1)设定换能器检测的声压满足如下条件：

$\tilde{p}(r_0,t)=-\frac{1}{4\mathrm{\π}}{\∫\∫}_{\mathrm{\Ω}}\frac{\▿\·(J\left(r^{`}\right)\×B\left(r^{`}\right))}{|r_0-r^{`}|}\mathrm{d\Ω}{|}_{T(r_0,r^{`})=t}$

其中r₀为换能器聚焦区轴线的延长线与换能器的交点，t为换能器接收的时间轴，为换能器所接收的声压，J为感应电流密度，B为稳恒磁场强度，r`为样本上某一点的位置，T(r₀,-r₀)为重叠的聚焦区轴线间的声传播时间。

(2)进行声压补偿

定义换能器接收信号的修正值S(r₀,t)满足：

$S(r_0,t)=\tilde{p}(r_0,t)\·t$

(3)求解对称位置换能器所检测的声压的相关系数

相对于位置r₀的换能器接收的声压修正后为S(r₀,t)，设定与r₀位置对称的信号采集点位置为-r₀，故所述信号采集点-r₀接收的声压修正后为S(-r₀,t)，考虑声的传播特性，位置r₀和位置-r₀的声压信号修正值在时间和幅度上均具有高度相关，定义满足：

$R_{r_0}\left(t\right)={\∫}_{-\∞}^{\∞}S(r_0,\mathrm{\τ})S(-r_0,-(\mathrm{\τ}-t))\mathrm{d\τ}$

为聚焦区域轴线相互重叠的一组换能器所接收信号的相关系数；

(4)求解两个换能器聚焦区轴线相互重叠部位的声穿越时间

由相关性原理及卷积定理可以得到，当t＝T(r₀,-r₀)时，式T(r₀,-r₀)可以得到最大值，故由下式得到两个换能器聚焦区轴线相互重叠部位的声穿越时间

$T(r_0,-r_0)=\arg \max \left(R_{r_0}\left(t\right)\right);$

(5)求解得到整个扫描过程的声穿越时间

对阶段1)扫描所得到的M×N对重叠的聚焦区轴线位置，依次重复步骤(1)～步骤(4)，得到M×N组穿越时间。

4.根据权利要求1所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，阶段3)所述的求解被测样本的声速分布图，具体包括：

(1)建立对称位置换能器连线的传播时间矩阵

对于扫描过程中换能器动态聚焦区的位置分布，定义传播时间矩阵T_NM的任何一个元素要满足如下公式：

T_i＝T(r_i,-r_i)；1≤i≤N×M

其中r_i,-r_i代表第i对聚焦区轴线延长线与换能器交点所处的位置；

(2)修正传播时间矩阵T_NM，得到声穿越时间T'_MN

设定耦合剂的所处区域作为背景，背景声速为c₀，对于扫描过程中换能器的动态聚焦区位置分布，得到对称位置换能器之间的距离等同于旋转直径，故排除背景声速c₀影响后的声穿越时间T'_MN的元素满足：

$T_{\mathrm{MN}}^{\′}=T_{\mathrm{MN}}-\frac{2r}{c_0}\·A_{\mathrm{MN}}$

A_MN为元素全为1的M×N矩阵；

由上式得到声穿越时间T'_MN，T'_MN是被测样本被换能器聚焦区轴线穿过部位的声穿越时间再减去相同距离耦合剂的声穿越时间；

(3)划分各单元网格

将一个感兴趣的包含被测样本的扫描空间的区域分成n×n的网格，网格大小为Δs，每个网格的声穿越时间为t_j；1≤j≤n²，网格大小Δs应当满足如下：

$\mathrm{\Δs}\≤\frac{c_{\min }}{N{f}_{\max }};N\≥10$

上式中f_max为磁声耦合成像自由空间中的最大声波频率，满足：

f_max≥2f_in

f_in为磁声耦合成像的激励信号频率，C_min为磁声耦合系统中的最小声速；

(4)求解声穿越时间

根据步骤(3)的网格，按如下公式进行迭代重建，求解每一个网格的声穿越时间值t_j：

$t_j^{(k+1)}=t_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}\frac{T_i^{\′}-\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}{t}_m}{\underset{m=1}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{im}}^2};1\≤i\≤N\×M,1\≤j\≤n^2$

k为迭代次数，λ为松弛因子，a为权重因子，当网格m对声穿越时间t_j有贡献时取值为1，否则为0；

(5)对结果再次修正之后得到每个网格实际声速

步骤(4)求解得到的声穿越时间为去除背景声速之后的声穿越时间，对结果再次修正之后得到每个网格实际声速c_j；1≤j≤n²满足：

$c_j=\frac{\mathrm{\Δs}}{t_j}+c_0;1\≤j\≤n^2.$

5.根据权利要求4所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，步骤(4)所述的迭代重建步骤包括：

a)给第j个网格的声穿越时间t_j赋初值：

$t_j^{\left(k\right)}=t_j^{\left(0\right)},1\≤j\≤n^2$

式中，初值选取依据人体组织的平均声速和耦合剂的背景声速，满足：

$t_j^{\left(0\right)}=\frac{\mathrm{\Δs}}{{\stackrel{\‾}{c}}_p}-\frac{\mathrm{\Δs}}{c_0}$

其中为人体组织的平均声速；

b)计算第i对换能器间的声穿越时间修正后的估计值：

$\widetilde{T_i^{\′}}=\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{x}_j^{\left(k\right)}$

c)计算误差：

$\mathrm{\Δi}=T_i^{\′}-\widetilde{T_i^{\′}}$

d)计算第j个未知量的修正值：

$C_{\mathrm{ij}}^V=\mathrm{\Δi}\frac{a_{\mathrm{ij}}}{\underset{j}{\overset{n^2}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}^2}$

e)对t_j的值进行修正：

$T_j^{(k+1)}=T_j^{\left(k\right)}+{\mathrm{\λ}}^{\left(k\right)}{C}_{\mathrm{ij}}^V$

f)k＝k+1；重复步骤b)到f)，直到完成所有声穿越时间方程，即完成一轮迭代；

g)以上一轮迭代作为初值，重复b)到f)进行新一轮迭代，直到得到符合收敛要求的结果。

6.根据权利要求1所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，阶段4)所述的求解原始声源分布，具体包括：

(1)声速重排及声速空间扩展

利用阶段3)所求得到的每个网格实际声速c_j；1≤j≤n²，进行声速重排，得到n×n的声速分布矩阵将所述声速分布矩阵向外扩充为一个m×m；m>n的自由空间，声速分布为C_m×m；扩充部分的网格的声速为c₀；

(2)反转声源求解及声压转置

对扫描过程中换能器接收到的声压信号为p(r,t)，其中r表示空间位置；依据时间反演理论，将采集的声压信号反转作为声源，定义用于反演的声源项为q(r,t)满足：

q(r,t)＝p(r,T_d-t)；1≤i≤N×M

T_d表示检测过程中，空间每一个点的检测时间，故得基于时间反演的波动方程如下：

${\▿}^2{p}_0(r,t)-\frac{1}{c{\left(r\right)}^2}\frac{{\∂}^2{p}_0(r,t)}{\∂t^2}=q(r,t)$

c(r)为声速分布C_m×m圆柱坐标形式，p₀(r,T_d)为需要求解的原始声源分布；

(3)声压数据重采样

定义时间间隔Δt： $\mathrm{\Δt}\≤\frac{1}{c_{\max }\sqrt{3*{\left(\frac{1}{\mathrm{\Δs}}\right)}^2}}$

其中，c_max为磁声耦合成像实验中的最大声速，

依据时间间隔Δt对反演的声源项q(r,t)进行重采样，采样结果用反演声源Q(t)表示：

$Q\left(t\right)=q^l=q\left(\mathrm{l\Δt}\right);1\≤l\≤\frac{T_d}{\mathrm{\Δt}};$

(4)时间和空间交替迭代求解原始声源

依据声波波动方程，结合磁声耦合成像声源产生机制，得到不均匀媒介中的声压和振动速度的表达式，如下：

$\frac{\∂p_0}{\∂t}=-\mathrm{\ρ}{C_{m\×m}}^2(\frac{\∂v_x}{\∂x}+\frac{\∂v_y}{\∂y})+q$

$\left\{\begin{array}{c}\frac{\∂v_x}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂x}\\ \frac{\∂v_y}{\∂t}=-\frac{1}{\mathrm{\ρ}}\frac{\∂p_0}{\∂y}\end{array}\right.$

其中，v_x,v_y表示质点二维平面x,y方向的振动速度，p₀表示声压，ρ表示物质密度，c表示声速，t表示波动持续时间。采用时域有限差分的方法，在时间和空间上交替离散上述表达式，利用反演声源Q求解下一个时刻的质点振动速度，并利用质点振动速度求解声压，依此交替；

(5)提取样本所在区域声源分布

依据自由空间大小，将原始声源分布p₀(r,T_d)，以自由空间中心为基准点，半径为r之外的所有区域声压置零，所得到的声压值P₀即为初始声源分布。

7.根据权利要求6所述的一种应用于声学非均匀媒介磁声耦合成像重建方法，其特征在于，步骤(4)所述的不均匀媒介中的声压和振动速度的表达式具体求解过程如下：

a)确定所述表达式求解的初始条件，所述的初始条件包括：磁通密度B，物质密度ρ，

声速媒介声速分布C_m×m和采集的声压信号Q(lΔt)；

b)离散化不均匀媒介中的声压和振动速度的表达式

$\begin{array}{c}{p}_0^{i+1}(i,j)=p_0^l(i,j)-\mathrm{\ρ}(i,j)c^2(i,j)\·\mathrm{\Δt}\·\frac{1}{\mathrm{\Δs}}\·\left\{\right[v_x^{l+1/2}(i+\frac{1}{2},j)-v_x^{l+1/2}(i-\frac{1}{2},j)]\\ +[v_y^{l+1/2}(i,j+\frac{1}{2})-v_y^{l+1/2}(i,j-\frac{1}{2})]\}+c^2(i,j)\·q^l\·\mathrm{\Δt}\end{array}$

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

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

c)确定初始时刻的声压

将过程a)中确定的初始条件代入过程b)中所给出的声压方程中，求解得到初始时刻t_a＝0的声压p₀；

d)将过程c)所求t_a时刻的声压p₀代入过程b)中所给出的振动速度方程中，求解时刻的振速v_x,v_y；

e)依据过程d)所求t_b时刻的振速v_x,v_y，将换能器检测的反演声源Q代入过程b)中所给出的声压方程中，求解时刻的声压p，存储数据；

f)如果t_a＝T_d，转到过程d)，否则结束；

当时t_a＝T_d时，依据上述步骤所求解得到的声压分布p(r,T_d)，即包含原始声源分布。