**4. Numerical Experiment**

The method is tested on the 3D Overthrust model. The model is 16 km, 5 km, and 3 km in size, and the migration velocity is obtained by Gaussian smoothing with a window of 40 m. The Ricker wavelet with the main frequency of 30 Hz is used in this experiment. As shown in Figure 3, the RTM result has an insufficient spatial resolution, and the amplitude of the reflectors attenuates rapidly as the depth increases.

**Figure 3.** The 3D data volume of RTM in the 3D Overthrust model.

Then we compute the PSF using the same wavelet. Based on the migration velocity, the sampling interval (400 m, 400 m, 400 m) of the global PSF is set for the scattering model. To give intuitive insight into PSF, we present the X-Z tangent plane for point spread functions in Figure 4. As the geological structure and the illumination condition change, we can easily find that the difference of global PSFs in morphology is not negligible.

Figure 5 shows the results from the proposed imaging-domain LSM. Compared with the RTM results, the LSM results show higher resolution, fewer migration noise, more balanced amplitude in deep and better focusing on diffracted structures (as shown in Figure 6). To compare the spatial resolution between these results, we compute the vertical wavenumber spectrum of Figures 3 and 5 for the whole model and present them in Figure 7. From the vertical wavenumber spectrum, we can find that the result from the LSM method has a broader frequency-band, especially for mid-high wave numbers. It means that the

proposed LSM method can estimate the Hessian operator and remove those Hessian effects on the RTM results, which thus provides a more accurate reflectivity model.

**Figure 4.** The 2D Display of the X-Z tangent plane of the point spread function.

**Figure 5.** *Cont.*

**Figure 5.** The 3D data volume of LSM in the 3D Overthrust model.

**Figure 6.** The 2D display of the X-Z tangent planes results from (**a**) RTM; (**b**) LSM.

**Figure 7.** The comparison of vertical-wave-number spectrum between RTM (red) and LSM (blue) from Figures 3 and 5. The unit of y-axis is decibel, and x-axis indicates the wavenumber component.

### **5. Field Case**

The target area is a part of the Tarim oil field in Northwest China, which is located in the middle of the North Tarim Uplift Belt to the east of the Halahatang depression. It has a large span, the target layer is deeply buried, the physical properties of the matrix are poor, and the type and scale of fracture caves are seriously affected by karst weathering. From the current development stage, the positioning accuracy of small and medium-sized karst cave cannot meet the requirements of resolution analysis of venting leakage around the well. In this work area, the geophone array recorded 33,000 shots at 5 m depth covering an area of 81 km2. The main conundrum to be tackled in this work area is high-resolution imaging of Ordovician carbonate fracture-cavity reservoirs, with the depth ranging from 5.5 to 8 km. The migration velocity is built by a 3D travel time tomography as shown in Figure 8. The maximum cut-off frequency of the wavelet is 70 Hz.

**Figure 8.** Migration velocity for Tarim oil field.

Figure 9 shows the results from RTM and LSM at the inline index of 200. Compared with the RTM results (Figure 9a), the LSM results (Figure 9b) provide higher spatial resolution, finer reflection events, more balanced amplitude, better focus ability, and more accurate fractured-cavity imaging. Figure 10 shows the zoomed region in Figure 9. We can find that the LSM result provides finer descriptions of fractured-cavity imaging and the boundaries between adjacent fractured-cavities, and the fractured-cavity bodies with small sizes in Ordovician carbonatite can be presented with higher resolution. Moreover, the deep fracture in the LSM result (Figure 10b) is better retrieved, which demonstrates that the proposed method is effective in fracture description. Figure 11 shows the vertical wavenumber spectrum corresponding to the imaging results. The wave-number distribution of the LSM result is much broader, which lifts up at the low and mid-high sides of wave numbers. Besides, LSM can clearly provide information on cavities in the fracture zones, provide a high-quality scientific basis to analyze the development of fracture-controlled fractured-cavities and give guidance to the interpretation of reservoirs. To give an intuitive sight of the caves in Figure 10 denoted by red arrows, we compute the amplitude attributes from the RTM and LSM results and present a regional transverse section at the depth of these caves in Figure 12. We can find that the LSM results can distinguish the adjacent caves and provide better focusing, highlighting the proposed method can provide a high resolution in the horizontal direction and give a great help for reservoir interpretation and development. However, due to the limited observation aperture and shallow desert area, conventional RTM cannot provide correct images. The logging result located at the arrows in Figure 12 can also support the migration results, highlighting the effectiveness of the proposed method.

(**a**) (**b**)

**Figure 9.** Comparison between migration results for CDP index of 1198. (**a**) results from RTM; (**b**) results from LSM.

**Figure 10.** Locally zoomed-in reservoir sections marked in Figure 9. (**a**) results from RTM; (**b**) results from LSM.

**Figure 11.** Comparison in vertical-wave-number spectrum (RTM: red line; LSM: blue line). The unit of y-axis is decibel, and x-axis indicates the wavenumber component.

**Figure 12.** Comparison between results of amplitude attributes. (**a**) zoomed results from RTM; (**b**) zoomed results from LSM. The target caves are corresponding to the ones marked by red arrows in Figure 10.

### **6. Conclusions**

Utilization of the point spread function and the global space-varying high-dimension deconvolution algorithm effectively solves the demand of excessive computation in the iterative data-fitting process of classical least-squares migration and alleviates the high dependency on the accuracy of seismic wavelet and velocity model. Both the theoretical analysis and the model experiment demonstrate that the global space-varying point spread function can effectively approximate the Hessian and the high-dimensional spatial deconvolution algorithm can replace the deburring effects from the Hessian operator. The application on a field 3D data suggests that the method is effective in improving the imaging resolution, especially for caves and fractures, and provide higher-quality seismic images for ultra-deep fractured-cavity reservoirs.

**Author Contributions:** Conceptualization, B.L.; methodology, M.S., Y.B.; software, M.S., C.X. and B.L.; validation, M.S., C.X.; formal analysis, B.L., Y.B.; investigation, M.S., B.L.; resources, M.S., C.X. and B.L.; writing—original draft preparation, M.S.; writing—review and editing, M.S., C.X. and B.L.; supervision, B.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Key Research and Development Program of China (No. 2017YFB0202900), and China Postdoctoral Science Foundation (Grant No. 2020M671539).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.
