Block Iterative STMV Algorithm and Its Application in Multi-Targets Detection

Daizhu Zhu, Haoquan Guo, Yuanao Wei, Kaiju Wang

Shanghai Marine Electronic Equipment Research Institute, Shanghai, China.

DOI: 10.4236/jamp.2020.87103                

Abstract

STMV beamforming algorithm needs inversion operation of matrix, and its engineering application is limited due to its huge computational cost. This paper proposed block iterative STMV algorithm based on one-phase regressive filter, matrix inversion lemma and inversion of block matrix. The computational cost is reduced approximately as 1/4 M times as original algorithm when array number is M. The simulation results show that this algorithm maintains high azimuth resolution and good performance of detecting multi-targets. Within 1 - 2 dB directional index and higher azimuth discrimination of block iterative STMV algorithm are achieved than STMV algorithm for sea trial data processing. And its good robustness lays the foundation of its engineering application.

Keywords

Beamforming, Steered Minimum Variance (STMV), Block Iterative, Computational Cost, Multi-Targets Detection

1. Introduction

The bearing spectrum estimation is very important in sonar and radar fields. The original algorithm of bearing spectrum estimation based on arrays is conventional beam-forming (CBF), whose azimuth resolution is restricted by space “Fourier threshold”, often termed “Rayleigh threshold” [1] [2]. There are many kinds of high-resolution bearing spectrum estimation algorithms since 1970s and Capon proposed the minimum variance distortion response (MVDR) beam-forming algorithm. MVDR has two manifestations when dealing with wide-band signals. One is incoherent signal-subspace processing method (ISM) proposed by Wax, whose azimuth resolution decreases when dealing with correlated sources. The other one is coherent signal-subspace processing method (CSM) proposed by Wang and Kaveh. Importing the beam-focused idea, this algorithm’s estimation performance is better than ISM and can deal with correlated sources. Based on CSM, Krolik and Swinger proposed the steered minimum variance (STMV) algorithm, which can obtain wide-band gains and can deal with correlated sources directly.

Although STMV algorithm has good performance, it is not widely used for a long time because its huge computational cost has a high command for hardware performance when solving inversion operation of matrix, which restricts its practical application. Swingler [3] proposed sub-array STMV algorithm. Papers [4] [5] [6] [7] make further research on sub-array STMV algorithm. Zhu [8] proposed iterative STMV algorithm based on one-phase regressive filter and matrix inversion lemma which can reduce the computational cost to 2/M of STMV algorithm (M refer to array number) and meanwhile preserves perfectly the characters of STMV algorithm. Inspired by this idea, this paper proposes block iterative STMV algorithm based on one-phase regressive filter, matrix inversion lemma and inversion of block matrix, which can reduce computational cost further. Simulation experiments show that this iterative algorithm maintains high azimuth resolution of STMV and advantage of multi-targets detection, besides reducing computational cost markedly. The processing results of actual experimental data are also verified this algorithm has a good stability.

2. Steered Minimum Variance Beamforming

2.1. Basic Principles of STMV Algorithm

The basic principle of STMV based on STCM is keeping the response at the steered bearing angle $\theta$ as 1 by rotating the beam and minimizing the total energy of output.

Assuming that each array output is $x\left(t\right)={[x_1(t),x_2(t),\cdots ,x_M(t)]}^T$. Correspondingly, frequency domain output is $X\left(f_k\right)$. For CBF, according to the setting azimuth, each array output is inserted the time-delay ${\tau }_m\left(\theta \right)=md\cos \left(\theta \right)/c$, $m=0,\mathrm{...},M-1$. The output after time-delay is

$Y\left(f_k,\theta \right)={\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)$ (1)

where M is the array number, $\theta$ is the bearing angle, d is the interval between the adjacent array, c is the acoustic speed, ${\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)=diag\left(\left[1,e^{j2\pi f_{k}d\cos \left(\theta \right)/c},\cdots ,e^{j2\pi f_k\left(M-1\right)d\cos \left(\theta \right)/c}\right]\right)$ is the steered matrix.

Assuming that the weight vector of each array is $w={[w_1,w_2,\cdots ,w_M]}^T$, the output of array is

$y\left(t,\theta \right)=w^H{\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)e^{j2\pi f_{k}t}}$ (2)

where, h is the highest frequency point, and l is the lowest one.

The average power of array output is

$\begin{array}{l}P\left(\theta \right)=E\left[y\left(t,\theta \right)y{\left(t,\theta \right)}^H\right]\\ =w^{H}E\left[\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)}\right){\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)}\right)}^H\right]w\end{array}$ (3)

Considering that $X\left(f_k\right)$ and $X\left(f_m\right)$ are uncorrelated when T is large enough and $k\ne m$. Therefore

$\begin{array}{c}P\left(\theta \right)=w^{H}E\left[\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)}\right){\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)}\right)}^H\right]w\\ =w^{H}E\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\displaystyle \underset{m=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)X\left(f_k\right)X^H\left(f_m\right){\stackrel{\dot{}}{T}}_p^H\left(f_m,\theta \right)}}\right)w\\ =w^H\left({\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)E\left[X\left(f_k\right)X^H\left(f_k\right)\right]{\stackrel{\dot{}}{T}}_p^H\left(f_k,\theta \right)}\right)w\\ =w^H{R}_{stmv}\left(\theta \right)w\end{array}$ (4)

Estimating the cross-spectral density matrix (CSDM) $R\left(f_k\right)$ using the frequency domain output $X_n\left(f_k\right)$ of N snapshots, shown as Equation (5)

$R\left(f_k\right)=E\left[X\left(f_k\right)X^H\left(f_k\right)\right]={\displaystyle \underset{n=1}{\overset{N}{\sum }}{X}_n\left(f_k\right)X_n^H\left(f_k\right)}$ (5)

Defining $R_{stmv}\left(\theta \right)$ as the steered covariance matrix (STCM) corresponding to the bearing angle $\theta$, shown as Equation (6)

$R_{stmv}\left(\theta \right)={\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)R\left(f_k\right){\stackrel{\dot{}}{T}}_p^H\left(f_k,\theta \right)}$ (6)

The relationship between STCM and CSDM is also given by the formula.

The best weight vector of beam-former is to keep the response as 1 at the specified direction, meanwhile minimize the power of array output, i.e.

$\{\begin{array}{l}w=\arg \left[\min \left(w^H{R}_{stmv}^{}w\right)\right]\\ w^H{L}_M=1\end{array}$ (7)

where, $L_M$ is $M\times 1$ vector of ones.

The weight vector can be obtained as Equation (8) by applying the method of Lagrange multiplier.

$w=\frac{R_{stmv}^{-1}\left(\theta \right)L_M}{L_M{}^H{R}_{stmv}^{-1}\left(\theta \right)L_M}$ (8)

Substituting Equation (8) into Equation (4),

$P_{stmv}\left(\theta \right)=\frac{1}{L_M{}^H{R}_{stmv}^{-1}\left(\theta \right)L_M}$ (9)

STCM comprehensively utilizes information of each frequency and the matrix is invertible if the product of frequency points and snapshots is bigger than order of the matrix.

2.2. Principles of Iterative STMV Algorithm

When STMV beams are formed, the average output in frequency domain can be obtained through summing up N outputs $X_n\left(f_k\right)$ in frequency domain.

$\bar{X}\left(f_k\right)={\displaystyle \underset{n=1}{\overset{N}{\sum }}{X}_n\left(f_k\right)}$ (10)

Then STCM $R_{stmv}\left(\theta \right)$ can be represented as

$R_{stmv}\left(\theta \right)={\displaystyle \underset{k=l}{\overset{h}{\sum }}{\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)\bar{X}\left(f_k\right){\bar{X}}_{}^H\left(f_k\right){\stackrel{\dot{}}{T}}_p^H\left(f_k,\theta \right)}$ (11)

Defining that $U={\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)\bar{X}\left(f_k\right)$, Equation (11) can be represented as

$R_{stmv}\left(\theta \right)={\displaystyle \underset{k=l}{\overset{h}{\sum }}U{U}^H}$ (12)

In this algorithm, the former data are utilized sufficiently and one-phase regressive filter is used to integrate the covariance matrix $R_{stmv}\left(\theta \right)$, then STCM ${\hat{R}}_{stmv}^{}(n)$ estimated by using the nth snapshot can be represented as

${\hat{R}}_{stmv}(n)=\alpha {\hat{R}}_{stmv}\left(n-1\right)+(1-\alpha ){\displaystyle \underset{k=l}{\overset{h}{\sum }}U{U}^H}$ (13)

By matrix inversion lemma,

${(A+x{y}^H)}^{-1}=A^{-1}-\frac{A^{-1}x{y}^H{A}^{-1}}{1+y^H{A}^{-1}x}$ (14)

Therefore, ${\hat{R}}_{stmv}^{-1}(n)$ can be represented as

$\begin{array}{l}{\hat{R}}_{stmv}^{-1}(n)={[}^{\alpha }\\ ={[}^{\alpha }\\ ={\tilde{R}}_{stmv}^{-1}(h-1)-\frac{(1-\alpha ){\tilde{R}}^{-1}{}_{stmv}(h-1)U\left(f_h\right)U{\left(f_h\right)}^H{\tilde{R}}^{-1}{}_{stmv}(h-1)}{1+(1-\alpha )U{\left(f_h\right)}^H{\tilde{R}}^{-1}{}_{stmv}(h-1)U\left(f_h\right)}\end{array}$ (15)

where ${\tilde{R}}^{}{}_{stmv}(h-1)=\alpha {\hat{R}}_{stmv}\left(n-1\right)+(1-\alpha ){\displaystyle \underset{k=l}{\overset{h-1}{\sum }}U{U}^H}$. Obviously, in order to solve

${\hat{R}}_{stmv}^{-1}(n)$, we should solve ${\tilde{R}}^{-1}{}_{stmv}(h-1)$ first. Noting that ${\tilde{R}}^{}{}_{stmv}(h-1)$ and ${\hat{R}}_{stmv}^{}(n)$ have the same form, hence in order to solve ${\tilde{R}}^{-1}{}_{stmv}(h-1)$, we should solve ${\tilde{R}}^{-1}{}_{stmv}(h-2)$ first. By analogy, we can solve ${\hat{R}}_{stmv}^{-1}(n)$ by iterative as long as ${\tilde{R}}^{-1}{}_{stmv}(l)$ is solved. ${\tilde{R}}^{-1}{}_{stmv}(l)$ can be represented as

$\begin{array}{l}{\tilde{R}}^{-1}{}_{stmv}(l)={[}^{\alpha }\\ =\frac{1}{\alpha }[{\hat{R}}^{-1}{}_{stmv}(n-1)-\frac{(1-\alpha ){\hat{R}}^{-1}{}_{stmv}(n-1)U\left(f_l\right)U^{\text{H}}\left(f_l\right){\hat{R}}^{-1}{}_{stmv}(n-1)}{\alpha +(1-\alpha )U^{\text{H}}\left(f_l\right){\hat{R}}^{-1}{}_{stmv}(n-1)U\left(f_l\right)}]\end{array}$ (16)

Utilizing Equations (15), (16) comprehensively, ${\hat{R}}_{stmv}^{-1}(n)$ can be estimated by $h-l$ iterations.

2.3. Principles of Block Iterative STMV Algorithm

2.3.1. Inversion Algorithm of Block Matrix and Complexity Analysis

In [9], the inversion formula of $2\times 2$ block matrix has been elaborated.

Lemma Assuming that $A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{12}^H& A_{22}\end{array}\right]$, $A_{11}$ and $A_{22}$ are m order Hermitian matrices, $A_{12}$ is m order matrix, then

$A^{-1}=\left[\begin{array}{cc}{A}_{11}^{-1}+A_{11}^{-1}{A}_{12}{\left(A_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}\right)}^{-1}{\left(A_{11}^{-1}{A}_{12}\right)}^H& -A_{11}^{-1}{A}_{12}{\left(A_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}\right)}^{-1}\\ {\left(-A_{11}^{-1}{A}_{12}{\left(A_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}\right)}^{-1}\right)}^H& {\left(A_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}\right)}^{-1}\end{array}\right]$ (17)

Noting that computing cost of multiplication is tens of times of add operation, so this paper only takes computational cost of multiplication into consideration. Assuming inversion computational cost of m order matrix, then inversion computational cost of $2\times 2$ block matrix is

$x_2=2C\left(m\right)+4m^3$ (18)

Inversion formula and computational cost of $n\times n$ block matrix is deduced below.

Theorem 1. Assuming that $A=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{12}^H& A_{22}& A_{23}\\ A_{13}^H& A_{23}^H& A_{33}\end{array}\right]$, $A_{11}$, $A_{22}$ and $A_{33}$ are m order Hermitian matrices, $A_{12}$, $A_{13}$, $A_{32}$ and $A_{33}$ are m order matrix, then

$A^{-1}=\left[\begin{array}{cc}{A}_{11}^{-1}+D{E}^{-1}{D}^H& D{E}^{-1}\\ {\left(D{E}^{-1}\right)}^H& E^{-1}\end{array}\right]$ (19)

where $D=\left(-A_{11}^{-11}{A}_{12},-A_{11}^{-11}{A}_{13}\right)$, $E=\left[\begin{array}{cc}{A}_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}& A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}\\ {\left(A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}\right)}^H& A_{33}-A_{13}^H{A}_{11}^{-1}{A}_{13}\end{array}\right]$, naming Equation (19) as “combination inversion formula”.

Proof: Noting

$\begin{array}{l}\left[\begin{array}{ccc}{I}_m& O_m& O_m\\ {\left(-A_{11}^{-11}{A}_{12}\right)}^H& I_m& O_m\\ {\left(-A_{11}^{-11}{A}_{13}\right)}^H& O_m& I_m\end{array}\right]\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{12}^H& A_{22}& A_{23}\\ A_{13}^H& A_{23}^H& A_{33}\end{array}\right]\left[\begin{array}{ccc}{I}_m& -A_{11}^{-11}{A}_{12}& -A_{11}^{-11}{A}_{13}\\ O_m& I_m& O_m\\ O_m& O_m& I_m\end{array}\right]\\ =\left[\begin{array}{ccc}{A}_{11}& O_m& O_m\\ O_m& A_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}& A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}\\ O_m& {\left(A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}\right)}^H& A_{33}-A_{13}^H{A}_{11}^{-1}{A}_{13}\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& O_{m\times 2m}\\ O_{m\times 2m}& E\end{array}\right]\end{array}$

Then,

$\begin{array}{c}{A}^{-1}=\left[\begin{array}{cc}{I}_m& D\\ O_{2m\times m}& I_{2m}\end{array}\right]{\left[\begin{array}{cc}{A}_{11}& O_{m\times 2m}\\ O_{2m\times m}& E\end{array}\right]}^{-1}\left[\begin{array}{cc}{I}_m& O_{m\times 2m}\\ D^H& I_{2m}\end{array}\right]\\ =\left[\begin{array}{cc}{I}_m& D\\ O_{2m\times m}& I_{2m}\end{array}\right]\left[\begin{array}{cc}{A}_{11}^{-1}& O_{m\times 2m}\\ O_{2m\times m}& E^{-1}\end{array}\right]\left[\begin{array}{cc}{I}_m& O_{m\times 2m}\\ D^H& I_{2m}\end{array}\right]\\ =\left[\begin{array}{cc}{A}_{11}^{-1}+D{E}^{-1}{D}^H& D{E}^{-1}\\ {\left(D{E}^{-1}\right)}^H& E^{-1}\end{array}\right]\end{array}$

where, $E^{-1}$ can be derived by Equation (17). Then inversion computational cost of $3\times 3$ block matrix is

$x_3=x_2+C\left(m\right)+11m^3$ (20)

Theorem 2. Assuming that $A=\left[\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{12}^H& A_{22}& A_{23}& A_{24}\\ A_{13}^H& A_{23}^H& A_{33}& A_{34}\\ A_{14}^H& A_{24}^H& A_{34}^H& A_{44}\end{array}\right]$, $A_{11}$, $A_{22}$, $A_{33}$ and $A_{44}$ are m order Hermitian matrices, $A_{12}$, $A_{13}$, $A_{14}$, $A_{23}$, $A_{24}$ and $A_{34}$ are m order matrix, then

$A^{-1}=\left[\begin{array}{cc}{A}_{11}^{-1}+D{E}^{-1}{D}^H& D{E}^{-1}\\ {\left(D{E}^{-1}\right)}^H& E^{-1}\end{array}\right]$ (21)

where $D=\left(-A_{11}^{-11}{A}_{12},-A_{11}^{-11}{A}_{13},-A_{11}^{-11}{A}_{14}\right)$,

$E=\left[\begin{array}{ccc}{A}_{22}-A_{12}^H{A}_{11}^{-1}{A}_{12}& A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}& A_{24}-A_{12}^H{A}_{11}^{-1}{A}_{14}\\ {\left(A_{23}-A_{12}^H{A}_{11}^{-1}{A}_{13}\right)}^H& A_{33}-A_{13}^H{A}_{11}^{-1}{A}_{13}& A_{34}-A_{13}^H{A}_{11}^{-1}{A}_{14}\\ {\left(A_{24}-A_{12}^H{A}_{11}^{-1}{A}_{14}\right)}^H& {\left(A_{34}-A_{13}^H{A}_{11}^{-1}{A}_{14}\right)}^H& A_{44}-A_{14}^H{A}_{11}^{-1}{A}_{14}\end{array}\right]$

The proof is the same as above, where, $E^{-1}$ can be derived by Equation (19). Then inversion computational cost of $4\times 4$ block matrix is

$x_4=x_3+C\left(m\right)+21m^3$ (22)

Illustration: In order to derive inversion formula of $n\times n$ block matrix, we can solve $A_{11}^{-1}$ first and then construct matrices D and E. Noting that E is $\left(n-1\right)\times \left(n-1\right)$ block matrix, $E^{-1}$ can be solved by applying inversion formula of $\left(n-1\right)\times \left(n-1\right)$ block matrix. By analogy, we can solve $A^{-1}$ by applying “combination inversion formula” finally. This process is basement of deriving inversion computational cost of block matrix.

In summary, we consider the general formula of inversion computational cost of $n\times n$ block matrix. Inversion computational cost of $n\times n$ block matrix is

$x_n=x_{n-1}+C\left(m\right)+\frac{5}{2}\left(n-1\right)m^3+\frac{3}{2}{\left(n-1\right)}^2{m}^3$ (23)

Combining Equations (18), (20) and (22), we can derive the general formula

$x_n=nC\left(m\right)+\frac{1}{2}\left({\left(n-1\right)}^3+4{\left(n-1\right)}^2+3\left(n-1\right)\right)m^3$ (24)

2.3.2. Principles of Algorithm

Assuming that U is segmented equally into K parts, calling $U_i$, where $i=1,\cdots ,K$. The length of each segment is $k=M/K$. Then STCM $R_{stmv}\left(\theta \right)$ is segmented equally into $K\times K$ parts, which order of each segment is k. Considering $R_{stmv}\left(\theta \right)$ is Hermitian matrix, then Equation (12) can be represented as

$R_{stmv}\left(\theta \right)=\left[\begin{array}{cccc}{\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_1^H}& {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_2^H}& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_K^H}\\ {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_2^H}\right)}^H& {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_2^H}& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_K^H}\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_K^H}\right)}^H& {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_K^H}\right)}^H& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_K{U}_K^H}\end{array}\right]$ (25)

Assuming that STCM ${\hat{R}}_{stmv}^{}(n)$ estimated by using the nth snapshot is segmented equally into $K\times K$ parts, calling ${\hat{R}}_{stmv}^{ij}(n)$, where $i,j=1,\mathrm{...},K$. Considering ${\hat{R}}_{stmv}^{}(n)$ is Hermitian matrix, then Equation (13) can be represented as

$\begin{array}{l}\left[\begin{array}{cccc}{\hat{R}}_{stmv}^{11}(n)& {\hat{R}}_{stmv}^{12}(n)& \cdots & {\hat{R}}_{stmv}^{1K}(n)\\ {\left({\hat{R}}_{stmv}^{12}(n)\right)}^H& {\hat{R}}_{stmv}^{22}(n)& \cdots & {\hat{R}}_{stmv}^{2K}(n)\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\hat{R}}_{stmv}^{1K}(n)\right)}^H& {\left({\hat{R}}_{stmv}^{2K}(n)\right)}^H& \cdots & {\hat{R}}_{stmv}^{KK}(n)\end{array}\right]\\ =\alpha \left[\begin{array}{cccc}{\hat{R}}_{stmv}^{11}(n-1)& {\hat{R}}_{stmv}^{12}(n-1)& \cdots & {\hat{R}}_{stmv}^{1K}(n-1)\\ {\left({\hat{R}}_{stmv}^{12}(n-1)\right)}^H& {\hat{R}}_{stmv}^{22}(n-1)& \cdots & {\hat{R}}_{stmv}^{2K}(n-1)\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\hat{R}}_{stmv}^{1K}(n-1)\right)}^H& {\left({\hat{R}}_{stmv}^{2K}(n-1)\right)}^H& \cdots & {\hat{R}}_{stmv}^{KK}(n-1)\end{array}\right]\\ +\left(1-\alpha \right)\left[\begin{array}{cccc}{\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_1^H}& {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_2^H}& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_K^H}\\ {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_2^H}\right)}^H& {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_2^H}& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_K^H}\\ ⋮& ⋮& \ddots & ⋮\\ {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_1{U}_K^H}\right)}^H& {\left({\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_2{U}_K^H}\right)}^H& \cdots & {\displaystyle \underset{i=l}{\overset{h}{\sum }}{U}_K{U}_K^H}\end{array}\right]\end{array}$ (26)

Applying inversion formula of block matrix equation to solving ${\hat{R}}_{stmv}^{-1}(n)$, the basic principle is

1) Solving ${\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}$ by iterative Equations (15) and (16);

2) Constructing matrices

$D=\left(-{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{12}(n),\mathrm{...},-{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{1K}(n)\right)$ (27)

$E=\left[\begin{array}{ccc}{\hat{R}}_{stmv}^{22}(n)-{\left({\hat{R}}_{stmv}^{12}(n)\right)}^H{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{12}(n)& \cdots & {\hat{R}}_{stmv}^{2K}(n)-{\left({\hat{R}}_{stmv}^{12}(n)\right)}^H{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{1K}(n)\\ ⋮& \ddots & ⋮\\ {\left({\hat{R}}_{stmv}^{2K}(n)-{\left({\hat{R}}_{stmv}^{12}(n)\right)}^H{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{1K}(n)\right)}^H& \cdots & {\hat{R}}_{stmv}^{KK}(n)-{\left({\hat{R}}_{stmv}^{1K}(n)\right)}^H{\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}{\hat{R}}_{stmv}^{1K}(n)\end{array}\right]$ (28)

3) Applying inversion formula of block matrix to solving $E^{-1}$ recursively;

4) Applying “combination inversion formula” to solving $R_{stmv}^{-1}(\theta )$.

3. Analysis of Complexity

3.1. Complexity of Classic STMV Algorithm

Assuming that array number is M, the bandwidth of signal is B and CSDM is estimated by the frequency domain output of N snapshots. The computational cost of STMV is mainly focused on two parts as follows:

1) The CSDM is estimated with Equation (5). It should be multiplied $M^2$ times at each frequency point, and $M^{2}NB$ times of multiplication operation is necessary in all;

2) The following calculation is necessary for each beam while estimating spatial spectrum:

a) The STCM is estimated with Equation (6). It should be multiplied $2M^2$ times at each frequency point, and $2M^{2}B$ times of multiplication operation is necessary in all;

b) The inversion operation of STCM needs $M^3$ times of multiplication operation, then it should be multiplied $M^2+M$ times while estimating spatial spectrum with Equation (9).

In summary, assuming the beam number is L, then the total times of multiplication operation is

$N{M}^{2}B+\left[M^3\left(2B+1\right)+M^2+M\right]L$ (29)

3.2. Complexity of Iterative STMV Algorithm

Assuming that array number is M, the bandwidth of signal is B. The computational cost of iterative STMV algorithm is mainly focused on two parts as follows:

1) It should be multiplied by $M^2$ times at each frequency point when utilizing $U={\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)\bar{X}\left(f_k\right)$ to solve U. $M^{2}B$ times of multiplication operation is necessary in all;

2) Each iteration operation needs $3M^2+M$ times of multiplication operation at each frequency point when solving ${\hat{R}}_{stmv}^{-1}(n)$ according to Equations (15) and (16). $\left(3M^2+M\right)B$ times of multiplication operation is necessary in all, then it should be multiplied $M^2+M$ times while estimating spatial spectrum with Equation (9).

In summary, assuming the beam number is L, then the total times of multiplication operation is

$\left[M^2\left(4B+1\right)+M\left(B+1\right)\right]L$ (30)

3.3. Complexity of Block Iterative STMV Algorithm

Assuming that array number is M, the bandwidth of signal is B and U is segmented equally into K parts, calling $U_i$ which length is $k=M/K$, where $i=1,\cdots ,K$. The computational cost of block iterative STMV algorithm is mainly focused on five parts as follows:

1) It should be multiplied by $M^2$ times at each frequency point when utilizing $U={\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)\bar{X}\left(f_k\right)$ to solve U. $M^{2}B$ times of multiplication operation is necessary in all;

2) Each iteration operation needs $3k^2+k$ times of multiplication operation at each frequency point when solving ${\left({\hat{R}}_{stmv}^{11}(n)\right)}^{-1}$ according to Equations (15) and (16). $\left(3k^2+k\right)B$ times of multiplication operation is necessary in all;

3) In order to construct matrices D and E, it needs to solve ${\hat{R}}_{stmv}^{12}(n),\mathrm{...},{\hat{R}}_{stmv}^{1K}(n)$, ${\hat{R}}_{stmv}^{22}(n),\mathrm{...},{\hat{R}}_{stmv}^{2K}(n)$, ${\hat{R}}_{stmv}^{33}(n),\mathrm{...},{\hat{R}}_{stmv}^{3K}(n)$, …, ${\hat{R}}_{stmv}^{KK}(n)$. $\left(\frac{K\left(K+1\right)}{2}-1\right)k^{2}B$ times of multiplication operation is necessary in all;

4) Applying inversion formula of block matrix to solving $E^{-1}$ recursively, It should be multiplied by $\left(K-1\right)C\left(k\right)+\frac{1}{2}\left({\left(K-2\right)}^3+4{\left(K-2\right)}^2+3\left(K-2\right)\right)k^3$ times;

5) Applying “combination inversion formula” to solving $R_{stmv}^{-1}(\theta )$, It should be multiplied by ${\left(K-1\right)}^2{k}^3+\left(K-1\right)k^3$ times. Then it should be multiplied $M^2+M$ times while estimating spatial spectrum with Equation (9).

In summary, assuming the beam number is L, then the total times of multiplication operation is

$\begin{array}{l}\left(M^{2}B+\left(3k^2+k\right)B+\left(\frac{K\left(K+1\right)}{2}-1\right)k^{2}B+\left(K-1\right)C\left(k\right)\right)L\\ +\left(\frac{1}{2}\left({\left(K-2\right)}^3+4{\left(K-2\right)}^2+3\left(K-2\right)\right)k^3+{\left(K-1\right)}^2{k}^3+\left(K-1\right)k^3+M^2+M\right)L\end{array}$ (31)

where, the computational cost of k order inverse matrix is $C\left(k\right)$, where $C\left(k\right)=k^3$. Considering $k=M/K$ meanwhile, Equation (31) can be simplified to

$\left(\left(\frac{3}{2}+\frac{1}{2K}+\frac{2}{K^2}\right)M^{2}B+\frac{1}{K}MB+\left(\frac{1}{2}-\frac{1}{2K^2}\right)M^3+M^2+M\right)L$ (32)

Illustration: In fact, considering ${\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)$ is diagonal matrix, the computational cost can be reduced to MB when utilizing $U={\stackrel{\dot{}}{T}}_p\left(f_k,\theta \right)\bar{X}\left(f_k\right)$ to solve U. Now, the total times of multiplication operation

$\left(\left(\frac{1}{2}+\frac{1}{2K}+\frac{2}{K^2}\right)M^{2}B+\left(1+\frac{1}{K}\right)MB+\left(\frac{1}{2}-\frac{1}{2K^2}\right)M^3+M^2+M\right)L$ (33)

Assuming that the bandwidth of signal $B=5\text{kHz}$ and a beam is formed every two degrees. That’s to say, it should estimate spatial spectrum at $L=91$ beams. The snapshot is $N=4$, $K=20$. Then, the ratio of computational cost between the STMV algorithm and block iterative algorithm for various array number M are shown in Figure 1.

The computational cost of iterative algorithm proposed in paper [8] is approximately $2/M$ times of computational cost of STMV algorithm when Equation (30) is divided by Equation (29). For the same of analysis, Equations (29),

(33) can be simplified to $2M^{3}BL$, $\frac{\left(M^{2}B+M^3\right)L}{2}+\frac{M^{2}BL}{2K}$. If we take $L\gg N$

and $B\gg M$ both into consideration, then the ratio of computational cost of two algorithms is

$\frac{\frac{\left(M^{2}B+M^3\right)L}{2}+\frac{M^{2}BL}{2K}}{2M^{3}BL}=\frac{1+\frac{M}{B}+\frac{1}{K}}{4M}\approx \frac{1}{4M}$ (34)

Figure 1. Ratio of computational cost of two algorithms.

The computational cost of the algorithm proposed in this paper is $1/4M$ times of the computational cost of original algorithm, and 1/8 times of the computational cost of iterative algorithm. The computational cost drops dramatically when $M\gg 5$.

4. Simulation and Analysis

4.1. Azimuth Resolution

The simulation conditions are listed as follows. Two sources signals and background noise defined over the same frequency band 0.05 fs - 0.25 fs are Gaussian white noise which are independent with each other. The signal-to-noise ratio (SNR) of each source is 0 dB, sampling frequency is fs and sound speed is 1500 m/s. The receiver is a uniformly spaced linear array of 48 sensors with inter-element spacing of d, corresponding to the maximum frequency. The length of each snapshot is 0. 8 s and the integration time of beam-former is 3.2 s, where $\alpha =2/3$. The bearing angles of two sources are assumed in two conditions: 80˚ and 88˚, 80˚ and 84˚. The spatial spectrum estimated by STMV algorithm and block iterative STMV algorithm ( $K=6$ ) are shown in Figure 2. In the condition of their bearing angles lied at 80˚ and 84˚, both of the above algorithms are simulated independently for 5 times and the results are shown in Figure 3.

It is seen from the simulation results that compared to STMV algorithm, there is no significant change both on azimuth resolution and main lobe width for block iterative STMV algorithm. The only change is that the side lobe fluctuate slightly. In summary, the azimuth resolution of block iterative STMV algorithm has no significant decrement while the computational cost reduced sharply.

4.2. Ability of Detecting Weak Source

The simulation conditions are listed as follows. There is an interference lied at the bearing angle 100˚ with SNR = 5 dB, and weak target source lied at the bearing angle 80˚. The simulation is in two conditions: SNR = −15 dB and SNR = −20 dB, other parameters are the same as former. The spatial spectrum estimated by STMV algorithm and block iterative STMV algorithm are shown in Figure 4. In the condition of SNR = −20 dB, both of the above algorithms are simulated independently for 5 times and the results are shown in Figure 5.

(a) (b)

Figure 2. Comparison of two algorithms on spatial spectrum. (a) Two sources lie at 80˚ and 88˚; (b) Two sources lie at 80˚ and 84˚.

(a) (b)

Figure 3. Comparison of two algorithms on spatial spectrum. (a) STMV algorithm; (b) Block iterative STMV algorithm.

(a) (b)

Figure 4. Comparison of two algorithms on spatial spectrum in different SNR. (a) SNR = 5 dB and −15 dB; (b) SNR = 4 dB and −20 dB.

(a) (b)

Figure 5. Comparison of two algorithms on spatial spectrum. (a) STMV algorithm; (b) Block iterative STMV algorithm.

It is seen from the simulation results that block iterative STMV algorithm inherits the ability to detect weak target source of STMV algorithm while decreasing the computational cost dramatically. Moreover, because the side lobe of spatial spectrum is stationary, the weak target source can be detected clearly. This algorithm can realize detection of weak target source in strong coherent interference.

5. The Sea Trial Data Analysis

The sea trial data is analyzed to verify the adaptation and stability of the block iterative STMV algorithm. The trial frame is shown in Figure 6, a uniformly spaced linear array of 48 sensors is towed behind the towboat sailing straightly, whose original course is 0˚. The target 2 originally lies at about bearing angle 75˚ and sails in the direction of 0˚. The target 4 sails in the direction of array. During the experiment, shipping is busy and compared to targets 2, 4, targets 1, 3 are far away from the towboat whose radiated noise is weak.

The spatial spectrum estimated with STMV algorithm and the block iterative STMV algorithm after the data filtered with a differential whiten filter is shown in Figure 7, Figure 8.

It is not hard to see that there are five targets from Figure 7, where the target at bearing angle about 20˚ is towboat. For the reason that target (i.e. target 3) at bearing angle about 100˚ is far away from the towboat, its bearing changes slowly. Target (i.e. target 1) at bearing angle about 80˚ whose bearing change fast has bearing crossing with target 3 at about 40 s, but these two targets can be discriminated during the short time before or after the crossing moment, which proves further the high azimuth resolution character of block iterative STMV algorithm. Besides, from Figure 8, we can see that there are two weak targets at bearing angle about 65˚, whose bearings are first close, then far away.

In order to further illustrate the ability to detect weak target, the spatial spectrum of the 27^th, the 45^th and the 126^st moment are shown as Figures 9-11. Although the SNR of target 2, and target 4 is about 20 dB lower than target 1 and target 3, these two weak targets still can be detected validly.

Figure 6. Trail frame.

Figure 7. Diagram of STMV.

Figure 8. Diagram of block iterative STMV.

Figure 9. Spatial spectrum at 27 s.

Figure 10. Spatial spectrum at 45 s.

Figure 11. Spatial spectrum at 126 s.

Multi-targets with strong and weak SNR can be detected with block iterative STMV algorithm in Figures 9-11. The directional index within 1 - 2 dB for targets and higher azimuth discrimination are improved than STMV, and the power on background azimuth without target is suppressed about 1dB. It’s due to the matrixes are added first in Equation (5) for STMV algorithm and the spectrum are added first in Equation (10) for block iterative STMV algorithm for the continuum sequence, the correlation of target’s radial signal is utilized for block iterative STMV algorithm and better performance is achieved. Both of the target signal background noise are Gaussian noise which without any correlation in simulation, so the same results can’t be obtained in Figure 4 and Figure 5.

6. Conclusion

First, the computational cost of STMV algorithm and iterative STMV algorithm is analyzed and their application in practical engineering is restricted due to huge computational cost. This paper proposed block iterative STMV algorithm based on one-phase regressive filter, matrix inversion lemma and inversion of block matrix and analyzed its computational cost, which is 1/4 M times of computational cost of STMV algorithm. This algorithm maintains high azimuth resolution and advantage of weak target detection of STMV algorithm showed by simulation results. The computational cost is reduced by about 140 times if this algorithm is used to analyze the sea trail data of a uniformly spaced linear array of 48 sensors. Meanwhile this algorithm can validly detect weak targets when several targets with very large power difference and approximate bearing angles presented. The performance of block iterative STMV algorithm which has the directional index within 1 - 2 dB for targets, the higher azimuth discrimination and the background azimuth power suppress about 1 dB is improved than STMV algorithm. It also has good robustness when processing sea trial data, which lays the foundation of its engineering application.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Van Trees, L. (2003) Optimum Array Processing. John Wiley & Sons Inc., New York, 48.
[2]	Wang, Y.L. (2004) Theory and Algorithm on Spatial Spectrum Estimation. Tsinghua University Press, Beijing. (In Chinese)
[3]	Swingler, D. (1999) A Low-Complexity MVDR Beamformer for Use with Short Observation Times. Proceeding of the IEEE, 47, 1154-1160. https://doi.org/10.1109/78.752616
[4]	Zhou, S.Z. and Du, X.M. (2009) Research and Application of Fast-Convergent Minimum Variance Distortionless Response Algorithm. Chinese Journal of Acoustics, 34, 515-520. (In Chinese)
[5]	Lu, Z.X., Zhou, S.Z. and Gao, Y. (2014) Research on SAVSTMV Beamforming Algorithm. Chinese Journal of Acoustics, 33, 477-480. (In Chinese)
[6]	Xu, G. and Zhou, S.Z. (2014) Development and Application of MVDR Adaptive Beamforming Technique in Underwater Acoustic. Chinese Journal of Acoustics, 33, 554-558. (In Chinese)
[7]	Qian, B., Wang, P. and Wang, X. (2013) Subarray Optimization and Partition Technology in Dimension Reduction Steered Minimum Variance. Journal of Harbin Engineering University, 34, 445-449. (In Chinese)
[8]	Zhu, D.Z., et al. (2010) Iterative Algorithm of Steered Minimum Variance and Its Application in Weak Targets Detection. Journal of Shanghai Jiaotong University (Science), 15, 694-701. https://doi.org/10.1007/s12204-010-1071-6
[9]	Zhang, X.D. (2004) The Analysis and Application of Matrix. Tsinghua University Press, Beijing. (In Chinese)