Generalized sidelobe canceler beamforming combined with Eigenspace-Wiener postfilter for medical ultrasound imaging

2.1Synthetic aperture (SA)

Synthetic aperture (SA) beamforming is a widely used non-adaptive beamforming algorithm. It can realize transmitting and receiving focusing of beam [23], which provides a good foundation for the adaptive beamforming. Generally, the SA beamforming is implemented as follow: Firstly, in each scan, only one element is used to emit ultrasonic pulses, and all the array elements receive the echo signal. Moreover, the received data are processed by DS beamforming to realize the receiving focus. Then we repeat the above process on each single transmission element in the array and obtain a series of low-resolution images (LRI). Finally, all the low-resolution images are added and averaged to produce a high-resolution image (HRI), which achieves the transmitting focus.

For a linear transducer array, we assume it consists of M elements, the high-resolution image of SA can be expressed as:

where ωi is the weight power. In general, the weight is preset, and we often use window functions for it [24]. Because of the theory of acoustic reciprocity, there is no difference between the transmit (Tx) aperture focusing and receive (Rx) aperture focusing where the weight is applied [25]. And previous researches show that the adaptive beamforming is also effective in the Tx aperture [26]. Therefore, we can introduce a new adaptive beamforming to Eq. (1) to calculate an adaptive weight for the better image quality.

2.2Sensor signal model

Considering a linear transducer array consisting of M elements, the beamformer output is expressed by [27]:

where k is time instance and the weight vector 𝝎⁢(k)=[ω1⁢(k),ω2⁢(k),…,ωM⁢(k)]T. ∗ denotes complex conjugate, (⋅)T and (⋅)H denote transpose and conjugate transpose, respectively. 𝑿⁢(k) represents the received signal which has achieved time compensation, 𝑿⁢(k)=[x1⁢(k),x2⁢(k),…,xM⁢(k)]T. For the received signal 𝑿⁢(k), it mainly consists of the following two components:

where 𝑺⁢(k) denotes the desired signals which contain the useful information of detected points. 𝑷⁢(k) denotes the interference and noise signals including sidelobe signals and system noise which will cause artifacts and decrease the image quality. Therefore, the vector 𝝎⁢(k) should retain desired signals as much as possible, and filter out interference and noise 𝑷⁢(k). Since the adaptive beamforming can adjust the weight vector to changes of signals, the weight for each imaging point is optimal. Thus, the quality of echo images is improved.

2.3Generalized sidelobe canceler (GSC)

In adaptive beamforming, the weight vector should minimize the interference and noise power while the desired signal can pass without loss. This constraint is formed by [28]:

where 𝑹 denotes the covariance matrix of interference and noise signal, and 𝒂 represents the direction of arrival of the desired signals. In this paper, 𝒂 becomes a unit vector because the echo signal has been delayed and focused.

Figure 1.

Structure of generalized sidelobe canceler.

Figure 1 shows the structure of GSC beamforming. It separates weight vector 𝝎 into two parts 𝝎q and 𝝎a. According to the constraint of Eq. (4), we set 𝝎q parallel to the desired signal 𝑺⁢(k), so the desired signal can pass through the road without distortion. In contrast, for the road where 𝝎a located, a blocking matrix B is added to isolate desired signal. Thus, 𝝎a is only concerned with interference and noise signals 𝑷⁢(k). Finally, the desired signal is retained in the beamformer output, and the interference signal is removed. The weight vector 𝝎 is given by:

Through the operations above, the constraints of Eq. (4) can be redefined as:

And the answer of Eq. (2.3) is:

The blocking matrix 𝑩 is orthogonal to the desired signal, so it should satisfy [29]:

In practice, 𝑹 is unavailable because the interference and noise signals is unknown. Thus, the sample covariance matrix 𝑹^ is often used instead. To get a reliable matrix, some preprocessing method should be applied including the spatial smoothing and diagonal loading [30]. Through the spatial smoothing, the whole array is divided into several subarrays by overlapping smoothing. Then the sub-covariance matrices are calculated and averaged as the estimation of the sample covariance matrix [31]. Moreover, the diagonal loading adds a tiny white noise to the sample covariance matrix to improve the regularity [32]. With the help of the preprocessing method, the sample covariance matrix is given by:

where 𝑮p⁢(k) denotes the subarray and 𝑮p⁢(k)=[xp⁢(k),xp+1⁢(k),…,xp+L-1⁢(k)]. P denotes the number of subarrays, and L is the length of each subarray. The relationship between P and L is P=M-L+1 (L⩽M/2). ε determines the white noise energy added to 𝑹^, and Δ is an adjustable parameter which ranges from 10 to 100 as usual [33].

2.4Wiener postfilter

In medical ultrasound imaging, the Wiener postfilter is introduced to minimize the mean square error between the total output power of the beamformer and the desired signal power [34]. The response of the filter is calculated by:

The optimal solution Eq. (13) is

where 𝝎H⁢𝑹⁢𝝎 denotes the beamformer output power, |S|2 is the desired signal power, and |P|2 represents noise signal power. According to Eq. (14), the filter response H could be constructed, if we obtain the desired signal power and noise signal power. Some estimation methods has been proposed [21, 22], but the complexity and accuracy are unsatisfactory.

3.1Eigenspace decomposition

For the GSC beamforming, the sample covariance matrix 𝑹^ satisfies the following properties:

where 𝑬=[𝒆1,𝒆2,…,𝒆L], and ∧=diag[λ1,λ2,…,λL], λ1⩾λ2⩾…⩾λL. This processing is called the eigenspace decomposition, in which the λi is the eigenvalue, and 𝒆i is the corresponding eigenvector. In medical ultrasound imaging, the desired signal is in the same direction as 𝒆1, so the desired signal power is related to the eigenvector with the largest eigenvalue [35]. Thus, the span of the eigenvectors can be divided into two parts based on the magnitude of eigenvalues. The one consists of the eigenvectors with larger eigenvalues. For it contains the desired signal 𝑺⁢(k), we call it the signal space 𝑬s. The other is constructed by the remaining eigenvectors, which is called the noise space 𝑬p. Since the desired signal is related to the eigenvector which has the biggest eigenvalue λ1, we introduce a parameter δ, and the signal space 𝑬s can be constructed by the eigenvectors whose eigenvalues satisfy λi>δ⁢λ1. Then, the rest eigenvectors are used to construct the noise space 𝑬p.

Through the eigenspace decomposition, the sample covariance matrix 𝑹^ can be expressed as:

where ∧s and ∧p represent diagonal matrices composed of eigenvalues of signal space and noise space, respectively. 𝑹^s and 𝑹^p denote the covariance matrices of the corresponding space.

3.2Eigenspace-Wiener postfilter

According to Eq. (2) and (16), the beamformer output is given by:

Equation (17) indicates that the output energy of the beamformer can be calculated by the sample covariance and adaptive weight vector. What is more, the output energy can be further decomposed into signal energy 𝝎H⁢𝑹^s⁢𝝎 and noise energy 𝝎H⁢𝑹^p⁢𝝎. Based on Eq. (16), the covariance matrices of signal and noise are obtained by the eigenspace decomposition. Thus, we can use 𝝎H⁢𝑹^s⁢𝝎 and 𝝎H⁢𝑹^p⁢𝝎 as the estimators of |S|2 and |P|2 in Eq. (14), respectively. Then, the Eigenspace-Wiener postfilter can be expressed as:

3.3GSC with the Eigenspace-Wiener postfilter

To improve the performance of GSC, the adaptive weight is optimized by the Eigenspace-Wiener postfilter. And the new optimized weight can be expressed as:

With the new adaptive weight, the beamformer output is:

Comparing the Eqs (17) and (20), the first terms of the two equations are same, which represent the desired signal power, but the second terms are different which represent noise power. Specifically, for the echo signal from the point targets, the desired signal is much bigger than noise signal, so the noise power in Eqs (17) and (20) are all small, but our method can obtain a smaller result. For the region of anechoic lesion, there are few desired signals in the echo and it is full of noise. Since the noise signal power is much bigger than the desired signal power, a quantity of artifacts and noise will occur in the echo images. However, noise power in GSC-ESW is very small because the denominator of the second term of Eq. (20) is much bigger than 1. Therefore, GSC-ESW can make the output closer to the desired signals no matter how the noise signals change. As the noise is further suppressed, the GSC-ESW can achieve better lateral resolution and contrast.

4.1Simulated point targets

In this section, there are several point targets in the detection space, and the simulated data are processed by different algorithms. Figure 2 shows the imaging results. Figure 3 is the lateral variation curve of point targets at the depth of 40 mm and 45 mm. As shown in Fig. 3, GSC-ESW obtains the narrowest main lobe and lowest sidelobe between the four algorithms. Taking z= 40 mm as an example, the FWHM and PSL of different algorithms are listed in Table 1. It shows that GSC-ESW can reduce the FWHM by 85.5%, 80.5%, and 38.9% compared to DS, SA, and GSC beamformer. At the same time, the PSL of GSC-ESW is 75.91, 63.17, and 46.89 dB lower than that of DS, SA, and GSC beamformer. These results show that the GSC-ESW can improve the lateral resolution, and reduce the sidelobe energy level.

Figure 3.

Lateral variation of different algorithms at (a) z= 40 mm; (b) z= 45 mm.

Figure 4.

Simulated circular cyst target images with different algorithms. (a) DS; (b) SA; (c) GSC; (d) GSC-ESW. The dynamic range of image is 60 dB.

Figure 5.

Lateral variation at z= 40 mm.

4.2Simulated cyst targets

In this section, the echo data of a circular cyst target in a high speckle noise environment is simulated by Field II. The circular cyst is located at (x,z) = (0, 40) mm, and the radius is 5 mm. Within the cyst, there is no scattering point, and several scatterers are put outside the cyst to simulate the background noise. Figure 4 shows the beamforming response of DS, SA, GSC, and GSC-ESW. Figure 5 presents the lateral variation at the center of the cyst. The CR of the images are displayed in Table 2. GSC-ESW leads to 123.6%, 47.7%, 84.4% improvement in CR compared to DS, SA, and GSC. In addition, a lower intensity inside the cyst is achieved by GSC-ESW, indicating that the noise signal is effectively suppressed. This experiment indicates that the GSC-ESW beamformer performs well in contrast enhancement.

Figure 6.

Real data images with different algorithms. (a) DS; (b) SA; (c) GSC; (d) GSC-ESW. The dynamic range of image is 60 dB.

4.3Real echo data experiment

In this section, we use the real echo data to test algorithms. Figure 6 shows the beamforming response of different algorithms. Since the actual system is not as ideal as the simulated environment, more noise and errors will appear in the echo data [38], which make it difficult for beamformers to obtain accurate images. From Fig. 6a, we can see that the cyst is completely drowned out by the noise signal, while GSC-ESW leads to a clearly distinguishable results in Fig. 6d. The SA and GSC are also obviously disturbed by noise in Fig. 6b and c. Since the DS is difficult to distinguish the target, we do not carry out quantitative analysis on it. The lateral variation of SA, GSC, and GSC-ESW is shown in Fig. 7. According to quantitative calculation, the CR of SA, GSC and GSC-ESW are 8.56 dB, 6.25 dB, and 31.27 dB, respectively. This experiment shows that the new beamformer can still perform well in practical application.

Figure 7.

Lateral variation of real data at z= 88 mm.