基于加权平均导数的频率-空间域正演模拟及GPU实现

传统基于旋转坐标系的频率-空间域正演模拟方法仅适用于方形网格,而实际生产中矩形网格广泛存在,本文提出一种适用性广的正演差分算子,不仅适用于方形网格而且适用于矩形网格.通过综合运用平均导数法、加速项加权平均、模拟退火法压制频散和减少单个波长所需网格点数,从而提高算法精度和减少计算量.在该方法的基础上采用不完全LU分解作为求解Helmholtz方程的预条件,并利用图形处理器加速计算速度,很大程度上提高了频率域正演的效率.

Frequency-space domain acoustic modeling based on an average-derivative and GPU implementation

Full-waveform inversion (FWI) can supply accurate velocity model for seismic imaging and is becoming a very important technology in seismic exploration. The frequency domain modeling is the base of frequency domain FWI and can be used for each iteration. In this case, accurate and fast frequency domain modeling is very necessary. The traditional frequency domain modeling methods can only be used in square grid and the memory and computation requirements are huge. We present a high efficient and adaptable frequency domain modeling technique in this paper.We design a 17-point scheme to perform the frequency domain modeling and use the average-derivative method to overcome the disadvantage that the traditional frequency domain modeling method can not be used in rectangular grid. Our operator can not only be applied in the condition of square grid but also rectangular grid, which is very useful in the real world situation because the rectangular grid is very popular in field work. Therefore, our operator is of greater applicability than the traditional methods. In order to reduce the dispersion and the number of grid points per shortest wavelength, we combine the mass acceleration term into the point calculation and get the optimal coefficients through simulated annealing algorithm. As the program Matlab provides the large scale optimization method, we use fmincon in this program to further optimize the coefficients. At last, the number of grid points per shortest wavelength is reduced to 2.5 comparing to 4 in 9-point method which is the most common frequency domain modeling method with phase velocity errors no larger than 1%. Therefore, the requirements of calculation amount and memory are reduced greatly. There are two main ways in solving Helmholtz equation which is the algorithm of frequency domain acoustic modeling, including direct method and iterative method. The direct method is LU decomposition while the matrix of LU decomposition in our algorithm usually is huge in 2D data and is difficult to solve using current computer. Therefore, iterative method is an option. However, in iterative method the convergent rate usually is very low, which leads to huge computation. We use the incomplete LU decomposition as the precondition and accelerate the speed of convergence. For the same test, the traditional method needs more than 20000 times to converge to a stable value while our method needs only 6 times, which reduces the amount of calculation greatly. In order to accelerate the process of frequency domain modeling, we apply the Graphics Processing Unit (GPU) because GPU has powerful computing capability. In the GPU, fair use of video memory is very crucial and can lead to higher efficiency. Using the GPU, we obtain twice speed of the traditional CPU method in modeling. There are two main methods in frequency domain acoustic modeling, including 9-point method and 25-point method. We compare our method with those of two methods in accuracy and efficiency. We use the same velocity model and the dx=dz=25 m. The dominant frequency of the Ricker wavelet is 20 Hz. In order to remove all disturbances, we do not apply any boundary condition. Because our method only needs 2.4 grids, to reach the same accuracy our operator can use relatively coarse mesh and it leads to less computing work. Meanwhile, the storage is reduced, too. For example, if the 9-point method needs nx×ny grids, the nonzero elements of coefficient matrix is 9×nx×nx×nz×nz, the 25-point method is 9.76×nx×nx×nz×nz and our method is 6.12×nx×nx×nz×nz. Therefore, our method produces a great savings of memory. In order to test our method in the complex velocity model, we applied it in the Marmousi model the sampling intervals of which are 12.5 m and 4 m. The grid of Marmousi is rectangular and our method can be used in this kind of model while 9-point and 25-point can not. We also give the result of time domain method which is second-order in time and twelfth-order in space. Comparing the result from time-domain method with the result from our technique, we find that there is no significant difference between them. Besides, the boundary condition is very important in modeling, we derive the perfectly matched layer equation based on our 17-point scheme and verify it in the Marmousi model.We propose a 17-point scheme which has the advantage of wide application and high accuracy to perform frequency domain acoustic modeling. It overcomes the disadvantage of traditional methods which can only be used in the square grid. After optimization, the number of grid points per shortest wavelength is reduced to 2.5, which reduces huge memory requirement and computing time. Because of the slow convergence, we apply the incomplete LU decomposition as the precondition in solving Helmholtz equation and accelerate the speed of convergence. As the GPU has the powerful computing capability, we use it to accelerate the process of modeling and save much time on computing.