CN105760892A - 一种改进的最小方差超声成像方法 - Google Patents

Abstract

本发明涉及一种改进的最小方差超声成像方法，属于超声成像技术领域；该方法首先对接收阵元的采样信号进行延时处理和前后向平滑处理，得到样本协方差矩阵估计；然后将前后向协方差矩阵估计进行特征值分解，构造信号子空间，同时利用椭球覆盖方向向量值域和限定方向向量模值对方向向量增加一对约束条件；在期望信号子空间中，基于最小方差准则，计算得到自适应波束形成权值；最后将自适应波束形成权值对经过前后向平滑处理的多路数据进行加权求和，从而得到一路自适应波束信号；该方法能够解决现有最小方差算法在图像分辨率、对比度以及对噪声鲁棒性等方面的问题，可以从整体上提高超声成像的质量。

Description

一种改进的最小方差超声成像方法

技术领域

本发明属于超声成像技术领域，涉及一种改进的最小方差超声成像方法。

背景技术

超声波因其无损、价格低廉且易于生成和控制等优点，在医学诊断领域得到了广泛的应用。如何提高超声图像质量是准确诊断疾病的前提，也是目前超声成像算法的研究重点。传统的延时叠加算法(DelayandSum,DAS)具有成像速度快的优点，是目前超声成像中最为广泛使用的一种波束形成算法。其主要原理是根据阵列与目标点的位置关系，对阵列中的每个阵元施加不同的延时，将接收的数据延时对齐后进行叠加。传统DAS算法成像质量较低，栅瓣等级高且对比度低。

国内外很多学者在DAS基础上引入了自适应加权算法提高图像质量。Capon算法也即最小方差算法(MinimumVariance,MV)是最常用的自适应加权算法，它是根据保持期望方向增益不变，并且使阵列输出能量最小的原则，计算出聚焦延时后信号的加权矢量。由于该方法是实时根据回波数据计算加权值，所以该算法相比于传统延时叠加算法可以有效的降低旁瓣等级，从而提高图像横向分辨率，也可以提高图像对比度，但该算法的缺点是稳健性远不如传统的延时叠加算法。

在超声成像中，根据最小方差准则计算出波束形成权值然后进行成像时，虽然图像分辨率和对比度较传统延时叠加算法有所提升，但算法鲁棒性下降，而且容易使有用信号相消，这在信噪比较低的情况下对图像质量有较大影响。因此在最小方差算法的基础上算法分辨率、对比度和鲁棒性都还有很大的提升空间。

综上所述，急需发明一种能够同时提高图像分辨率、对比度，并且还能保持算法的稳健性和抗噪能力，以全面提高超声图像的整体质量。

发明内容

有鉴于此，本发明的目的在于提供一种改进的最小方差超声成像方法，该方法能够同时提高图像分辨率、对比度以及波束形成鲁棒性的算法，有效克服了传统自适应波束形成算法稳健性低，不能显著提高图像对比度和分辨率等问题，全面提高了超声图像的整体质量。本发明的目的是研究一种改进的最小方差算法，提高超声成像质量。

为达到上述目的，本发明提供如下技术方案：

一种改进的最小方差超声成像方法，该方法包括以下步骤：

S1：对超声阵元接收的回波信号进行放大处理和A/D转换，以获得超声成像所需要的回波数据；

S2：将接收阵列依次划分为具有一个重叠阵元的子阵，然后对相应划分子阵接收的回波信号进行前后向平滑处理，以获得样本协方差矩阵；

S3：对样本协方差矩阵进行特征分解，构建信号子空间；

S4：利用椭球覆盖方向向量值域和限定方向向量模值增加方向向量一对约束条件；

S5：在期望信号子空间中，结合方向向量，根据最小方差准则，计算得到波束形成权值；

S6：将自适应波束形成权值对采样信号数据进行加权求和得到自适应波束信号。

进一步，在步骤S2中，具体包括：

S21：对采样信号进行聚焦延时处理，得到聚焦延时处理之后的信号x(k)，x(k)表示为x(k)＝[x₁(k),x₂(k),…,x_N(k)]，其中N表示超声阵列的阵元个数，k表示为对应采样深度的采样时刻；

S22：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向协方差矩阵估计

$\tilde{R}\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{R}_1\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H$

公式中表示第l个子阵的前向输出向量，为的共轭转置；

S23：定义为后向重叠向量，其中l＝1,2,…,N；类比S22，通过以下计算公式，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{\tilde{x}}_d^l\left(k\right){\tilde{x}}_d^l{\left(k\right)}^H$

公式中表示第l个子阵的后向输出向量，表示的共轭转置；

S24：通过计算前向协方差矩阵估计和后向协方差矩阵估计的求和平均得到前后向协方差矩阵估计

${\tilde{R}}_{FB}\left(k\right)=\frac{1}{2}\left(\tilde{R}\right(k)+{\tilde{R}}^b(k\left)\right)$

进一步，在步骤S3中，通过以下公式对进行特征分解：

${\tilde{R}}_{FB}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_i{e}_i{e}_{\stackrel{\·}{t}}^H-E_M{\mathrm{\Λ}}_M{E}_M^H$

其中，λ_i为的特征值，且按降序排列，λ₁≥λ₂≥…≥λ_N，e_i为λ_i对应的特征向量，为e_i的共轭转置，特征向量矩阵E_M＝[e₁…e_M]；为E_M的共轭转置，特征值矩阵Λ_M＝diag[λ₁…λ_M]；将矩阵划分为期望信号子空间及与之正交的噪声子空间：

${\tilde{R}}_{FB}=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H.$

其中Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量组成的信号子空间，E_n为较小特征值对应特征向量组成的噪声子空间，E_s ^H，E_n ^H为其对应共轭转置。

进一步，在步骤S4中，增加一对方向向量约束条件，具体包括以下两个条件：

1)：当方向向量存在偏差时，利用椭球覆盖方向向量的值域，即增加约束条件：

||a-a₁||²≤ε

其中a为假定的方向向量，a₁为期望信号方向向量，ε为误差边界；

2)：在此基础上，增加一个方向向量模值约束条件来提高本算法的鲁棒性：

||a₁||²＝M

其中，M为假定方向向量模值；

进一步，所述步骤S5中计算得到波束形成权值，具体步骤如下：

S51：将该算法简化为解决优化问题：

$\min f\left(a\right)=a^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^{H}a$ s.t||a-a₁||²≤ε,||a₁||²＝M

其中，a^H为方向向量的转置，为Λ_s的逆矩阵；利用拉格朗日算子法，求出期望信号方向向量的估计值：

$a_1=(N-\mathrm{\ϵ}/2)\frac{{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}{a^H{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}$

其中I为单位阵，再次使用拉格朗日算法可以求出算子γ的上限：

$\mathrm{\γ}\≤\frac{{\mathrm{\ξ}}_s}{{\left(M\mathrm{\ζ}\right)}^{1/2}-1}$

其中算子ξ_s为，ξ_s＝1/λ_s，λ_s为中较大特征值；算子ζ为，M为方向向量模值，ε为误差边界；

S52：通过最小方差原理计算自适应波束形成权值：

$w=E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1{\left(a_1^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1\right)}^{-1}$

公式中，a₁为期望信号方向向量，w为自适应波束形成权值，为Λ_s的逆矩阵。

进一步，所述步骤S6中，将改进的最小方差算法计算得到波束形成权值对所述的聚焦延时后的采样信号通过以下公式进行加权求和，计算得到自适应波束信号：

$y\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{w}^H\left(k\right)x_d^l\left(k\right)$

其中，y(k)表示计算得到的自适应波束信号，w^H表示w的共轭转置，表示第l个子阵的输出向量。

进一步，所述子阵元数目L的取值为L≤N/2。

本发明的有益效果在于：本发明采用了一种改进的最小方差超声成像算法，该算法将接收信号分解为期望信号子空间和噪声子空间，可在一定程度上解决信号相消的问题，针对期望信号子空间利用MV原则求解加权向量，另外在求解加权矢量时，对一般固定的方向矢量增加了一对约束条件，使算法对噪声的鲁棒性进一步增加。因此本发明提出的算法在图像分辨率、对比度和对噪声的鲁棒性等方面都有较大提高，克服了传统自适应波束形成算法对噪声鲁棒性低，不能显著提高图像对比度等问题。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法的流程图；

图2为前后向空间平滑算法示意图；

图3为全发全收下4种成像算法的点目标仿真结果；

图4为全发全收下3种算法在不同深度处横向分辨率；

图5为全发全收下4种成像算法的吸声斑仿真结果；

图6为合成孔径下4种成像算法点目标仿真结果；

图7为合成孔径下4种成像算法吸声斑仿真结果；

图8为40mm深度处波束截面图；

图9为4种成像算法实验数据成像结果；

图10为4种成像算法40mm深度处波束截面图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明的算法流程图，如图所示，本发明提供一种在超声成像中的改进最小方差算法，包括以下步骤：

步骤S1：对超声阵元接收的回波信号进行放大处理与A/D转换，以获得超声成像所需的回波数据；

步骤S2：将接收阵列依次划分为具有一个重叠阵元的子阵，然后对相应划分子阵接收的回波信号进行前后向平滑处理，以获得样本协方差矩阵。图2给出了前后向空间平滑算法示意图，具体包括以下步骤：

S21：对采样信号进行聚焦延时处理，得到聚焦延时处理之后的信号x(k)，x(k)表示为x(k)＝[x₁(k),x₂(k),…,x_N(k)]，其中N表示超声阵列的阵元个数，K表示为对应采样深度的采样时刻；

S22：把N个阵元依次划分为阵元数目为L的子阵，本发明的实例中子阵阵元数目L的取值为M/2，L的取值上限为M/2，当L＝M/2时图像的分辨率最高，稳健性较差。考虑到本发明已采用前后向空间平滑滤波提高算法稳健性，因此取L＝M/2，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向协方差矩阵估计

$\tilde{R}\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{R}_l\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H$

公式中表示第l个子阵的前向输出向量，为的共轭转置；

S23：定义为后向重叠向量，其中l＝1,2,…,N。类比S2，通过以下计算公式，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{\tilde{x}}_d^{l}d\left(k\right){\tilde{x}}_d^l{\left(k\right)}^H$

公式中表示第l个子阵的后向输出向量，表示的共轭转置；

S24：得到前向协方差矩阵估计和后向协方差矩阵估计后，通过求和平均得到前后向协方差矩阵估计

${\tilde{R}}_{FB}\left(k\right)=\frac{1}{2}\left(\tilde{R}\right(k)+{\tilde{R}}^b(k\left)\right)$

步骤S3：对样本协方差矩阵进行特征分解，构建信号子空间，可以通过以下公式对进行特征分解；

${\tilde{R}}_{FB}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_i{e}_i{e}_i^H-E_M{\mathrm{\Λ}}_M{E}_M^H$

其中，λ_i为的特征值，且按降序排列，λ₁≥λ₂≥…≥λ_N，e_i为λ_i对应的特征向量，为e_i的共轭转置，特征向量矩阵E_M＝[e₁…e_M]；为E_M的共轭转置，特征值矩阵Λ_M＝diag[λ₁…λ_M]；将矩阵划分为期望信号子空间及与之正交的噪声子空间：

${\tilde{R}}_{FB}=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H$

其中Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量组成的信号子空间，E_n为较小特征值对应特征向量组成的噪声子空间，E_s ^H，E_n ^H为其对应共轭转置。

用大于最大特征值0.5倍的特征值所对应的特征向量组成信号子空间E_s，而其余特征值对应的特征向量组成噪声子空间E_n。信号子空间中特征向量的个数的选取决定了主瓣宽度和旁瓣等级。反之特征向量的选取也与主瓣宽度和旁瓣等级有关，主瓣信号的能量主要集中在较大的特征值所对应的特征向量中，旁瓣信号主要集中在较小特征值所对应的特征值向量中。因此，一般用大于最大特征值的δ倍特征值所对应的特征向量组成信号子空间，δ取值范围为0到1之间，本实例中δ取0.5。

步骤S4：利用椭球覆盖方向向量值域和限定方向向量模值增加方向向量一对约束条件，具体包括以下步骤：

S41：由于方向向量存在偏差，利用椭球覆盖方向向量的值域，即增加约束条件：

||a-a₁||²≤ε

其中a为假定的方向向量，a₁为期望信号方向向量，ε为误差边界；

S42：在此基础上，增加一个方向向量模值约束条件来提高本算法的鲁棒性：

||a₁||²＝M

其中，M为假定方向向量模值；

步骤S5：在期望信号子空间内，结合方向向量，根据最小方差准则，计算得到波束形成权值，具体包括以下步骤；

S51：将本实例问题简化为解决优化问题：

$\min f\left(a\right)=a^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^{H}a$ s.t||a-a₁||²≤ε,||a₁||²＝M

其中，a^H为方向向量的转置，为Λ_s的逆矩阵；利用拉格朗日算子法，求出期望信号方向向量的估计值：

$a_1=(N-\mathrm{\ϵ}/2)\frac{{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}{a^H{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}$

其中I为单位阵，再次使用拉格朗日算法可以求出算子γ的上限：

其中算子ξ_s为，ξ_s＝1/λ_s，λ_s为中较大特征值；算子ζ为，M为方向向量模值，ε为误差边界；

S52：通过最小方差原理计算自适应波束形成权值：

$w=E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1{\left(a_1^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1\right)}^{-1}$

公式中，a₁为期望信号方向向量，w为自适应波束形成权值，为Λ_s的逆矩阵。

步骤S6：通过以下公式将自适应波束形成权值对采样信号数据进行加权求和得到自适应波束信号：

$y\left(k\right)=\frac{1}{N-L+1}\underset{i=1}{\overset{N-L+1}{\Σ}}{w}^H\left(k\right)x_d^l\left(k\right)$

其中，y(k)表示计算得到的自适应波束信号，w^H表示w的共轭转置，表示第l个子阵的输出向量。

为了验证本发明的有效性，在本实例中，利用FieldII对医学成像中常用的点散射目标和斑散射目标进行成像并对体膜进行实际数据采集。FieldII是基于线性系统空间响应原理，它的仿真结果与实际成像很接近，已被国际上广泛认同为仿真超声系统的标准。设置14个目标点，纵向间隔为5mm，横向间隔为2mm，深度分布在35mm～65mm深度处，并在回波数据中加入20dB的噪声，图像成像的动态范围为60dB。同时设定半径为3mm的吸声斑，深度在37mm～43mm之间，并在回波数据中加入一定强度的高斯白噪声，设定成像的动态范围为60dB。分别用全发全收和合成孔径进行4种算法成像，并比较各种成像算法的分辨率、对比度和对噪声的鲁棒性。体膜数据采集中心频率为f₀＝3.5MHz，采样频率为f_s＝25MHz。阵元个数N＝16，阵元间距为0.78mm，所成图像动态范围为60dB，采用4种成像算法并比较成像效果。

图3给出了全发全收下4种成像算法的点目标仿真结果。图4给出了在全发全收下三种算法在不同深度处横向分辨率。从图3直观来看，DAS算法成像效果质量最差，旁瓣等级最高，横向分辨率不如传统MV算法，ESBMV算法相较于传统MV算法有一定程度的提高，近场区域旁瓣等级更低，本文发明的IMV算法横向分辨率最好。同时，可以看出本文发明的IMV算法成像质量在MV以及ESBMV算法基础上有较大改进，图像上噪声明显较少，对噪声的鲁棒性有较大改善。图4中，(a)-6dB处散射点的分辨率，(b)-20dB处散射点的分辨率，从图4可以看出，随着深度增加，3种成像算法在分辨率上均有不同程度的降低，但本文发明的IMV算法分辨率仍高于MV算法。在不同深度时，本文发明的IMV算法-6dB处分辨率在MV算法基础上提高了一倍左右，远优于DAS算法，-20dB处本文发明的IMV算法仍优于MV。

图5给出了全发全收下4种成像算法的吸声斑仿真结果，表1给出了全发全收下4种算法对比度。从图5的结果可以看出本文发明IMV算法图像中的噪声可见度最低，中心暗斑处的噪声最少，对噪声鲁棒性最好。从表1可以看出MV算法虽然中心暗斑处功率相较于DAS算法提高了近5.8dB左右，但背景区域平均功率也有所上升，对比度在DAS基础上只提高了1dB左右；ESBMV算法的中心暗斑处功率及背景功率在MV基础上均有所提高，但总体对比度在MV基础上有提高。本文发明算法中心平均功率最大，在DAS、MV和ESBMV的基础上分别提高了12.8dB、7dB、5dB左右，同时背景区域的平均功率相比于DAS只上升了3dB左右，且低于MV与ESBMV的背景功率，对比度在DAS基础上提高了近10dB，在MV及ESBMV基础上分别提高了9dB、8dB左右。

表1全发全收下4种算法的对比度

成像算法	中心平均功率/dB	背景平均功率/dB	对比度/dB6 -->
				DAS	-33.74	-19.45	14.29
MV	-39.53	-24.16	15.37
				ESBMV	-41.66	-25.33	16.33
IMV	-46.62	-22.35	24.27

图6给出了合成孔径下4种成像算法点目标仿真结果，图7给出了合成孔径下4种成像算法吸声斑仿真结果，表2给出了合成孔径下4种算法对比度。从图6、7可以看出，在合成孔径模式下，各种算法的成像质量相对于全发全收模式均有不同程度的提高。在点目标情形下，体现为对噪声的鲁棒性有进一步的提高。在吸声斑仿真时，表现在中心暗斑处功率的上升，从图像上看，本专利发明IMV算法图像噪声可见度进一步降低，中心暗斑功率相对于DAS、MV、ESBMV最大，噪声鲁棒性最好。

从表2可以看出，在合成孔径下各种算法中心暗斑处的平均功率都有进一步的提高，本专利所提算法提高了近16dB，但背景区域平均功率也上升最多。MV算法中心平均功率略优于全发全收，同时背景区域平均功率有所下降，但对比度在DAS基础上并无改善。ESBMV算法中心功率略优于DAS、MV算法，背景区域功率相较于全发全收模式有所降低，但对比度提高不明显。本专利所提算法中心功率最大，在DAS、MV、ESBMV基础上分别提高了20dB、21dB、19dB左右，分辨率在DAS基础上提高了9dB，在ESBMV基础上提高了近8.5dB。

表2合成孔径下4种算法的对比度

成像算法	中心平均功率/dB	背景平均功率/dB	对比度/dB
				DAS	-41.86	-20.16	21.70
MV	-41.23	-19.83	21.40
				ESBMV	-43.40	-20.92	22.48
IMV	-62.31	-31.35	30.96

图8给出了40mm深度处波束截面图。图8中，(a)PW模式波束截面图(b)，SA模式波束截面图，从图8(a)可以看出，在全发全收模式下，MV、ESBMV算法旁瓣等级相比于DAS有一定程度的下降，但本专利所提的IMV算法的旁瓣等级最低。从图8(b)可以看出，合成孔径模式下，MV算法的中心暗斑平均功率与DAS相差不大，ESBMV算法在MV基础上有一定提高，本专利所提的IMV算法的旁瓣等级最低，且中心暗斑处的功率最优。

图9给出了4种成像算法实验数据成像结果，图10给出了4种成像算法40mm深度处波束截面图。由图9可以看出，本专利所提IMV算法的成像质量明显优于DAS，MV以及ESBMV，所成图像的对比度有明显提高，而且噪声被明显抑制。相比于MV及ESBMV算法，IMV算法在远场区域噪声的抑制能力更强。从图10中可以看出，MV与ESBMV算法两者的分辨率相当，但在传统的DAS基础上均有一定程度的提高。本专利所提IMV算法在三者中表现最优，分辨率和对比度均有提高，同时可以看出所提算法噪声抑制能力相较DAS、MV、ESBMV算法有较大提高。因此，本专利所提IMV算法具有更高分辨率、对比度以及对噪声的鲁棒性。

最后说明的是，以上优选实例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (7)

1.一种改进的最小方差超声成像方法，其特征在于：该方法包括以下步骤：

S1：对超声阵元接收的回波信号进行放大处理和A/D转换，以获得超声成像所需要的回波数据；

S2：将接收阵列依次划分为具有一个重叠阵元的子阵，然后对相应划分子阵接收的回波信号进行前后向平滑处理，以获得样本协方差矩阵；

S3：对样本协方差矩阵进行特征分解，构建信号子空间；

S4：利用椭球覆盖方向向量值域和限定方向向量模值增加方向向量一对约束条件；

S5：在期望信号子空间中，结合方向向量，根据最小方差准则，计算得到波束形成权值；

S6：将自适应波束形成权值对采样信号数据进行加权求和得到自适应波束信号。

2.根据权利要求1所述的一种改进的最小方差超声成像方法，其特征在于：在步骤S2中，具体包括：

S21：对采样信号进行聚焦延时处理，得到聚焦延时处理之后的信号x(k)，x(k)表示为x(k)＝[x₁(k),x₂(k),…,x_N(k)]，其中N表示超声阵列的阵元个数，k表示为对应采样深度的采样时刻；

S22：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向协方差矩阵估计

$\tilde{R}\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{R}_l\left(k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H$

公式中表示第l个子阵的前向输出向量，为的共轭转置；

S23：定义为后向重叠向量，其中l＝1,2,…,N；类比S22，通过以下计算公式，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{N-L+1}\underset{i=1}{\overset{N-L+1}{\Σ}}{\tilde{x}}_d^l\left(k\right){\tilde{x}}_d^l{\left(k\right)}^H$

公式中表示第l个子阵的后向输出向量，表示的共轭转置；

S24：通过以下计算公式计算前向协方差矩阵估计和后向协方差矩阵估计的求和平均，得到前后向协方差矩阵估计

${\tilde{R}}_{FB}\left(k\right)=\frac{1}{2}\left(\tilde{R}\right(k)+{\tilde{R}}^b(k\left)\right).$

3.根据权利要求2所述的一种改进的最小方差超声成像方法，其特征在于：

在步骤S3中，通过以下公式对进行特征分解：

${\tilde{R}}_{FB}=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_i{e}_i{e}_i^H=E_M{\mathrm{\Λ}}_M{E}_M^H$

其中，λ_i为的特征值，且按降序排列，λ₁≥λ₂≥…≥λ_N，e_i为λ_i对应的特征向量，为e_i的共轭转置，特征向量矩阵E_M＝[e₁…e_M]；为E_M的共轭转置，特征值矩阵Λ_M＝diag[λ₁…λ_M]；将矩阵划分为期望信号子空间及与之正交的噪声子空间：

${\tilde{R}}_{FB}=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H$

其中Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量组成的信号子空间，E_n为较小特征值对应特征向量组成的噪声子空间，E_s ^H，E_n ^H分别为E_s和E_n对应的共轭转置。

4.根据权利要求3所述的一种改进的最小方差超声成像方法，其特征在于：在步骤S4中，增加一对方向向量约束条件，具体包括以下两个条件：

1)：当方向向量存在偏差时，利用椭球覆盖方向向量的值域，即增加约束条件：

||a-a₁||²≤ε

其中a为假定的方向向量，a₁为期望信号方向向量，ε为误差边界；

2)：在此基础上，增加一个方向向量模值约束条件来提高算法的鲁棒性：

||a₁||²＝M

其中，M为假定方向向量模值。

5.根据权利要求4所述的一种改进的最小方差超声成像方法，其特征在于：所述步骤S5中计算得到波束形成权值，具体步骤如下：

S51：将该算法简化为解决优化问题：

$\min f\left(a\right)=a^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^{H}a$ s.t||a-a₁||²≤ε,||a₁||²＝M

其中a^H为方向向量a的转置，为Λ_s的逆矩阵；利用拉格朗日算子法，求出期望信号方向向量的估计值：

$a_1=(N-\mathrm{\ϵ}/2)\frac{{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}{a^H{(E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H+\mathrm{\γ}I)}^{-1}a}$

其中，I为单位阵，再次使用拉格朗日算法可以求出算子γ的上限：

$\mathrm{\γ}\≤\frac{{\mathrm{\ξ}}_s}{{\left(M\mathrm{\ζ}\right)}^{1/2}-1}$

其中算子ξ_s为，ξ_s＝1/λ_s，λ_s为中较大特征值；算子ζ为，M为方向向量模值，ε为误差边界；

S52：通过最小方差原理计算自适应波束形成权值：

$w=E_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1{\left(a_1^H{E}_s{\mathrm{\Λ}}_s^{-1}{E}_s^H{a}_1\right)}^{-1}$

公式中，a₁为期望信号方向向量，w为自适应波束形成权值，为Λ_s的逆矩阵。

6.根据权利要求5所述的一种改进的最小方差超声成像方法，其特征在于：所述步骤S6中，将改进的最小方差算法计算得到波束形成权值对所述的聚焦延时后的采样信号通过以下公式进行加权求和，计算得到自适应波束信号：

$y\left(k\right)=\frac{1}{N-L+1}\underset{i=1}{\overset{N-L+1}{\Σ}}{w}^H\left(k\right)x_d^l\left(k\right)$

其中，y(k)表示计算得到的自适应波束信号，w^H表示w的共轭转置，表示第l个子阵的输出向量。

7.根据权利要求2至6中任一项所述的一种改进的最小方差超声成像方法，其特征在于：所述子阵元数目L的取值为L≤N/2。