Exploring the subsurface structure and stratification of Mars advances our understanding of Martian geology, hydrological evolution and palaeoclimatic changes, and has been a main task for past and continuing Mars exploration missions^1–10. Utopia Planitia, the smooth plains of volcanic and sedimentary strata that infilled the Utopia impact crater, has been a prime target for such exploration as it is inferred to have hosted an ancient ocean on Mars^11–13. However, 45 years have passed since Viking-2 provided ground-based detection results. Here we report an in situ ground-penetrating radar survey of Martian subsurface structure in a southern marginal area of Utopia Planitia conducted by the Zhurong rover of the Tianwen-1 mission. A detailed subsurface image profile is constructed along the roughly 1,171 m traverse of the rover, showing an approximately 70-m-thick, multi-layered structure below a less than 10-m-thick regolith. Although alternative models deserve further scrutiny, the new radar image suggests the occurrence of episodic hydraulic flooding sedimentation that is interpreted to represent the basin infilling of Utopia Planitia during the Late Hesperian to Amazonian. While no direct evidence for the existence of liquid water was found within the radar detection depth range, we cannot rule out the presence of saline ice in the subsurface of the landing area.

As the largest basin on Mars, Utopia Planitia has both experienced and recorded variations of the Martian palaeoclimate. Layered subsurface structures have been identified by ground-penetrating radar in southern Utopia Planitia but lateral variations of the subsurface, potentially linked to the Martian palaeoclimatic evolution, have not been investigated. Here we report the lateral frequency-variation patterns of Zhurong radar reflections and interpret them as buried polygonal terrain below a depth of 35 m. Sixteen polygonal wedges were identified within ∼1.2 km distance, suggesting a wide distribution of such terrain under Utopia Planitia. The contrast above and below ∼35 m depth represents a notable transformation of aqueous activity or thermal conditions in the Late Hesperian–Early Amazonian. The interpreted buried polygons, possibly generated by freeze–thaw cycles, imply that there was a strong palaeoclimatic variability at low-to-mid latitudes (∼25° N), potentially due to the high obliquity of ancient Mars.

The traditional methods for ground-penetrating radar (GPR) wavefield simulation suffer from numerical dispersion phenomena and difficulties for models with complex geometry. We propose to use the physics-informed neural networks (PINNs) to solve the GPR wave equation and model the wavefield propagation. Deep fully connected networks are used to approximate the solution to the equations. Automatic differentiation with back propagation is taken to calculate the partial derivatives. The loss function combining the wave equation, boundary condition, and temporal constraints is constructed and minimized to train the PINNs. As the solver is mesh-free, this simulation method avoids numerical dispersion artifacts and is very flexible in implementation for inversion. The numerical examples indicate that the PINNs solver for GPR wavefield simulation has high accuracy and efficiency. We also explore the possibility of inverting the electrical parameters from the observed wavefields. The results demonstrate the promising potential of the proposed approach in GPR forward modeling.

In the implementation of migration of multiples, reverse time migration (RTM) is superior to other migration algorithms because it can handle steeply dipping structures and offer high-resolution images of the complex subsurface. However, the RTM results using two-way wave equation contain high-amplitude, low-frequency noise and false images generated by improper wave paths in migration velocity model with sharp velocity interfaces or strong velocity gradients. To improve the imaging quality in RTM of multiples, we separate the upgoing and downgoing waves in the propagation of source and receiver wavefields. A complex function involved with the Hilbert transform is used in wavefield decomposition. Our approach is cost effective and avoids the large storage of wavefield snapshots required by the conventional wavefield separation technique. We applied migration of multiples with wavefield decomposition on a simple two-layer model and the Sigsbee 2B synthetic data set. Our results demonstrate that the proposed approach can improve the image generated by migration of multiples significantly.

The separation of upgoing and downgoing wavefields is an important technique in the processing of vertical seismic profiling data and ocean bottom cable data. It is also used in reverse time migration (RTM) based on the two-way wave equation to suppress low-frequency, high-amplitude noises and false images. Therefore, we model upgoing and downgoing wavefields directly in the wavefield propagation. There are several methods to obtain separated wavefields. The methods using the Fourier transform require storage of the wavefields, which is not practical due to the extremely high disk-space requirements. Methods using Poynting vectors have an ambiguity problem when crossing a peak or a trough of the wavefields. To improve the accuracy and stability of the modeled upgoing and downgoing wavefields in a complicated velocity model, we evaluate an efficient forward-modeling approach purely based on the Hilbert transform in 3D acoustic wavefield simulation. This method is implemented by the Hilbert transform along the time and depth axis, instead of the Fourier transform. We explicitly derive the formulas for upgoing and downgoing wavefield propagation and attach reproducible source codes. Applications to synthetic models indicate that this method can forward propagate upgoing and downgoing wavefields effectively and improve the imaging quality in migration. This method has various potential applications, e.g., 3D seismic imaging with high computation efficiency.

To accurately image the geologic structures from walkaway vertical seismic profile (VSP) data, it is necessary to estimate the subsurface velocity field with high resolution and enhanced illumination of deep reservoirs. Because full-waveform inversion (FWI) suffers from cycle skipping when the initial model is far from the true model, we have adopted frequency-dependent reflection wave-equation traveltime inversion (FRWT) to generate the background velocity model for VSP migration. The upgoing reflection data are separated from the original shot gathers, and dynamic warping is used to evaluate the traveltime differences between the observed data and the calculated data. Different frequency bands of the data are inverted in sequence to reconstruct the velocity model with higher resolution. We also implement wavefield decomposition on the gradient field to extract the contributions of reflection components and improve the updated model. The inverted results obtained from FRWT can be used as the initial model for conventional FWI or the velocity model for reverse time migration. The experiments on synthetic data and field data demonstrate that our approach can effectively recover the background velocity model from walkaway VSP data.