CN109187771A - 一种融合特征值分解的低复杂度最小方差超声成像方法 - Google Patents
一种融合特征值分解的低复杂度最小方差超声成像方法 Download PDFInfo
- Publication number
- CN109187771A CN109187771A CN201811243839.2A CN201811243839A CN109187771A CN 109187771 A CN109187771 A CN 109187771A CN 201811243839 A CN201811243839 A CN 201811243839A CN 109187771 A CN109187771 A CN 109187771A
- Authority
- CN
- China
- Prior art keywords
- matrix
- array
- covariance matrix
- algorithm
- beam domain
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 74
- 230000004927 fusion Effects 0.000 title claims abstract description 11
- 239000011159 matrix material Substances 0.000 claims abstract description 146
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 129
- 239000013598 vector Substances 0.000 claims abstract description 44
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 39
- 238000000034 method Methods 0.000 claims description 28
- 238000005070 sampling Methods 0.000 claims description 19
- 238000009499 grossing Methods 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 13
- 230000009467 reduction Effects 0.000 claims description 13
- 230000003044 adaptive effect Effects 0.000 claims description 10
- 238000003491 array Methods 0.000 claims description 9
- 238000012285 ultrasound imaging Methods 0.000 claims description 8
- 230000009466 transformation Effects 0.000 claims description 6
- 230000017105 transposition Effects 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 4
- 230000021615 conjugation Effects 0.000 claims description 3
- 230000003321 amplification Effects 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 claims description 2
- 238000001914 filtration Methods 0.000 claims description 2
- 238000003199 nucleic acid amplification method Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 21
- 239000000284 extract Substances 0.000 abstract 1
- 239000000523 sample Substances 0.000 description 32
- 238000010521 absorption reaction Methods 0.000 description 10
- 238000002474 experimental method Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000001514 detection method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 241000519995 Stachys sylvatica Species 0.000 description 1
- 210000001015 abdomen Anatomy 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
- 230000002792 vascular Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
Landscapes
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明涉及一种融合特征值分解的低复杂度最小方差超声成像方法,属于超声成像领域。首先,利用离散余弦变换将回波数据转换到维数较少的波束域,然后对样本协方差矩阵进行特征值分解提取信号子空间,并选取最大的特征值和其对应的特征向量,其余特征值在保证样本协方差矩阵迹不变的情况下取相同值,将矩阵的求逆运算简化为向量的乘法运算。本发明提出的算法能够使运行时间明显少于基于特征值分解的最小方差算法,且对噪声具有良好的鲁棒性,成像效果明显优于传统的延时叠加算法、最小方差算法和波束域最小方差算法。
Description
技术领域
本发明属于超声成像领域,涉及一种融合特征值分解的低复杂度最小方差超声成像方法。
背景技术
超声成像中应用最为广泛的,也是最简单的波束形成技术即延时叠加算法(DelayAnd Sum,DAS),它是根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算,然后对延时后的数据对齐叠加。传统DAS算法复杂度低,成像速度快,但由于其采用固定窗函数加权导致主瓣宽度增加,分辨率较低。
近年来,为了提高波束形成算法的对比度和分辨率,自适应算法得到越来越多的研究。1969年Capon提出的最小方差(Minimum Variance,MV)波束形成算法是目前使用最为广泛的自适应算法。该方法依据保持期望方向增益不变,且使阵列输出能量达到最小的原则,通过动态地计算聚焦延时后的信号加权矢量,再将该矢量与输入信号相乘,提高了图像对比度和分辨率,但该算法的缺点是稳健性远不如传统的延时叠加算法,而且容易使有用信号相消,这在信噪比较低的情况下对图像质量有较大影响。因此,在最小方差算法的基础上算法分辨率、对比度和鲁棒性都还有很大的提升空间。
此外,自适应算法在鲁棒性和运行效率不及DAS算法,造成这些问题的主要原因是自适应算法涉及矩阵求逆和矩阵乘法运算,导致算法复杂度较高。设一个n×n维的矩阵,求逆的复杂度为O(n3),而传统DAS算法只有O(n);基于特征值分解的最小方差算法虽然具有较好的成像效果,但存在复杂度高、鲁棒性差的问题。
综上所述,急需发明一种融合特征值分解的低复杂度最小方差方法,在提高图像分辨率、对比度,能够保持算法稳健性、降低算法复杂度的波束形成算法,以全面整体提高超声成像质量。
发明内容
有鉴于此,本发明的目的在于提供一种融合特征值分解的低复杂度最小方差超声成像方法,该方法引入离散余弦变换构造转换矩阵,降低了样本协方差矩阵维数,并在矩阵求逆方面进行了改进,通过对特征值的特殊选取将矩阵求逆运算转换为向量乘法运算,降低矩阵求逆的复杂度等级,进一步提高了成像速度。
为达到上述目的,本发明提供如下技术方案:
一种融合特征值分解的低复杂度最小方差超声成像方法,包括以下步骤:
S1:对超声阵元接收的回波信号进行放大滤波处理,AD转换和延时处理,以获得超声回波数据;
S2:通过离散余弦变换得到转换矩阵,将接收阵列依次划分为L个具有重叠阵元的子阵,然后对相应接收子阵的回波信号进行前后向空间平滑处理;
S3:利用波束域转换矩阵将空间平滑后的子阵列回波信号转换到低维波束域,采用对角加载技术增加算法稳定性,得到波束域样本协方差矩阵估计;
S4:对波束域样本协方差矩阵进行特征值分解,提取信号子空间;
S5:噪声子空间对应的特征值在保证协方差矩阵迹不变的情况下取相同值,简化对角加载后的样本协方差矩阵的逆矩阵;
S6:利用简化后的样本协方差逆矩阵计算融合特征值分解的低复杂度的最小方差波束形成算法的最优权矢量;
S7:使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,得到自适应波束信号。
进一步,在步骤S2中通过离散余弦变换构造转换矩阵,对回波数据进行空间平滑处理,具体包括以下步骤:
S21:通过离散余弦变换构造(1+p)×L维波束域转换矩阵:
其中,矩阵T满足TTH=I,I为单位阵;Tm,n表示矩阵T第m行,第n列的值,L为每个子阵的阵元数,p表示协方差矩阵的降秩参数,满足p+1≤L降低样本协方差矩阵维数,[·]H为共轭转置运算;
S22:将N个阵元分为(N-L+1)个子阵列,其中每个子阵有L个阵元,设表示第l个阵元域子阵列接收的回波信号:
其中,N为阵元数,l表示第l个子阵列,表示第l个阵元在第k个采样点的回波数据,以此类推和分别表示第l+1个阵元和第l+L-1个阵元在第k个采样点的回波数据,[·]T表示矩阵转置运算。
进一步,在步骤S3中,利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,采用对角加载技术增加算法稳定性,得到波束域样本协方差矩阵估计具体包括以下步骤:
S31:利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,以第l个子阵列为例:
其中,为第l个子阵列对应的波束域回波数据,维数为(1+p)×1;表示第l个阵元域子阵列接收的回波信号,维数为L×1;T(1+p)×L表示波束域转换矩阵,维数为(1+p)×L;得到波束域子阵列数据后,利用公式求得波束域样本协方差矩阵Rb,其中为xb的转置共轭,E表示求取矩阵的期望;
S32:通过以下计算公式对波束域样本协方差矩阵进行对角加载,得到对角加载后的协方差矩阵增加算法稳定性:
其中,ε为对角加载系数,满足δ为常数,满足取
S33:,通过以下计算公式,得到波束域最小方差的最优权矢量为:
其中,为对角加载后的协方差矩阵,为的逆矩阵;ab=Ta为波束域方向向量,为ab的转置共轭,其中子阵的阵元数L=32,降秩参数p=8。
进一步,在步骤S4中,通过下式对波束域样本协方差矩阵进行特征值分解,提取信号子空间:
其中,Es=[e1,e2,···,eq]为信号子空间,q表示信号子空间的维数;En=[eq+1,eq+2,···,ep+1]为噪声子空间;为Es的共轭转置,为En的共轭转置;Λs=diag{λ1,λ2,···,λq},Λn=diag{λq+1,λq+2,···,λp+1};λi(i=1,2,···,p+1)为样本协方差矩阵的p+1个特征值,满足λ1≥λ2≥···≥λp+1,ei(i=1,2,···,p+1)为特征值λi对应的特征向量。
进一步,在步骤S5中,噪声子空间对应的特征值在保证协方差矩阵迹不变的情况下取相同值,简化对角加载后的样本协方差矩阵的逆矩阵;具体步骤如下:
S51:噪声子空间对应的特征值在协方差矩阵迹不变的情况下取相同值,保证超声回波信号能量的恒定,即:
其中表示求取矩阵的迹,即矩阵所有对角元素之和,q表示信号子空间的维数p为降秩参数,为对角加载后的协方差矩阵;
令则对角加载后的样本协方差矩阵的逆矩阵化简为:
其中,ei表示信号子空间和噪声子空间的向量,为ei的转置共轭矩阵;α-1表示I表示单位矩阵;
S52:进一步化简S51中的样本协方差矩阵的逆矩阵的求解过程,取q=1将矩阵的求逆运算转换为一次向量的乘法运算,如下式所示:
进一步,在步骤S6中,通过以下计算公式,将S52中计算的简化后的样本协方差矩阵的逆矩阵代入步骤S33中的波束域最小方差的最优权矢量中得到改进后的权矢量wib,将该权矢量向信号子空间进行投影,得到最优权矢量wibmv:
其中,Es为信号子空间,为Es的转置共轭;通过S5和S6步骤将降维后的矩阵求逆运算转换为向量的乘法运算。
进一步,在步骤S7中,使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,得到自适应波束信号:
其中,y(k)表示计算得到的自适应波束信号,N表示阵元数,L为每个子阵的阵元数;融合特征值分解的低复杂度最小方差超声成像方法得出的最优权矢量wibmv,表示wibmv的共轭转置,表示第l个子阵的输出向量,k表示第个k个采样点。
本发明的有益效果在于:本发明采用了融合特征值分解的低复杂度的最小方差波束形成超声成像方法,该方法首先通过离散余弦变换构造转换矩阵T,并对回波数据进行空间平滑处理;然后利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,降低算法复杂度,再通过对波束域样本协方差矩阵进行特征值分解,简化对角加载后的样本协方差矩阵的逆矩阵,再次降低计算复杂度;利用简化后的样本协方差矩阵计算特征分解的最优权矢量;使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,最终得到自适应波束信号。因此,本发明所提方法在提高图像分辨率、对比度的同时,能够保持算法稳健性、降低算法复杂度、减少运算时间加快成像速度,解决了基于特征值分解的最小方差算法虽然具有较好的成像效果,但复杂度高、鲁棒性差的问题。
附图说明
为了使本发明的目的、技术方案和有益效果更加清楚,本发明提供如下附图进行说明:
图1为本发明所述方法的流程图;
图2为前后向空间平滑算法示意图;
图3为5种算法点目标成像结果对比图;
图4为轴向距离40mm处5种算法分辨率对比图;
图5为轴向距离60mm处5种算法分辨率对比图;
图6为在20dB白噪声背景下5种算法点目标成像对比图;
图7为不同中心频率5种算法点目标成像结果对比图;
图8为5种算法吸声斑目标对比图;
图9为5种算法加噪声吸声斑目标图;
图10为5种算法geabr_0数据成像对比图;
图11为5种算法geabr_0数据成像12mm处分辨率对比图。
具体实施方式
下面将结合附图,对本发明的优选实施例进行详细的描述。
图1为本发明的方法流程图,如图所示,本发明提供一种融合特征值分解的低复杂度最小方差超声成像方法,包括以下步骤:
步骤S1:对回波信号进行放大和AD转换并进行延时聚焦处理,得到聚焦延时处理之后的信号x(k),x(k)表示为x(k)=[x1(k),x2(k),...,xN(k)],其中N表示超声阵列的阵元个数,k表示为对应采样深度的采样时刻。
步骤S2:通过离散余弦变换得到转换矩阵T,将接收阵列依次划分为一个具有重叠阵元的子阵,然后对相应接收子阵的回波信号进行前后向空间平滑处理。图2给出了前后向空间平滑算法示意图,具体包括以下步骤:
S21:通过离散余弦变换构造(1+p)×L维波束域转换矩阵:
矩阵T满足TTH=I,其中I为单位阵,Tm,n表示矩阵T第m行,第n列的值,L为每个子阵的阵元数,p表示协方差矩阵的降秩参数,满足p+1≤L降低样本协方差矩阵维数。
S22:将N个阵元分为(N-L+1)个子阵列,其中每个子阵含有L个阵元,设表示第l个阵元域子阵列接收的回波信号:
其中,N为阵元数,l表示第l个子阵列,表示第l个阵元在第k个采样点的回波数据,以此类推和分别表示第l+1个阵元和第l+L-1个阵元在第k个采样点的回波数据,[·]T表示矩阵转置运算。
步骤S3:利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,以第l个子阵列为例:
其中,为第l个子阵列对应的波束域回波数据,维数为(1+p)×1;表示第l个阵元域子阵列接收的回波信号,维数为L×1;T(1+p)×L表示波束域转换矩阵,维数为(1+p)×L;得到波束域子阵列数据后,利用公式求得波束域样本协方差矩阵Rb,其中为xb的转置共轭,E表示求取矩阵的期望;
S32:通过以下计算公式对波束域样本协方差矩阵进行对角加载,得到对角加载后的协方差矩阵增加算法稳定性:
其中,ε为对角加载系数,满足δ为常数,满足取
S33:,通过以下计算公式,得到波束域最小方差的最优权矢量为:
其中,为对角加载后的协方差矩阵,为的逆矩阵;ab=Ta为波束域方向向量,为ab的转置共轭;算法将矩阵求逆运算量由O(L3)降低为O((1+p)3),其中子阵的阵元数L=32,降秩参数p=8。
步骤S4:通过下式对波束域样本协方差矩阵进行特征值分解:
其中,Es=[e1,e2,···,eq]为信号子空间,q表示信号子空间的维数;En=[eq+1,eq+2,···,ep+1]为噪声子空间;为Es的共轭转置,为En的共轭转置。Λs=diag{λ1,λ2,···,λq},Λn=diag{λq+1,λq+2,···,λp+1};λi(i=1,2,···,p+1)为样本协方差矩阵的p+1个特征值,满足λ1≥λ2≥···≥λp+1,ei(i=1,2,···,p+1)为特征值λi对应的特征向量。
步骤S5:使噪声子空间对应的特征值在协方差矩阵的迹不变的情况下取相同值,保证超声回波信号能量的恒定的情况下,简化对角加载后的样本协方差矩阵的逆矩阵。具体步骤如下:
S51:噪声子空间对应的特征值在协方差矩阵迹不变的情况下取相同值,保证超声回波信号能量的恒定,即:
令则对角加载后的样本协方差矩阵的逆矩阵化简为:
S52:进一步化简S51中的样本协方差矩阵的逆矩阵的求解过程,取q=1将矩阵的求逆运算转换为一次向量的乘法运算,复杂度由O(L3)降低为O((p+1)2),如下式所示:
步骤S6:通过以下计算公式,将S52中计算的简化后的样本协方差矩阵的逆矩阵代入步骤S33中的波束域最小方差的最优权矢量中得到改进后的权矢量wib,将该权矢量向信号子空间进行投影,得到最优权矢量wibmv:
通过S5和S6步骤将降维后的矩阵求逆运算转换为向量的乘法运算,其复杂度也由O(L3)降低为O((1+p)2)。
步骤S7:使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,得到自适应波束信号:
其中,y(k)表示计算得到的自适应波束信号,N表示阵元数,L为每个子阵的阵元数;融合特征值分解的低复杂度最小方差超声成像方法得出的最优权矢量wibmv,表示wibmv的共轭转置,表示第l个子阵的输出向量,k表示第个k个采样点。
FieldII是丹麦理工大学基于声学原理开发的一款超声实验仿真平台,其在理论研究上获得了广泛的认可和使用。为验证所提算法的有效性,利用Field II对超声成像中常用的点散射目标和吸声斑目标进行成像并利用实际实验数据进行成像对比实验。在点目标仿真实验中,设置了15个目标点,轴向距离分布在30mm~80mm,每隔5mm设定目标点,其中40mm和60mm处分布3个目标点其余位置各有1个目标点。仿真采用线性阵列发射超声信号,发射为定点聚焦,而接收设为动态聚焦,成像的动态范围是60dB,成像方式为线扫描成像。同时,设一中心在25mm,半径为3mm的圆形区域吸声斑,外部随机分布着100000个散射点,仿真在接收回波中不加噪声和加入一定噪声的情况,并设定成像动态范围为80dB。实验所采用的阵元中心频率为2.5MHz,阵元数目为64个,间距为0.11mm,采样频率为50MHz,声速为1540m/s,设成像动态范围为60dB。对上述三个实验目标采用延时叠加算法(DAS),最小方差算法(MV),波束域最小方差算法(BMV),基于特征值分解的最小方差算法(ESBMV)以及融合特征值分解的低复杂度最小方差超声成像算法(IBMV)进行对比成像实验。
图3给出了5种算法点目标成像效果对比图,从图3中可以看出DAS算法成像质量最差,分辨率最低,相比于其他4种算法横向伪影最多,两个散射点已经相互干扰难以区分。其余自适应算法明显优于DAS算法,其中MV和BMV成像效果基本一致,IBMV和ESBMV成像效果直观上基本一致。
图4给出了轴向距离40mm处5种算法分辨率对比图,图5给出了轴向距离60mm处5种算法分辨率对比图。从图4、图5中可以看出,相比传统DAS算法,自适应算法在主瓣宽度上有较大提升,40mm处改善更为明显;从图3也可以直观看出,自适应算法对旁瓣峰值的抑制同样较为明显,其中MV和BMV算法无论在40mm或是60mm处主瓣宽度都相差无几,由分辨率图也可以看出旁瓣峰值基本一致,分辨率曲线几乎重合。对于IBMV算法,无论主瓣宽度和旁瓣峰值都优于MV和BMV算法,在40mm处,主瓣宽度优于ESBMV,但旁瓣峰值略高于ESBMV。
图6给出了在20dB白噪声背景下5种算法点目标成像对比图,使仿真更加贴近真检测环境。由图6可以看出,加入较强的噪声后,背景区域明显观察出噪声白斑的存在,但不同算法仍可以较为清晰的观察出目标点,说明融合特征值分解的低复杂度最小方差超声成像算法(IBMV)对比的算法对噪声均具有一定的鲁棒性,可以在噪声存在的情况下保持成像效果。分辨率曲线趋势与无噪声情况下基本一致。
图7给出了不同中心频率5种算法点目标成像结果对比图,在实际医用B型超声检测过程中,常因为检测部位的不同而选用不同中心频率的超声探头,应用于腹部的凸阵探头频率通常为2.5MHz、3.5MHz和5.0MHz,腔内探头一般是6.5MHz,血管探测的线性探头通常为7.5MHz,高频的线性探头能达到10MHz和12MHz。图7选取了4中有代表性的探头频率进行点目标成像仿真,声速设为1540m/s。由图7可以直观看出,IBMV算法优于ESBMV具有较少的侧向伪影;MV和BMV成像效果基本一致,优于传统DAS算法,不及ESBMV和IBMV算法;随着中心频率的增加,不同算法的成像效果均有提高,纵向伪影明显降低;与其他算法相比,IBMV能够在不同中心频率下保持更好的成像效果,适用于多种检测对象,具有较强的稳定性。
图8给出了5种算法吸声斑目标对比图。从图8看出,DAS成像效果最差,吸声斑区域有较多的噪声干扰,背景区域与中心暗斑成像对比度较差;自适应算法相比DAS算法成像效果有明显提升,其中ESBMV和IBMV算法优于MV和BMV算法,可以清晰地区分出吸声斑区域和背景噪声区域,成像对比度较大。为了更加直观的评估吸声斑的成像效果,引入对比度(contrast ratio,CR)和背景区域方差来衡量成像结果,其中对比度是指中心暗斑平均功率与背景区域平均功率之差的绝对值;选取中心暗斑区域为整个吸声斑区域,背景区域为整个散射点分布区域;背景区域方差表征算法的鲁棒性,其值越小算法鲁棒性越好,具体计算结果如表1所示:
表1算法对比度对比表
由表1可见,DAS算法成像效果最差,对比度仅为24.03dB,由于算法运算简单,所以鲁棒性最强,背景区域方差为7.18dB;而MV和BMV对比度有所提高,但由于涉及到矩阵的求逆运算,复杂度相比DAS有较大提升,所以算法的稳健性不及DAS算法;ESBMV和IBMV算法在MV的基础上增加了矩阵的特征值分解和权值投影运算,所以鲁棒性相比MV较差,但对比度在原有基础上提升明显,且IBMV算法的稳定性优于ESBMV。
图9给出了5种算法加噪声吸声斑目标图。从表2可以看出,中心暗斑平均功率和背景区域平均功率相比表1均有一定程度的提升,对比度下降明显;由于传统DAS算法鲁棒性最强,所以在有噪声情况下对比度优于自适应算法,本文提出的IBMV在对比度和稳健性方面均优于ESBMV算法。
表2加噪声对比度对比表
算法复杂度分析:
表3成像时间对比表
由表5可以看出,传统DAS算法无论是点目标成像还是吸声斑目标成像运行时间都远远优于自适应算法;BMV相比于MV算法由于降低了样本协方差矩阵的维数,运行时间有所提升,而ESBMV虽然在成像分辨率和对比度方面提高明显,但运行时间却远远多于其他算法,大约是MV运行时间的3倍。相比之下,IBMV算法虽然引入了矩阵特征值分解运算,但在矩阵维数和矩阵求逆方面做了较大改进,在成像效果与ESBMV算法基本一致的情况下,明显缩短了运行时间,大约为ESBMV运行时间的50%,且算法鲁棒性也略优于ESBMV;相比于MV算法,IBMV的成像分辨率和对比度有显著提升,但运行时间稍显不足。综合看来,本文提出的IBMV算法能够实现运行时间没有较大增加的情况下,大幅度提升成像效果;在成像效果基本一致的情况下,大幅度提升运行效率。从算法复杂度角度考虑,矩阵求逆和特征值分解复杂度为O(n3)、特征值排序和向量乘法复杂度为O(n2)、权值投影复杂度为O(n3+n2);传统DAS算法对所有阵元信号进行简单的叠加运算,复杂度为O(N);MV算法涉及子阵列回波信号协方差矩阵的求逆运算,复杂度为O(L3);ESBMV算法中包含协方差矩阵的求逆运算、特征值分解和排序运算以及权值投影运算,复杂度为O(3L3+2L2);BMV算法对协方差矩阵进行降维,其复杂度为O((p+1)3);IBMV算法将降维后的协方差矩阵求逆转换为向量乘法运算,同时涉及特征值分解和排序运算以及权值投影,复杂度为O(3(p+1)2+2(p+1)3);其中N=64是本文仿真的阵元数;L=32是空间平滑后子阵列的阵元数;p表示协方差矩阵的降秩参数,提取波束域信号低频部分保证算法的成像效果,实验验证p=8满足条件。
图10给出了5种算法geabr_0数据成像对比图;图11给出了5种算法geabr_0数据成像12mm处分辨率对比图;实验采用的具体参数如下:阵元中心频率为3.33MHz,阵元数目为64个,阵元间距为0.2413mm,采样频率为17.76MHz,声速为1500m/s;图8为上述5种算法的实验成像图,成像的动态范围为60dB。从图10可以看出,传统DAS算法成像效果最差,点目标受背景噪声干扰严重;ESBMV和IBMV算法成像结果明显优于其它算法,背景噪声区域较暗,点目标亮度较亮,相比于背景区域较为明显,成像对比度较大。为了更加直观的观察出不同算法的成像分辨率,截取轴向距离为12mm处的数据,做出不同算法的成像分辨率对比结果,如图11所示。由图9可见,MV和BMV算法成像分辨率相当,且主瓣宽度和旁瓣等级均低于传统DAS算法;ESBMV和IBMV算法在MV的基础上进一步抑制了旁瓣等级,提高了成像对比度;实验结果与之前仿真结果类似,验证了所提算法的有效性。
最后说明的是,以上优选实施例仅用以说明本发明的技术方案而非限制,尽管通过上述优选实施例已经对本发明进行了详细的描述,但本领域技术人员应当理解,可以在形式上和细节上对其作出各种各样的改变,而不偏离本发明权利要求书所限定的范围。
Claims (7)
1.一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:该方法包括以下步骤:
S1:对超声阵元接收的回波信号进行放大滤波处理,AD转换和延时处理,以获得超声回波数据;
S2:通过离散余弦变换得到转换矩阵,将接收阵列依次划分为L个具有重叠阵元的子阵,然后对相应接收子阵的回波信号进行前后向空间平滑处理;
S3:利用波束域转换矩阵将空间平滑后的子阵列回波信号转换到低维波束域,采用对角加载技术增加算法稳定性,得到波束域样本协方差矩阵估计;
S4:对波束域样本协方差矩阵进行特征值分解,提取信号子空间;
S5:噪声子空间对应的特征值在保证协方差矩阵迹不变的情况下取相同值,简化对角加载后的样本协方差矩阵的逆矩阵;
S6:利用简化后的样本协方差逆矩阵计算融合特征值分解的低复杂度的最小方差波束形成算法的最优权矢量;
S7:使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,得到自适应波束信号。
2.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S2中通过离散余弦变换构造转换矩阵,对回波数据进行空间平滑处理,具体包括以下步骤:
S21:通过离散余弦变换构造(1+p)×L维波束域转换矩阵:
其中,矩阵T满足TTH=I,I为单位阵;Tm,n表示矩阵T第m行,第n列的值,L为每个子阵的阵元数,p表示协方差矩阵的降秩参数,满足p+1≤L降低样本协方差矩阵维数,[·]H为共轭转置运算;
S22:将N个阵元分为(N-L+1)个子阵列,其中每个子阵有L个阵元,设表示第l个阵元域子阵列接收的回波信号:
其中,N为阵元数,l表示第l个子阵列,表示第l个阵元在第k个采样点的回波数据,以此类推和分别表示第l+1个阵元和第l+L-1个阵元在第k个采样点的回波数据,[·]T表示矩阵转置运算。
3.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S3中,利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,采用对角加载技术增加算法稳定性,得到波束域样本协方差矩阵估计具体包括以下步骤:
S31:利用波束域转换矩阵T将空间平滑后的子阵列回波信号转换到低维波束域,以第l个子阵列为例:
其中,为第l个子阵列对应的波束域回波数据,维数为(1+p)×1;表示第l个阵元域子阵列接收的回波信号,维数为L×1;T(1+p)×L表示波束域转换矩阵,维数为(1+p)×L;得到波束域子阵列数据后,利用公式求得波束域样本协方差矩阵Rb,其中为xb的转置共轭,E表示求取矩阵的期望;
S32:通过以下计算公式对波束域样本协方差矩阵进行对角加载,得到对角加载后的协方差矩阵增加算法稳定性:
其中,ε为对角加载系数,满足δ为常数,满足取
S33:,通过以下计算公式,得到波束域最小方差的最优权矢量为:
其中,为对角加载后的协方差矩阵,为的逆矩阵;ab=Ta为波束域方向向量,为ab的转置共轭,其中子阵的阵元数L=32,降秩参数p=8。
4.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S4中,通过下式对波束域样本协方差矩阵进行特征值分解,提取信号子空间:
其中,Es=[e1,e2,···,eq]为信号子空间,q表示信号子空间的维数;En=[eq+1,eq+2,···,ep+1]为噪声子空间;为Es的共轭转置,为En的共轭转置;Λs=diag{λ1,λ2,···,λq},Λn=diag{λq+1,λq+2,···,λp+1};λi(i=1,2,···,p+1)为样本协方差矩阵的p+1个特征值,满足λ1≥λ2≥···≥λp+1,ei(i=1,2,···,p+1)为特征值λi对应的特征向量。
5.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S5中,噪声子空间对应的特征值在保证协方差矩阵迹不变的情况下取相同值,简化对角加载后的样本协方差矩阵的逆矩阵;具体步骤如下:
S51:噪声子空间对应的特征值在协方差矩阵迹不变的情况下取相同值,保证超声回波信号能量的恒定,即:
其中表示求取矩阵的迹,即矩阵所有对角元素之和,q表示信号子空间的维数p为降秩参数,为对角加载后的协方差矩阵;
令则对角加载后的样本协方差矩阵的逆矩阵化简为:
其中,ei表示信号子空间和噪声子空间的向量,为ei的转置共轭矩阵;α-1表示I表示单位矩阵;
S52:进一步化简S51中的样本协方差矩阵的逆矩阵的求解过程,取q=1将矩阵的求逆运算转换为一次向量的乘法运算,如下式所示:
6.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S6中,通过以下计算公式,将S52中计算的简化后的样本协方差矩阵的逆矩阵代入步骤S33中的波束域最小方差的最优权矢量中得到改进后的权矢量wib,将该权矢量向信号子空间进行投影,得到最优权矢量wibmv:
其中,Es为信号子空间,为Es的转置共轭;通过S5和S6步骤将降维后的矩阵求逆运算转换为向量的乘法运算。
7.根据权利要求1所述的一种融合特征值分解的低复杂度最小方差超声成像方法,其特征在于:在步骤S7中,使用融合特征值分解的低复杂度的最小方差波束形成权值对采样信号进行加权求和,得到自适应波束信号:
其中,y(k)表示计算得到的自适应波束信号,N表示阵元数,L为每个子阵的阵元数;融合特征值分解的低复杂度最小方差超声成像方法得出的最优权矢量wibmv,表示wibmv的共轭转置,表示第l个子阵的输出向量,k表示第个k个采样点。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811243839.2A CN109187771B (zh) | 2018-10-24 | 2018-10-24 | 一种融合特征值分解的低复杂度最小方差超声成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811243839.2A CN109187771B (zh) | 2018-10-24 | 2018-10-24 | 一种融合特征值分解的低复杂度最小方差超声成像方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109187771A true CN109187771A (zh) | 2019-01-11 |
CN109187771B CN109187771B (zh) | 2020-12-04 |
Family
ID=64943118
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811243839.2A Expired - Fee Related CN109187771B (zh) | 2018-10-24 | 2018-10-24 | 一种融合特征值分解的低复杂度最小方差超声成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109187771B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110501711A (zh) * | 2019-08-15 | 2019-11-26 | 重庆大学 | 一种基于乘幂法的低复杂度最小方差超声成像方法 |
Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102499712A (zh) * | 2011-09-30 | 2012-06-20 | 重庆大学 | 一种基于特征空间的前后向自适应波束形成方法 |
KR20130124210A (ko) * | 2012-05-03 | 2013-11-13 | 삼성메디슨 주식회사 | 수신 빔포밍을 수행하는 초음파 시스템 및 방법 |
CN103536316A (zh) * | 2013-09-22 | 2014-01-29 | 华中科技大学 | 一种空时平滑相干因子类自适应超声成像方法 |
EP2846169A1 (en) * | 2013-09-10 | 2015-03-11 | Seiko Epson Corporation | Ultrasonic measurement apparatus, ultrasonic imaging apparatus, and ultrasonic measurement method |
CN105223567A (zh) * | 2015-09-28 | 2016-01-06 | 中国科学院声学研究所 | 一种应用于超声成像的稳健宽带自适应波束形成方法 |
CN105760892A (zh) * | 2016-03-10 | 2016-07-13 | 重庆大学 | 一种改进的最小方差超声成像方法 |
CN106510761A (zh) * | 2016-12-12 | 2017-03-22 | 重庆大学 | 一种信噪比后滤波与特征空间融合的最小方差超声成像方法 |
CN106814141A (zh) * | 2017-01-04 | 2017-06-09 | 天津大学 | 一种基于正交匹配追踪的相控阵超声信号压缩方法 |
CN107817297A (zh) * | 2017-11-23 | 2018-03-20 | 西安电子科技大学 | 一种基于超声回波rf数据的超声成像处理方法及处理系统 |
CN108309352A (zh) * | 2018-03-28 | 2018-07-24 | 东北大学 | 一种余弦变换域超声成像方法 |
-
2018
- 2018-10-24 CN CN201811243839.2A patent/CN109187771B/zh not_active Expired - Fee Related
Patent Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102499712A (zh) * | 2011-09-30 | 2012-06-20 | 重庆大学 | 一种基于特征空间的前后向自适应波束形成方法 |
KR20130124210A (ko) * | 2012-05-03 | 2013-11-13 | 삼성메디슨 주식회사 | 수신 빔포밍을 수행하는 초음파 시스템 및 방법 |
EP2846169A1 (en) * | 2013-09-10 | 2015-03-11 | Seiko Epson Corporation | Ultrasonic measurement apparatus, ultrasonic imaging apparatus, and ultrasonic measurement method |
CN103536316A (zh) * | 2013-09-22 | 2014-01-29 | 华中科技大学 | 一种空时平滑相干因子类自适应超声成像方法 |
CN105223567A (zh) * | 2015-09-28 | 2016-01-06 | 中国科学院声学研究所 | 一种应用于超声成像的稳健宽带自适应波束形成方法 |
CN105760892A (zh) * | 2016-03-10 | 2016-07-13 | 重庆大学 | 一种改进的最小方差超声成像方法 |
CN106510761A (zh) * | 2016-12-12 | 2017-03-22 | 重庆大学 | 一种信噪比后滤波与特征空间融合的最小方差超声成像方法 |
CN106814141A (zh) * | 2017-01-04 | 2017-06-09 | 天津大学 | 一种基于正交匹配追踪的相控阵超声信号压缩方法 |
CN107817297A (zh) * | 2017-11-23 | 2018-03-20 | 西安电子科技大学 | 一种基于超声回波rf数据的超声成像处理方法及处理系统 |
CN108309352A (zh) * | 2018-03-28 | 2018-07-24 | 东北大学 | 一种余弦变换域超声成像方法 |
Non-Patent Citations (3)
Title |
---|
ZHILIANG BAI,ET AL: "Ultrasonic Phased Array Compressive Imaging in Time and Frequency Domain: Simulation,Experimental Verification and Real Application", 《SENSORS》 * |
孟德明 等: "融合特征空间最小方差波束形成和广义相干系数的超声成像方法", 《生物医学工程研究》 * |
王平 等: "超声成像中基于特征空间的前后向最小方差波束形成", 《声学学报》 * |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110501711A (zh) * | 2019-08-15 | 2019-11-26 | 重庆大学 | 一种基于乘幂法的低复杂度最小方差超声成像方法 |
Also Published As
Publication number | Publication date |
---|---|
CN109187771B (zh) | 2020-12-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106510761B (zh) | 一种信噪比后滤波与特征空间融合的最小方差超声成像方法 | |
CN110501423B (zh) | 一种基于频域分段的高分辨率最小方差超声成像方法 | |
Xu et al. | Spatio-temporally smoothed coherence factor for ultrasound imaging [correspondence] | |
Zhao et al. | Subarray coherence based postfilter for eigenspace based minimum variance beamformer in ultrasound plane-wave imaging | |
Matrone et al. | Depth-of-field enhancement in filtered-delay multiply and sum beamformed images using synthetic aperture focusing | |
CN108670304B (zh) | 一种基于改进dmas算法的超声平面波成像方法 | |
CN102764139B (zh) | 基于特征空间分析和区域判别的医学超声波束形成方法 | |
CN108836389B (zh) | 平面波相关点相干自适应波束合成成像方法 | |
CN109164453A (zh) | 一种融合高度相干滤波器的最小方差超声成像方法 | |
EP3278734B1 (en) | Beamforming device, ultrasonic imaging device, and beamforming method allowing simple spatial smoothing operation | |
Agarwal et al. | Improving spatial resolution using incoherent subtraction of receive beams having different apodizations | |
CN111208213A (zh) | 融合交替乘子迭代的谱寻求子带最小方差超声成像算法 | |
Salari et al. | User parameter-free minimum variance beamformer in medical ultrasound imaging | |
CN108309352A (zh) | 一种余弦变换域超声成像方法 | |
Matrone et al. | Ultrasound synthetic aperture focusing with the delay multiply and sum beamforming algorithm | |
CN111856474A (zh) | 一种基于子阵的空时域条件相干系数超声成像方法 | |
CN110501711A (zh) | 一种基于乘幂法的低复杂度最小方差超声成像方法 | |
CN111278363B (zh) | 超声成像设备、系统及其超声造影成像的图像增强方法 | |
CN109187771B (zh) | 一种融合特征值分解的低复杂度最小方差超声成像方法 | |
Kozai et al. | Optimization of frequency and plane-wave compounding by minimum variance beamforming | |
CN112120730A (zh) | 一种基于混合子空间投影的广义旁瓣相消超声成像方法 | |
CN113625286B (zh) | 基于相干特征的强稳健性截断相干系数超声波束形成方法 | |
CN114415187A (zh) | 基于子空间斜投影的强稳健性最小方差超声波束形成方法 | |
Ziksari et al. | Fast beamforming method for plane wave compounding based on beamspace adaptive beamformer and delay-multiply-and-sum | |
CN113679418A (zh) | 一种基于压缩信噪比系数加权的超声成像方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201204 |
|
CF01 | Termination of patent right due to non-payment of annual fee |