## 2 Methods
The phase relation along the antidiagonal only depends on the channel position \(qw\), indicating that the direct waves from the optical absorbers manifest themselves as a particular coherence on the antidiagonals of the matrix. Using the coherent property, we can extract the direct wave by the following steps. Firstly, the original matrix **K** is rotated 45\({}^{\circ}\) counterclockwise so that the coherent property occurs on the columns. The elements are divided into two parts based on their symmetry characteristic and are expanded to two new matrices: \(\textbf{A}_{1}\) and \(\textbf{A}_{2}\)[26]. The matrices \(\textbf{A}_{1}\) and \(\textbf{A}_{2}\) are called indifferently **A** as they are filtered in the same way. A filtering matrix \(\textbf{F}=\textbf{C}\textbf{C}^{\text{t}}\) is then applied to extract the coherent part from matrix **A** as: \(\textbf{A}^{\text{F}}=\textbf{FA}=\textbf{FA}^{\text{C}}+\textbf{FA}^{\text{R}}\), where \(\textbf{C}=[\exp(jky_{1}^{2}/2Z),\)\(\ldots,\)\(\exp(jky_{m}^{2}/2Z),\)\(\ldots,\)\(\exp(jky_{N_{n}}^{2}/2Z)]\) with \(u=1,\)\(2,\)\(\ldots,\)\(N_{A}\), and \(y_{u}=[x_{m}+x_{n}+(N_{A}-1)w]\) is the characteristic spacing of direct waves. The coherent part remains unchanged during the filtering process, as \(\textbf{FA}^{\text{C}}=\textbf{A}^{\text{C}}\) and the random part decreases by a factor \(\sqrt{(\text{N}\,+\,1)/2}\). Finally, a filtered matrix \(\textbf{K}^{\text{F}}\) with the original coordinate can be obtained by applying the inverse process of the first step to matrices \(\textbf{A}_{1}\)F and \(\textbf{A}_{2}\)F.  
Using a time reversal operator, we can recover the image of the optical absorbers from the matrix \(\textbf{K}^{\text{F}}\) corresponding to the single sector area [28; 29]. The singular value decomposition is applied to the matrix \(\textbf{K}^{\text{F}}=\textbf{U}^{\text{F}}\textbf{A}^{\text{F}}\textbf{V}^{ \text{F}\dagger}\), and each optical absorber is mainly associated with one non-zero singular \(\lambda_{i}^{\text{F}}\) contained in diagonal matrix \(\boldsymbol{\Lambda}\)[30]. Therefore, by numerically back-propagating the singular vector, the image of the \(i\)th singular valued at a depth of \(Z=cT\) can be achieved by the time reversal operator \(\textbf{I}_{i}=\lambda_{1}^{\text{F}}|\textbf{G}\times\textbf{V}|\), where **G** with the component \(g_{m}\)\(=\exp(jkr_{m})/\sqrt{r_{m}}\) is the propagating operator between the array and the focal plane in a homogeneous media, and \(r_{m}\) is the distance between the \(m\)th element in the array and the point in the focal plane. By superimposing the imaging results of multiple sectors, we can further improve the imaging quality and reduce the influence of the inhomogeneous aberrating effect of the scattering layer. For the wanted imaging point A, we calculated the image value \(I_{m}\)\({}^{\text{A}}\) for all the linear arrays with \(m=1,\)\(2,\)\(\ldots,\)\(N_{S}\). Taking the first sector area as an example, the image value for points corresponding to the array can be obtained by applying the time reversal operator. If point A is within the imaging area, we can find two imaging points, B and C, closest to point A in the focal plane and then perform a weighted evaluation based on the distance from point A to points B and C:  
\[I_{1}^{\rm A}=\frac{I_{1}^{\rm B}d_{\rm AC}+I_{1}^{\rm C}d_{\rm AB}}{d_{\rm AC }+d_{\rm AB}} \tag{3}\]  
where \(d_{\rm AC}\) and \(d_{\rm AB}\) represent the distances between A, B and A, C, respectively. If point A is outside the imaging area, the imaging value is then recorded as 0. When we calculate the imaging value of point A for all sector areas, its final value is the result of taking the summation and average to all the non-zero values: \(I^{\rm A}=\sum_{m=1}^{N_{\rm s}}I_{m}^{\rm A}/(N_{s}-N_{z})\), where \(N_{\rm Z}\) is the number of zero values. As such, we can fully use the valid information received by all the transducers of the annular array for image reconstruction of each point in the ROI.