## 2 Extracting Green's Function  
For the array-invariant method, extracting the Green's function with blind deconvolution [8; 9; 11] from the received signals is necessary. Assume that an unknown source signal is \(s(t)\). The frequency spectrum of \(s(t)\) is \(S(\omega)=|S(\omega)|e^{i\Phi_{s}(\omega)}\), where \(\omega\) is the radian frequency and \(\Phi_{s}(\omega)\), is the corresponding phase. In this study, we focus on using a horizontal array in shallow water to realize passive localization. Consider a horizontal array with \(N\) elements, as shown in Figure 1. The received signal of the \(n\)th element can be expressed as  
\[P_{n}(\omega)=G(\mathbf{r}_{n},\mathbf{r}_{s},\omega)|S(\omega)|e^{i\Phi_{s}( \omega)}, \tag{1}\]  
where \(\mathbf{r}_{n}\) and \(\mathbf{r}_{s}\) are the position vectors of the \(n\)th element of the array and the unknown source, respectively, and \(G(\mathbf{r}_{n},\mathbf{r}_{s},\omega)\) is the Green's function between the source at \(\mathbf{r}_{s}\) and the \(n\)th element at \(\mathbf{r}_{n}\). In a ray path approximation [11], the Green's function \(G(\mathbf{r}_{n},\mathbf{r}_{s},\omega)\) can be written as  
\[G(\mathbf{r}_{n},\mathbf{r}_{s},\omega)=\sum_{k=1}^{K}A_{k}e^{i\omega\{\tau_{ n}(\theta_{k},\phi)-T(\theta_{k},\phi)\}}, \tag{2}\]  
where \(\phi\) is the horizontal azimuth shown in Figure 1, \(\theta_{k}\) is the grazing angle of the \(k\)th ray, \(\tau_{n}(\theta_{k},\phi)\) is the local time delay of the \(k\)th ray at the \(n\)th element, \(T(\theta_{k},\phi)\) is travel times of the \(k\)th ray, and \(A_{k}\) is the amplitude of the \(k\)th ray.  
Figure 1: Coordinate system defined with the azimuth angle \(\phi\), beam angle \(\widetilde{\theta}\), and the grazing angle \(\theta\) of incoming ray arrival [10].  
For the horizontal array in Figure 1, the local time delays depend mainly on the horizontal azimuth due to a small grazing angle, thus \(\tau_{n}(\theta_{k},\phi)\approx\tau_{n}(\phi)\)[11]. Subsequently, by steering the beam to \(\phi\), the output of the conventional plane-wave beamforming is identified by  
\[\begin{array}{rcl}F(\omega;\theta_{k},\phi)&=&\sum\limits_{n=1}^{N}e^{-i\omega \tau_{n}(\theta_{k},\phi)}P_{n}(\omega)\\ &=&|S(\omega)|e^{i\Phi_{k}(\omega)}\sum\limits_{n=1}^{N}\sum\limits_{k_{n}=1} ^{K}A_{k_{n}}e^{i\omega\left(\tau_{n}(\theta_{k_{n}},\phi)-T_{k_{n}}\right)}e^{- i\omega\tau_{n}(\theta_{k},\phi)}\\ &\approx&|F(\omega;\theta_{k},\phi)|e^{i\left(\Phi_{k}(\omega)-i\omega\tau_{k }\right)}\end{array} \tag{3}\]  
Blind convolution uses the beamforming output phase, \(\psi(\omega,\phi)=\Phi_{s}(\omega)-\omega T(\theta_{k},\phi)\), and removes the phase component \(\Phi_{s}(\omega)\) of the unknown source signal from \(P_{n}(\omega)\) by a phase rotation, expressed as:  
\[\hat{G}(\mathbf{r}_{n},\mathbf{r}_{s},\omega)=\frac{P_{n}(\omega)}{\sqrt{\sum \limits_{n=1}^{N}|P_{n}(\omega)|}^{2}}e^{-i\psi(\omega,\phi)}=\frac{G(\mathbf{ r}_{n},\mathbf{r}_{s},\omega)}{\sqrt{\sum\limits_{n=1}^{N}|G(\mathbf{r}_{n}, \mathbf{r}_{s},\omega)|}^{2}}e^{i\omega T(\theta_{k},\phi)} \tag{4}\]  
The beam-time migration in the existing array-invariant method [11] can be obtained by processing the extracted \(\hat{G}(\mathbf{r}_{n},\mathbf{r}_{s},\omega)\) with conventional plane-wave beamforming.