## 1 Introduction  
As an important branch of geophysics, seismic exploration uses seismic waves excited on the Earth's surface from seismic sources, such as explosions, and receives seismic waves reflected on the Earth's surface from underground to image underground structures. It is widely applied in the exploration of oil, gas, minerals, and the probing of urban underground structures [1]. As an important step in seismic data processing, migration image processing obtains image profiles of underground structures by using acquired seismic data, which are then used to explain geological structures. The imaging methods of exploratory seismic data include the Kirchhoff method, the one-way wave equation method, and reverse time migration (RTM) [2, 3], of which RTM currently enjoys the highest precision. In RTM calculations, although the velocity model of underground structures is one of the vital factors influencing imaging precision, a pre-estimated initial velocity model is often used to migrate, then the initial velocity model is updated with the imaging results, thus generating an iterative process in the calculation of migration and velocity models due to the precise velocity of underground models being unknown to us. For such iterations, an intermediary data volume is needed, which we refer to as surface offset gathers (SOGs) [4, 5]. Currently, RTM-based SOGs are a research hotspot. Among these, the shot-domain RTM, the simplest method to output SOGs, divides receiver data into a few parts by offsetting and then migrating them, respectively, and uses the mean offset to identify different results so as to generate SOGs. In this method, it is necessary to perform a migration calculation on data from every shot, so the calculation intensity is too large for it to be implemented on industry scale [6]. The improved algorithm by [PERSON] et al. only calculates a minority of sections of different offsets, and then supplements the necessary data in offset gathers through interpolation [7]. Another method to output SOGs is the wave equation method, called wave-equation Kirchhoff migration, which is named as such because the Born equation it uses is similar to the Kirchhoff integration equation. The difference between them is that in the calculation, the wave equation, rather than ray tracing, is used in the Green's function for propagation between shots and receivers [8; 9]. The attribute migration maps the geometrical and physical attributes in seismic data, such as travel time, angles of reflection and incidence, and shot-receiver distance, to the migration results through migration calculations performed twice, and extracts attributes from the ratios of the results from imaging performed twice. This was proposed by [PERSON] [10] and applied in the calculation of three-dimensional (3D) stereotomography [11]. [PERSON] et al. and [PERSON] et al. discussed the possibility of applying the method to output the attributes of the shot-receiver distances [12; 13]. In this study, based on attributive migration, we investigated how to apply this method to the calculation of RTM SOGs. The seismic migration calculation enjoys high computational intensity, and both industry and academia have carried out much research on the acceleration of seismic imaging. Owing to their multiple kernels and bandwidth, as well as advantages of low power consumption and small occupied space, graphic processing units (GPUs) have been introduced in general computations [14]. The application of NVIDIA GPUs in seismic image processing has significantly improved computing efficiency, for instance, in time migration and RTM and waveform inversion. [15; 16; 17], especially RTM, because it has the highest calculation accuracy in the seismic migration but is also the most time-consuming.