## 2 Methods  
As shown in Figure 1, an imaging target is placed in the region of interest (ROI), which is surrounded by a skull-like scattering layer. Under the illumination of the pulse laser, the optical absorber will radiate ultrasound waves to the surrounding medium as a result of the PA effect. After the acoustic waves propagate through the random scattering layer, a passive ultrasound annular array with 576 elements detects the ultrasound signal. The annular array is then divided into \(N_{\text{S}}\) sectors for data processing and imaging. Each sector has \(N=45\) elements and can be approximated as a linear array so that the filtered imaging method proposed in our previous studies can be applied for the image reconstruction [26; 27]. For each sector array, the recorded \(N\)-channel signals are arranged in a vector \(\mathbf{H}(t)=[H_{1}(t),H_{2}(t),\,\ldots,H_{n}(t),\,\ldots,H_{N}(t)]\) where \(n=1,2,\,\ldots\,N\). Figure 1: The schematic of the scenario considered in this study. For any target placed in the region of interest (ROI) within a scattering layer, its propagation distances to the \(n\)th transducers of a sector array, and the equivalent linear array are \(l_{1}\) and \(l_{2}\), respectively. The corresponding delay is then given by \(d_{n}\) = \((l_{1}-l_{2})F_{\text{S}}/c\), where \(c\) is the sound speed and \(F_{\text{S}}\) is the sampling frequency. The sector array can be approximated to a linear array by calculating the corresponding delay. For any imaging target, its propagation distances to the \(n\)th transducers of the two arrays are \(l_{1}\) and \(l_{2}\), respectively. The corresponding delay is \(d_{n}=(l_{1}-l_{2})F_{s}/c\), where \(c\) is the sound speed and \(F_{s}\) is the sampling frequency. The time domain signals for the linear array are then \(\textbf{J}(t)=[H_{1}(t-d_{1}),H_{2}(t-d_{2}),\)\(\ldots\)\(H_{n}(t-d_{n}),\)\(\ldots\)\(H_{N}(t-d_{N})]\). The signals can be written in their frequency domain by applying a Fourier transform to \(\textbf{P}(Tf)=[P_{1}(Tf),P_{2}(Tf),\)\(\ldots\)\(P_{n}(Tf),\)\(\ldots\)\(P_{N}(Tf)]\), where \(P_{n}(Tf)\) is the short-time Fourier transform of the \(n\)th channel signal in the time window \([T-\Delta t/2,\)\(T+\Delta t/2]\), \(T\) is the time of flight, and \(\Delta t\) is the width of the window. \(\textbf{P}(Tf)\) contains both the direct wave \(\textbf{P}^{\text{D}}(Tf)\) and the scattering wave \(\textbf{P}^{\text{S}}(Tf)\). When the concentration of the scatterers is sufficiently high, the scattering wave occupies a dominant position and the imaging quality is greatly degraded due to the randomness of the scattering process. To reduce the influence of the randomly scattering, a matrix filtering method can be applied to extract the direct wave from the received signals of a linear array [26]. Firstly, a response matrix can be constructed as:  
\[\textbf{K}=\textbf{P}\textbf{P}^{\text{T}}=\underbrace{\textbf{P}^{\text{D}}( \textbf{P}^{\text{D}})^{\text{T}}}_{\text{Coherence},\textbf{K}^{\text{C}}}+ \underbrace{[\textbf{P}^{\text{D}}(\textbf{P}^{\text{S}})^{\text{T}}+\textbf{ P}^{\text{S}}(\textbf{P}^{\text{D}})^{\text{T}}+\textbf{P}^{\text{S}}( \textbf{P}^{\text{S}})^{\text{T}}]}_{\text{Random},\textbf{K}^{\text{R}}} \tag{1}\]  
Depending on whether or not contain the term \(\textbf{P}^{\text{S}}\), the matrix **K** can be divided into two parts: \(\textbf{K}^{\text{R}}\) and \(\textbf{K}^{\text{C}}\), where \(\textbf{K}^{\text{R}}\) is a random term due to the random nature of the scatterer distribution and the scattering paths and \(\textbf{K}^{\text{C}}\) is independent of the scattering paths and shows a deterministic relation of the phase between its elements along the antidiagonal:  
\[\beta_{q}=\frac{k_{m-q,m+q}^{\text{C}}}{k_{m,m}^{\text{C}}}=\exp[jk(qw)^{2}/Z] \tag{2}\]  
where \(k\) is the wave number, \(w\) is the pitch size, and \(Z\) is the depth of the time window.