一种合成孔径数据的高信噪比成像方法

1.本发明属于检测方法技术领域，具体涉及一种合成孔径数据的高信噪比成像方法。


背景技术：

2.超声在临床诊断中是一种重要的成像工具，因为它能够提供实时的解剖和功能图像，同时不会引起任何副作用。随着电子信息技术的发展，合成孔径(synthetic transmit aperture，sta)数据已成为研究热点，它包含相控阵中每一对发射接收对的响应，该数据集已被广泛用于各种场景，如组织结构成像、流量估计和组织运动补偿等。sta数据集通常是通过单个振元依次激发和全部振元接收得到的，它采用双向动态聚焦来进行成像，使得横向空间分辨率和对比度大大提高，因为它采用了真正的波传播路径。此外，sta数据集可以从其他波束形成策略中恢复，一般它们采用多元素按照一定的权重同时激发，这样信噪比(snr)和帧率可以得到改善。
3.sta数据中包含了来自相控阵探头的关于测量区域的最大可能的独立信息，需要适当的成像算法对它进行处理以重建反射体的位置。延迟求和法(das)是医学成像中的一种标准方法，通过将波束作为射线来估计传播时间，并根据时间延迟和通过最小方差优化得到的权重来叠加数据以重建目标。这种方法在鲁棒性和信噪比方面有一定的优势，但当测量对象为多层结构时，需要耗时的路径计算。
4.此外，机器学习的快速发展也为超声成像提供了新的选择。然而，大多数研究人员将机器学习作为对das结果进行预处理或后处理的手段，用于超声数据集图像重建的神经网络(nn)非常有限。在此，hyun采用nn框架处理sta数据得到成像结果，它可以抑制斑点伪像；youn也初步尝试通过卷积神经网络直接处理超声信号来检测散射体。由于sta数据集包含大量的信号，需要一个庞大的网络来获得从输入到输出的正确和具有一定普适性的映射，这取决于庞大而多样的训练数据集，所以巨大的工作量目前阻碍了这项技术的推广。
5.全矩阵捕获(fmc)是从sta方法演变而来的，并已被广泛用于工业检测中。由于fmc和sta的数据集相同，fmc的新成像方法可以转化为sta，包括相移迁移(psm)方法，逆时偏移(rtm)方法，以及基于模型的方法。psm方法将接收到的信号视为测量区域表面的波场，用傅里叶变换分离各个谐波分量，通过相移算子重建测量区域的波场，之后实施成像条件来获得目标的聚焦图像。在测量多层结构的缺陷时，psm表现出很大的效率优势。rtm方法使用全波方程来重建测量区域的波场，然后执行成像条件来产生聚焦图像。频域rtm能够同时解决不同的激励，这使得计算复杂性的显著降低。这种方法不受声速模型复杂性的限制，可以获得高分辨率的聚焦图像，但它对声速模型的准确程度非常敏感。基于模型的方法使用矩阵方程来建立fmc数据集和测量目标之间的关系，并采用反演算法可以恢复高分辨的图像。尽管如此，由于fmc或sta数据集的规模，矩阵方程特别大，使该方法的计算成本很高，难以应用。
6.上述成像方法均假定测量区域的声速是均匀的。然而，生物体内的声速通常是轻
微不均匀的，其中亚分辨率散射的随机干扰会导致成像结果中出现斑点伪影。斑点通常被认为是恶化成像结果的噪声，因为它降低了目标的分辨率和信噪比。因此，为了获得高信噪比的目标图像，需要在很大程度上抑制斑点噪声，然而现有技术中还没有能够很好地抑制斑点噪声的成像方法，得到的目标图像信噪比较低。


技术实现要素：

7.为解决现有技术中存在的问题，本发明提供一种合成孔径数据的高信噪比成像方法，该方法既能消除sta数据中的大部分噪声，还能抑制成像过程中的斑点伪影，有效提高成像结果的信噪比和分辨率。
8.一种合成孔径数据(sta)的高信噪比成像方法，包括以下步骤：
9.(1)输入合成孔径数据，将合成孔径数据偏移到目标层的上表面，得到目标层数据；
10.(2)对目标层数据进行奇异值分解(svd)，提取出信号子空间，并计算每个频率的交叉谱矩阵(cms)；
11.(3)对每个频率的交叉谱矩阵进行反卷积，得到每个频率对应的成像结果，将各频率对应的成像结果进行叠加，得到最终的成像结果。
12.上述步骤中，通过对特征值的分析来划分信号子空间和噪声子空间，其中，任意划分方法均可使用。对于均匀介质，则实际操作时无需将合成孔径数据进行偏移，直接对测量得到的合成孔径数据进行后续处理即可。
13.作为优选，步骤(2)中，对目标层数据进行奇异值分解的形式为：
[0014][0015]
其中，d(ω)为频率为ω时的目标层数据；下标sig和n分别代表信号子空间和噪声子空间；ω为频率；u
sig
(ω)、v
sig
(ω)分别表示信号子空间的两组特征向量；∑
sig
(ω)为信号子空间的特征值；un(ω)、vn(ω)分别为噪声子空间的两组特征向量；∑n(ω)为噪声子空间的特征值；上标h为转置和复数共轭。
[0016]
合成发射孔径(sta)利用相控阵中的每个振元依次将球面波发射到目标区域，同时相控阵中的所有元件接收来自介质的反射波。对于n个振元的线阵，这种采集方式下得到的采集信号包含三个维度，这里表示为其中，n
t
为时域中的采样点数。因此，sta数据是一个响应矩阵，它包含了所有振元之间的响应关系。如图2中(a)所示，波从位于rs＝(xs,0)处的振元发射，在介质中传播遇到反射体后发生反射，反射波被位于rr＝(xr,0)处的振元接收。
[0017]
声波在发射振元rs和接收振元rr之间的传播可以在频域中建模得到：
[0018][0019]
其中，m为反射体的个数，反射体m位于rm＝(xm,zm)处，且反射系数为εm；f(ω)为激发信号的频谱；ds(θs,ω)和dr(θr,ω)分别为发射振元和接收振元的方向系数；θs和θr分别
为发射振元和接收振元的方位角；ω为频率。
[0020]
其中，方向系数定义为：
[0021][0022]
其中，a为振元宽度，c为声速。
[0023]
g是格林方程，由下式得到：
[0024][0025]
其中，h0为零阶第一类汉克尔方程；i为虚数单位；格林方程g(rm,rr,ω)为在rm点处对rr处点源的响应。
[0026]
在频域中所有振元的n
×
n响应矩阵可以写为：
[0027][0028]
其中：
[0029]gm
＝[d1(θ1,ω)g(r1,rm,ω),d2(θ2,ω)g(r2,rm,ω),
…
,dn(θ1,ω)g(rn,rm,ω)]
t
[0030]
(5)
[0031]
其中，h为转置和复数共轭；t为转置。
[0032]
g为一个n
×
m矩阵，表示为：
[0033]
g＝[g1,g2,
…
,gm]
ꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0034]
m为反射体的个数；n为振元个数。
[0035]
因为成像区域内对应的不同方位的方向系数相差不大，因此可以忽略振元的方向系数，则g为一个正交矩阵。
[0036]
q为一个m
×
m维矩阵，写作：
[0037][0038]
因此，测量得到的响应矩阵d(ω)可以建模为：
[0039]
d(ω)＝[d(xs,xr,ω)]＝gqgh n(ω)
ꢀꢀꢀ
(8)
[0040]
其中，d(ω)为一个n
×
n矩阵，它是在频率为ω时d(xs,xr,ω)中所有接收-发射振元对的集合，d(xs,xr,ω)可以由d(xs,xr,t)在时间维度进行一维傅里叶变换得到。
[0041]
这里假设反射点数目小于振元数目，即m《n。因此，p(ω)为一个秩为m的理想响应矩阵，但是该理想响应矩阵受到了噪声的污染后变为d(ω)，后者是满秩的。在理论上，d(ω)的前m个特征值要远大于其余的n-m个。
[0042]
测量得到的响应矩阵可以写为以下特征值分解形式：
[0043][0044]
其中，下标sig和n分别代表信号子空间和噪声子空间，根据特征值的大小来区分信号子空间和噪声子空间；∑
sig
(ω)、∑n(ω)分别为信号子空间和噪声子空间的特征值，
二者均为m
×
m维对角矩阵，对应于m个反射点；u
sig
(ω)和v
sig
(ω)分别为信号子空间的两组特征向量，为n
×
m维矩阵，二者在理论上应该是一样的，然而实际测量中的噪声会使它们产生差异，此处忽略v
sig
(ω)，将u
sig
(ω)作为信号子空间进行后续交叉谱矩阵的计算；un(ω)、vn(ω)分别为噪声子空间的两个特征向量，同样的为n
×
m维矩阵。
[0045]
因此，svd能够将散射信号从噪声中分离出来，可以得到下面的公式：
[0046][0047]
对比公式(10)中的三项，v
sig
(ω)、∑
sig
(ω)和u
sig
(ω)分别可以解释为入射波从激发阵元发射的信号，反射体的反射系数以及反射波的接收信号。
[0048]
因此，u
sig
(ω)＝[u
sig1
,
…
,u
sigm
]可被视为传感器接收到的从散射点主动发出的波的信号，如图2(b)所示。其中，u
sigm
为信号子空间u
sig
(ω)中的信号向量，m∈[1，m]，m为反射体的个数。
[0049]
因此，通过svd来提取u
sig
(ω)不仅滤除了大多数的噪声，而且将被动的异质体检测问题转变为主动的声源定位问题，进而可以利用交叉谱矩阵(cms)进行成像。
[0050]
作为优选，步骤(2)中，交叉谱矩阵cm(ω)的计算公式为：
[0051][0052]
式中，m为反射点的个数；ω为频率；u
sigm
为信号子空间中的信号向量；上标h为转置和复数共轭。
[0053]
波束成形是主动声源定位的重要方法。声波在二维均匀介质中的波动方程可以写为：
[0054][0055]
其中，p(x,z,t)为点(x,z)处的声压场，q(t)为位于(x0,z0)的点源，这时时域的波场为：
[0056][0057]
其中，距离为时间延迟为时间延迟为
[0058]
因此，相控阵振元收到的信号为表面的波场p(xr,0,t)，声源分布可以通过反转时间延迟重建。
[0059]
在(x,z)点处可能的声源可以通过下式估计：
[0060][0061]
其中，(xj,0)为第n个振元的位置；tj＝rj/c。
[0062]
对进行傅里叶变换后，得到波束成形在频域的表达式如下：
[0063][0064]
其中，h(x,z,ω)和y(x,z,ω)为n
×
1向量，分别写作：
[0065]
h(x,z,ω)＝[2πr1exp(iωt1),
…
,2πrnexp(iωtn)]
t
ꢀꢀꢀ
(15)
[0066]
y(x,z,ω)＝[p(x1,0,ω),
…
,p(xn,0,ω)]
t
ꢀꢀꢀ
(16)
[0067]
是复数，它的能量可以由下式计算：
[0068][0069]
其中，上标*表示共轭；c(ω)为包含测量信号的n
×
n矩阵，c(ω)为交叉谱矩阵(cms)，其中任一坐标为(i,j)的元素c
ij
也可以由下式得到：
[0070]cij
＝p(xi,0,ω)p
*
(xj,0,ω)
ꢀꢀꢀꢀ
(18)
[0071]
对于sta数据，u
sig
(ω)通过svd方法提取得到，并被当作是在时间t＝0时从异质体(反射体)主动发出的波经过相控阵接收得到的信号，这里定义成像条件如下：
[0072]
i(x,z)＝∫i(x,z,ω)dω＝∫hh(x,z,ω)cm(ω)h(x,z,ω)dω
ꢀꢀꢀ
(19)
[0073]
其中：
[0074][0075]
其中，u
sigm
为信号子空间u
sig
(ω)中的信号向量。
[0076]
根据提取的u
sig
(ω)先计算出交叉普矩阵cm(ω)，再利用成像条件以及cm(ω)进行成像，得到频率为ω时对应的成像结果。
[0077]
作为优选，步骤(3)中，利用非稳态点扩散函数(psf)对得到的每个频率的交叉谱矩阵进行反卷积。
[0078]
相控阵对点源的成像称为点扩散函数(psf)，它与硬件参数有关，如元件尺寸、元件间距和元件接收性能。相控阵的成像结果是psf和真实图像的卷积结果。
[0079]
针对一个点源λ＝(x
λ
,z
λ
)，声源引入一个csm为：
[0080]cλ
(ω)＝y(x
λ
,z
λ
,ω)yh(x
λ
,z
λ
,ω)
ꢀꢀꢀ
(21)
[0081]
因此，点扩散函数psf可以写作：
[0082]iλ
(x,z,ω)＝hh(x,z,ω)c
λ
(ω)h(x,z,ω)
ꢀꢀꢀ
(22)
[0083]
其中，psf是非稳态的，因为它是光源位置的函数。psf模糊了图像，使相邻的散射点无法区分，因此有必要采用反卷积来消除psf的影响并提高分辨率。
[0084]
作为优选，步骤(3)中，对每个频率的交叉谱矩阵进行反卷积时，当满足如下条件时，终止迭代：
[0085][0086]
其中，为第j 1次迭代得到的交叉谱矩阵所包含的信息。
[0087]
稀疏反卷积方法可以有效抑制斑点伪影，提高信噪比。clean算法是一种二维正交匹配追踪(omp)，其功能类似于l1范数的稀疏去卷积。clean算法假设源由散射点组成，它在
每次迭代时找到振幅最大的源，并用psf反卷积去除该点对成像结果的贡献。
[0088]
如图3所示，由真实图像(由交叉谱矩阵成像得到)和psf之间的卷积产生的波束成形图像可以被看作是相应源位置上psf的叠加。波束成形图像被定义为一个"脏"图像，其中三个源几乎是不可分离的。然后将最大振幅点提取到“净”图像中，并从“脏”图像中减去该点的相应psf。上述步骤重复进行，直到收敛重建出真实图像(即“净”图像)。
[0089]
在实施时，通过对cms进行更新来去除psf的作用：
[0090][0091]
其中，是在第j次迭代时的cms；γ为控制参数，它决定了去除的功率的大小；为在第j张“脏”图像(第j次迭代对应的脏图)中的像素强度的最大值，和分别为峰值点(最大振幅点)的水平和竖直坐标。
[0092]
然后，可以得到第j 1次迭代时的“脏”图像i
j 1
(x,z,ω)：
[0093][0094]
进一步，“净”图像i
clean
收集到峰值点信息：
[0095][0096]
当前迭代得到的csm包含信息比前一次迭代得到的csm包含信息多时，迭代终止，该终止迭代条件写作下式：
[0097][0098]
式中，为第j次迭代得到的csm包含的信息。
[0099]
作为优选，步骤(3)中，将各频率对应的成像结果进行叠加的计算公式如下：
[0100]
i(x,z)＝∫i
clean
(x,z,ω)dω
[0101]
其中，i(x,z)为最终的成像结果；i
clean
(x,z,ω)为每个频率对应的成像结果；ω为频率；dω为频域的离散步长。
[0102]
上述成像方法采用逐频迭代的方式，计算出选择频带范围[ω
min
,ω
max
]内各个离散频率处对应的声压场，其中频带范围由全矩阵数据决定，频域离散步长dω由采样频率决定。
[0103]
作为优选，步骤(1)中，先对合成孔径数据滤除直达波，再进行偏移。滤除直达波的方法可以是任意滤波方法，这里给出一种可行方法，可以采用窗函数来选择有效的信号区段。
[0104]
作为优选，步骤(1)中，将合成孔径数据偏移到目标层的上表面时，先将合成孔径数据转换到频率波数域，进行外推后，进行傅里叶逆变换，得到时间空间域的目标层数据；
[0105]
在频率波数域，按照如下公式对合成孔径数据进行偏移：
[0106][0107]
其中，为目标层数据；为频率波
数域的合成孔径数据；z
tl-1
为目标层上表面深度；ω为频率；为激发振元在水平方向上的波数；为接收振元在水平方向上的波数；dj为第j层的厚度；为第j层在竖直方向上的波数；i为虚数单位。
[0108]
非均匀介质中的层状结构会导致严重的声折射和混响，并使得图像质量恶化。因此，开发了sta数据的波场外推方法，以提高成像方法在层状结构介质中的成像质量。
[0109]
测量数据d(xs,xr,t)可以看作表面的波场，激发振元和接收振元分别都在空间中占有2个维度，所以sta数据的波场有五个维度，表示为p(xs,zs,xr,zr,t)，其中，xs和zs分别为激发振元的水平和竖直坐标，xr和zr分别为接收振元的水平和竖直坐标。
[0110]
表面的波场表示为：
[0111]
p(xs,0,xr,0,t)＝d(xs,xr,t)
ꢀꢀꢀꢀꢀꢀꢀ
(27)
[0112]
将表面波场通过一个三维傅里叶变换得到在频率波数域的波场为：
[0113][0114]
其中，和分别为激发振元和接收振元的水平波数；t为时间。
[0115]
表面波场的向下偏移过程如图4(a)和(b)所示，成像目标位于第二层介质中，为了进行波数成形，必须得到在第二层的上界面z1处的波场。
[0116]
首先，接收波场按照下式偏移：
[0117][0118]
其中，d1为第一层的厚度；为接收振元的竖直波数，定义为：
[0119][0120]
式中，c1为第一层的声速；ω为频率。
[0121]
偏移后，接收振元和发射振元的界面图如图4(c)和(d)所示，其中接收振元被偏移到界面z1处，而激发阵元仍然在测量表面上。
[0122]
进一步，将激发振元按照下式偏移到界面z1处：
[0123][0124]
其中，为激发振元竖直方向波数，定义为：
[0125][0126]
这时，如图4(e)和(f)中的横截面图所示，激发和接收振元均位于接口处。
[0127]
总的来说，表面上的波场可以直接偏移到深度z1：
[0128][0129]
其中，kz为向下偏移的总波数，定义为：
[0130][0131]
然后，将深度z1处的波场转换回空间域，得到了等价的在此深度处测量得到的sta数据：
[0132][0133]
可以用于对第二层中的异质体成像。
[0134]
与现有技术相比，本发明的有益效果为：
[0135]
本发明的合成孔径数据的高信噪比成像方法，利用奇异值分解(svd)来滤除sta数据中的大部分噪声提取出信号子空间，并将被动异质性检测问题转化为主动源定位问题。根据提取的信号子空间计算交叉谱矩阵并用于波束成形，同时利用非稳态点扩散函数(psf)对交叉谱矩阵进行反卷积，进一步抑制斑点伪影，大大改善了信噪比，提高了成像质量。另外还引入了相位偏移，使该成像方法能够满足分层介质的成像需求，适用性广。
附图说明
[0136]
图1为本发明实施例的成像方法流程图；
[0137]
图2中：(a)为合成孔径数据采集时的波传播示意图；(b)为信号子空间等价为主动声源定位问题的波传播是示意图；
[0138]
图3为针对波束成形图像进行psf反卷积来重建真实图像；
[0139]
图4中：(a)为sta数据对应同一激发振元的截面示意图；(b)为sta数据对应同一接收振元的截面示意图；(c)为向下偏移接收振元到深度z1后sta数据对应同一激发振元的截面示意图；(d)为向下偏移接收振元到深度z1后sta数据对应同一接收振元的截面示意图；(e)为向下偏移接收和发射振元到深度z1后sta数据对应同一激发振元的截面示意图；(f)为向下偏移接收和发射振元到深度z1后sta数据对应同一接收振元的截面示意图；
[0140]
图5中：(a)为猪肉直接测量示意图；(b)为(a)中直接测量得到的sta数据；(c)为猪肉楔块辅助测量示意图；(d)为(b)中楔块辅助测量得到的sta数据；(a)中直接测量数据由不同方法进行图像重建的结果：(e)为das、(f)为psm、(g)为rtm、(h)为本发明实施例中的方法；(b)中楔块辅助测量数据由不同方法进行图像重建的结果：(i)为das、(j)为psm、(k)为rtm、(l)本发明实施例中的方法。
具体实施方式
[0141]
如图1所示，一种合成孔径数据(sta)的高信噪比成像方法，包括以下步骤：
[0142]
(1)预处理：输入合成孔径数据d(xs,xr,t)，并采用窗函数来选择有效的信号区段，滤除直达波，得到
[0143]
(2)将预处理后的合成孔径数据先进行三维傅里叶变换到频率波数域，得到将通过下式进行偏移，得到目标层(tl层)在频率波数域的数据推导过程如下：
[0144]
非均匀介质中的层状结构会导致严重的声折射和混响，并使得图像质量恶化。因此，开发了sta数据的波场外推方法，以提高成像方法在层状结构介质中的成像质量。
[0145]
测量数据d(xs,xr,t)可以看作表面的波场，激发振元和接收振元分别都在空间中占有2个维度，所以sta数据的波场有五个维度，表示为p(xs,zs,xr,zr,t)，其中，xs和zs分别为激发振元的水平和竖直坐标，xr和zr分别为接收振元的水平和竖直坐标。
[0146]
表面的波场表示为：
[0147]
p(xs,0,xr,0,t)＝d(xs,xr,t)
ꢀꢀꢀꢀꢀ
(27)
[0148]
将表面波场通过一个三维傅里叶变换得到在频率波数域的波场为：
[0149][0150]
其中，和分别为激发振元和接收振元的水平波数；t为时间。
[0151]
表面波场的向下偏移过程如图4(a)和(b)所示，设定成像目标位于第二层介质中，为了进行波数成形，必须得到在第二层的上界面z1处的波场。
[0152]
首先，接收波场按照下式偏移：
[0153][0154]
其中，d1为第1层的厚度；为接收振元的竖直波数，定义为：
[0155][0156]
式中，c1为第一层的声速。
[0157]
偏移后，接收振元和发射振元的界面图如图4(c)和(d)所示，其中接收振元被偏移到界面z1处，而激发阵元仍然在测量表面上。
[0158]
进一步，将激发振元按照下式偏移到界面z1处：
[0159][0160]
其中，为激发振元竖直方向波数，定义为：
[0161][0162]
这时，如图4(e)和(f)中的横截面图所示，激发和接收振元均位于接口处。
[0163]
总的来说，表面上的波场可以直接偏移到深度z1：
[0164][0165]
其中，kz为向下偏移的总波数，定义为：
[0166][0167]
然后，将深度z1处的波场转换回空间域，得到了等价的在此深度处测量得到的sta数据：
[0168]
[0169]
可以用于对第二层中的异质体成像。
[0170]
由上述推导可得，在频率波数域的外推公式为：
[0171][0172]
式中，和分别为激发振元和接收振元的水平波数；z
tl-1
为目标层上表面深度；ω为频率；dj为第j层的厚度；为第j层在竖直方向上的波数；i为虚数单位。
[0173]
后将进行傅里叶逆变换得到目标层的合成孔径数据(目标层数据)
[0174][0175]
(3)对进行奇异值分解(svd)，提取出信号子空间。
[0176]
合成发射孔径(sta)利用相控阵中的每个振元依次将球面波发射到目标区域，同时相控阵中的所有元件接收来自介质的反射波。对于n个振元的线阵，这种采集方式下得到的采集信号包含三个维度，这里表示为其中，n
t
为时域中的采样点数。因此，sta数据是一个响应矩阵，它包含了所有振元之间的响应关系。如图2中(a)所示，波从位于rs＝(xs,0)处的振元发射，在介质中传播遇到反射体后发生反射，反射波被位于rr＝(xr,0)处的振元接收。
[0177]
声波在发射振元rs和接收振元rr之间的传播可以在频域中建模得到：
[0178][0179]
其中，m为反射体的个数，反射体位于rm＝(xm,zm)处，且反射系数为εm；f(ω)为激发信号的频谱；ds(θs,ω)和d
t
(θr,ω)分别为发射振元和接收振元的方向系数；θs和θr分别为发射振元和接收振元的方位角；ω为频率。
[0180]
其中，方向系数定义为：
[0181][0182]
其中，a为振元宽度，c为声速。
[0183]
g是格林方程，由下式得到：
[0184][0185]
其中，h0为零阶第一类汉克尔方程；i为虚数单位；格林方程g(rm,rr,ω)表示在rm点处对rr处点源的响应。
[0186]
在频域中所有振元的n
×
n响应矩阵可以写为：
[0187][0188]
其中：
[0189]gm
＝[d1(θ1,ω)g(r1,rm,ω),d2(θ2,ω)g(r2,rm,ω),
…
,dn(θ1,ω)g(rn,rm,ω)]
t
[0190]
(5)
[0191]
其中，h为转置和复数共轭；t为转置。
[0192]
g为一个n
×
m矩阵，表示为：
[0193]
g＝[g1,g2,
…
,gm]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0194]
m为反射体的个数；n为振元个数。
[0195]
因为成像区域对应的不同方位的方向系数相差不大，因此可以忽略振元的方向系数，则g为一个正交矩阵。
[0196]
q为一个m
×
m维矩阵，写作：
[0197][0198]
因此，测量得到的响应矩阵d(ω)可以建模为：
[0199]
d(ω)＝[d(xs,xr,ω)]＝gqgh n(ω)
ꢀꢀꢀ
(8)
[0200]
其中d(ω)为一个n
×
n矩阵，它是在频率为ω时d(xs,xr,ω)中所有接收-发射振元对的集合，d(xs,xr,ω)可以由d(xs,xr,t)在时间维度进行一维傅里叶变换得到。
[0201]
这里假设反射点数目小于振元数目，即m《n。因此，p(ω)为一个秩为m的理想响应矩阵，但是该理想响应矩阵受到了噪声的污染后变为d(ω)，后者是满秩的。在理论上，d(ω)的前m个特征值要远大于其余的n-m个。
[0202]
测量得到的响应矩阵可以写为以下特征值分解形式：
[0203][0204]
其中，下标sig和n分别代表信号子空间和噪声子空间，根据特征值的大小来区分信号子空间和噪声子空间；∑
sig
(ω)、∑n(ω)分别为信号子空间和噪声子空间的特征值，二者均为m
×
m维对角矩阵，对应于m个反射点；u
sig
(ω)和v
sig
(ω)分别为信号子空间的两组特征向量，为n
×
m维矩阵，二者在理论上应该是一样的，然而实际测量中的噪声会使它们产生差异此处忽略v
sig
(ω)，将u
sig
(ω)作为信号子空间进行后续交叉谱矩阵的计算；un(ω)、vn(ω)分别为噪声子空间的两组特征向量，同样的为n
×
m维矩阵。
[0205]
因此，svd能够将散射信号从噪声中分离出来，可以得到下面的公式：
[0206][0207]
对比公式(10)中的三项，v
sig
(ω)、∑
sig
(ω)和u
sig
(ω)分别可以解释为入射波从激发阵元发射的信号，反射体的反射系数以及反射波的接收信号。
[0208]
因此，u
sig
(ω)＝[u
sig1
,
…
,u
sigm
]可被视为传感器接收到的从散射点主动发出的波的信号，如图2(b)所示。其中，u
sigm
为信号子空间u
sig
(ω)中的信号向量，m∈[1，m]，m为反射体的个数。
[0209]
因此，通过svd来提取u
sig
(ω)不仅滤除了大多数的噪声，而且将被动的异质体检测问题转变为主动的声源定位问题。
[0210]
(4)根据提取的信号子空间，计算每个频率的交叉谱矩阵。
[0211]
波束成形是主动声源定位的重要方法。声波在二维均匀介质中的波动方程可以写为：
[0212][0213]
其中，p(x,z,t)为点(x,z)处的声压场，q(t)为位于(x0,z0)的点源，这时时域的波场为：
[0214][0215]
其中，距离为时间延迟为时间延迟为
[0216]
因此，相控阵振元收到的信号为表面的波场p(xr,0,t)，声源分布可以通过反转时间延迟重建。在(x,z)点处可能的声源可以通过下式估计：
[0217][0218]
其中，(xj,0)为第n个振元的位置；tj＝rj/c。
[0219]
对进行傅里叶变换后，得到波束成形在频域的表达式如下：
[0220][0221]
其中，h(x,z,ω)和y(x,z,ω)为n
×
1向量，分别写作：
[0222]
h(x,z,ω)＝[2πr1exp(iωt1),
…
,2πr
n exp(iωtn)]
t
ꢀꢀꢀ
(15)
[0223]
y(x,z,ω)＝[p(x1,0,ω),
…
,p(xn,0,ω)]
t
ꢀꢀꢀ
(16)
[0224]
是复数，它的能量可以由下式计算：
[0225][0226]
其中，上标*表示共轭；c(ω)为包含测量信号的n
×
n矩阵，c(ω)为交叉谱矩阵(cms)，其中任一坐标为(i,j)的元素c
ij
也可以由下式得到：
[0227]cij
＝p(xi,0,ω)p
*
(xj,0,ω)
ꢀꢀꢀꢀꢀꢀ
(18)
[0228]
对于sta数据，u
sig
(ω)通过svd方法提取得到，并被当作是在时间t＝0时从异质体(反射体)主动发出的波经过相控阵接收得到的信号，这里定义成像条件如下：
[0229]
i(x,z)＝∫i(x,z,ω)dω＝∫hh(x,z,ω)cm(ω)h(x,z,ω)dω
ꢀꢀꢀ
(19)
[0230]
其中：
[0231][0232]
其中，cm(ω)为频率为ω时的交叉普矩阵；u
sigm
为信号子空间u
sig
(ω)中的信号向量。
[0233]
根据提取的u
sig
(ω)先计算出交叉普矩阵cm(ω)，在利用成像条件以及cm(ω)进行成像，得到频率为ω时对应的成像结果。
[0234]
(5)利用非稳态点扩散函数(psf)对得到的每个频率的交叉谱矩阵进行反卷积，得到每个频率对应的高信噪比成像结果。
[0235]
相控阵对点源的成像称为点扩散函数(psf)，它与硬件参数有关，如元件尺寸、元件间距和元件接收性能。相控阵的成像结果是psf和真实图像的卷积结果。
[0236]
针对一个点源λ＝(x
λ
,z
λ
)，声源引入一个csm为：
[0237]cλ
(ω)＝y(x
λ
,z
λ
,ω)yh(x
λ
,z
λ
,ω)
[0238]
因此，点扩散函数psf可以写作：
[0239]iλ
(x,z,ω)＝hh(x,z,ω)c
λ
(ω)h(x,z,ω)
[0240]
其中，psf是非稳态的，因为它是光源位置的函数。psf模糊了图像，使相邻的散射点无法区分，因此有必要采用反卷积来消除psf的影响并提高分辨率。
[0241]
稀疏反卷积方法可以有效抑制斑点伪影，提高信噪比。clean算法是一种二维正交匹配追踪(omp)，其功能类似于l1范数的稀疏去卷积。clean算法假设源由散射点组成，它在每次迭代时找到振幅最大的源，并用psf反卷积去除该点对成像结果的贡献。
[0242]
如图3所示，由真实图像(由交叉谱矩阵成像得到)和psf之间的卷积产生的波束成形图像可以被看作是相应源位置上psf的叠加。波束成形图像被定义为一个"脏"图像，其中三个源几乎是不可分离的。然后将最大振幅点提取到“净”图像中，并从“脏”图像中减去该点的相应psf。上述步骤重复进行，直到收敛重建出真实图像(即“净”图像)。
[0243]
在实施时，通过对cms进行更新来去除psf的作用：
[0244][0245]
其中，是在第j次迭代时的cms；γ为控制参数，它决定了去除的功率的大小；为在第j张“脏”图像(第j次迭代对应的脏图)中的像素强度的最大值，和分别为峰值点(最大振幅点)的水平和竖直坐标。
[0246]
然后，可以得到第j 1次迭代时的“脏”图像i
j 1
(x,z,ω)：
[0247][0248]
进一步，“净”图像i
clean
收集到峰值点信息：
[0249][0250]
当前迭代得到的csm包含信息比前一次迭代得到的csm包含信息多时，迭代终止，该终止迭代条件写作下式：
[0251][0252]
式中，为第j次迭代得到的csm包含的信息。
[0253]
(6)将各频率对应的成像结果进行叠加，得到最终的成像结果。
[0254]
其中叠加的计算公式为：
[0255]
i(x,z)＝∫i
clean
(x,z,ω)dω
[0256]
其中，i(x,z)为最终的成像结果；i
clean
(x,z,ω)为每个频率对应的成像结果；dω为频域的离散步长。
[0257]
检测实验
[0258]
进行了体外实验以评估所提出的方法。采用猪大排肉作为基础组织，并将三根碳棒(φ＝1mm)插入该组织作为目标(反射体)。猪肉有微弱的异质性，其整体声速标定为1650m/s。首先，如图5中(a)所示，一个中心频率为2.25mhz、振元间隔为0.75mm的64元阵列换能器(汕头超声)直接贴在猪肉表面，并使用64/64oem-pa(aos.ltd,america)来收集图5中(b)中的sta数据集，其中采样频率为50mhz，时间范围为32μs。此外，一个声速为2337m/s的楔块(sc63-nl-z20)被安装在相控阵上以协助测量，如图5中(c)所示。然后收集第二个sta数据集，时间范围增加到42μs，如图5中(d)。
[0259]
为了定量分析成像结果，首先引入信噪比指数snr，它反映了图像中的测量目标区域(df)和背景区域(n-df)之间的对比度，它可以对信噪比进行定量表征，计算公式如下：
[0260][0261]
其中，和分别为目标区域和背景区域的像素幅值，mean和std分别为平均值和标准差。
[0262]
这两个sta数据集由四种成像方法处理，分别为das、psm、rtm和本实施例方法，成像结果如图5中(e)-(l)所示。在直接测量的成像结果中，只有图5中(e)中由das重建的图像不能识别最左边的目标。psm和rtm的成像结果中的背景伪影非常严重，在图5中(f)和(g)中表现为长条状和斑点状。在图5中(h)中，本实施例方法能清晰地重建出三个目标，并抑制了大部分的伪影。
[0263]
从楔块辅助测量的成像结果来看，图5中(i)中最左边的目标仍然不能被das识别。如图5中(j)和(k)所示，psm和rtm的成像结果中的斑点伪影在一定程度上被抑制，但条纹伪影仍然污染了图像。图5中(l)中，本实施例方法在抑制了所有伪影后重建了目标的清晰图像。
[0264]
表1图5中成像结果对应的信噪比和计算时间
[0265][0266]
表1列出了实验中不同成像方法的信噪比和计算时间。除rtm外，所有方法在楔块辅助测量中结果的信噪比都高于直接测量，因此楔块在这种情况下是有帮助的。本实施例方法的成像结果的信噪比远远高于其他方法。虽然psm重建的图像信噪比最低，但其计算时间比其他方法低得多，而rtm需要太多的计算时间。此外，本实施例方法比das需要的时间少得多。因此，本实施例的成像方法不仅提高了成像结果的信噪比，在效率上也非常有竞争力。