## 3 2D Deconvolution for Obtaining Better Beam-Time Migration
### 2D Deconvolution  
The local time delays \(\tau_{n}(\theta_{k},\phi)\) in Equation (2) can be written as \(\tau_{n}(\sin\widetilde{\theta}_{k})\), which is the function of the beam angle measured \(\sin\widetilde{\theta}_{k}\), and \(T(\theta_{k},\phi)\) can be abbreviated as \(T_{k}\). Then, the expression of the Green's function in (2) can be rewritten as  
\[G(\mathbf{r}_{n},\mathbf{r}_{s},\omega)=\sum\limits_{k=1}^{K}A_{k}e^{i\omega \left(\tau_{n}(\sin\widetilde{\theta}_{k})-T_{k}\right)} \tag{7}\]Let \(\mathbf{G}(\mathbf{r}_{s},\omega)\!=[G(\mathbf{r}_{1},\mathbf{r}_{s},\omega),G( \mathbf{r}_{2},\mathbf{r}_{s},\omega),\cdots G(\mathbf{r}_{N},\mathbf{r}_{s}, \omega)]\). We define the 2D rotation matrix as  
\[\mathbf{\Psi}=[\mu(t_{0}),\mu(t_{1}),\ldots\mu(t_{M})], \tag{8}\]  
where \(\mu(t_{m})\!=[\mu_{1}(t_{m}),\mu_{2}(t_{m}),\cdots\mu_{N}(t_{m})]^{T}\), and \(\mu_{n}(t_{m})=(1/N)e^{j\omega(\tau_{n}(\sin\varphi)-t_{m})}\). Assuming that we process \(M\) samples, \(t_{0}\) and \(t_{M}\) are the start and end time of the sampling segment, respectively. Considering (7), after the time-beam rotation, the output can be written as  
\[ \begin{array}{l}Y(\sin\varphi,t)=\sum\limits_{\omega_{\min}^{ \omega_{\max}}}^{\omega_{\max}}\mathbf{G}(\mathbf{r}_{s},\omega_{q})\mathbf{ \Psi}=\sum\limits_{\omega_{\min}^{\omega_{\max}}}^{\omega_{\max}}\sum\limits _{k=1}^{K}\sum\limits_{n=0}^{N}e^{-j\omega_{q}((\tau_{n}(\sin\varphi)-\tau_{ n}(\sin\tilde{\theta}_{k}))+(T_{k}-t))}\\ =\sum\limits_{k=1}^{K}\frac{A_{k}}{N}\omega_{\omega}^{\omega_{\max}}\sum\limits _{n=1}^{N}e^{-j\omega_{q}((\frac{nt(\sin\varphi-\sin\tilde{\theta}_{k})}{c}+( T_{k}-t))}\end{array}  \tag{9}\]  
where \(\omega_{\min}\) and \(\omega_{\max}\) are the highest and lowest frequencies of the Green's function \(\mathbf{G}(\mathbf{r}_{s},\omega_{q})\) in the discrete frequency domain, respectively. The power of \(Y(\sin\varphi,t)\) is  
\[B(\sin\varphi,t)=\left|Y(\sin\varphi,t)\right|^{2}=\sum\limits_{k=1}^{K}\frac {\left|A_{k}\right|^{2}}{N^{2}}\bigg{|}\sum\limits_{\omega_{q}=\omega_{\min} }^{\omega_{\max}^{\omega_{\max}^{\omega_{\max}^{\omega_{\max}^{\omega_{\max} }}}}}\sum\limits_{n=1}^{N}e^{-j\omega_{q}((\frac{nt(\sin\varphi-\sin\tilde{ \theta}_{k})}{c}+(T_{k}-t))})\bigg{|}^{2}. \tag{10}\]  
After a simple derivation, Equation (10) can be written as a 2D convolution of the beam-time pattern, denoted as PSF, with the source power distribution \(S\). \(PSF\) depends only on the difference in bearing and time, so \(PSF\) is shift-invariant in the steering angle and time, and can be written as  
\[B(\sin\varphi,t)=\int_{t_{0}}^{t_{M}}\int_{-1}^{1}PSF((\sin\varphi-\sin\hat{ \theta}),(t-\hat{t}))\cdot S(\sin\hat{\theta},\hat{t})d\sin\hat{\theta}d\hat{t}, \tag{11}\]  
where  
\[PSF((\sin\varphi-\sin\hat{\theta}),(t-\hat{t}))=\left|\frac{1}{N}\sum\limits_ {\omega_{q}=\omega_{\min}}^{\omega_{\max}}\sum\limits_{n=0}^{N}e^{-j\omega_{q }((\frac{nt(\sin\varphi-\sin\hat{\theta})}{c}+(t-\hat{t}))}\right|^{2}, \tag{12}\]  
and  
\[S(\sin\hat{\theta},\hat{t})=\sum\limits_{k=1}^{K}\left|A_{k}\right|^{2}\!\delta ((\sin\hat{\theta}-\sin\widetilde{\theta}_{k})+(\hat{t}-T_{k})). \tag{13}\]  
It should be noted that \(S(\sin\hat{\theta},\hat{t})\) separates the multiple arrivals, in terms of both the bearing and time domains, much better than \(B(\sin\varphi,t)\) does. If we estimate \(S(\sin\hat{\theta},\hat{t})\) from \(B(\sin\varphi,t)\) using 2D deconvolution, high-resolution beam-time migration can be obtained, which helps improve the localization accuracy. Some well-known deconvolution algorithms, such as the Richardson-Lucy (R-L) [12; 13; 14; 15; 16] algorithm, have been widely applied in various fields. We propose the use of the R-L algorithm to deconvolve \(B(\sin\varphi,t)\) both the bearing and time domains to improve the localization accuracy of the array-invariant method.