## 1 Introduction  
The Monte Carlo method has widely been used for geophysical source parameter estimation in remote sensing [1, 2, 3]. The method repeats nonlinear optimization from multiple random starting points to derive the best-fit source parameters that minimize an objective function indicating the misfit between geodetic measurement and the dislocation model. Therefore, it enables avoiding local minima and evaluating the uncertainties of parameters [4]. Interferometric synthetic radar (InSAR) has been one of the most common methods for surface displacement measurement. It compares the phase offset of two or more complex-valued SAR images obtained from different times and locations [5]. InSAR can derive precise geophysical information with centimetric or millimetric accuracy. For this reason, InSAR data have been widely employed in remote sensing, including earthquakes, volcanic activities, landslides [6], thaw-derived slope failure [7], or glacial ice movement [8] investigations. The dislocation model mathematically calculates the surface deformation of earthquakes or volcanic activities. The rectangular dislocation model, also known as the Okadamodel, is a popular model for earthquake source modeling, which assumes finite rectangular source and isotropic half-space [9; 10]. Other models include the prolates spheroid model [11] and the spherical source model [12] for magma source modeling. Most nonlinear optimization algorithms iteratively search local minima using gradient or Hessian approximation of the objective function. However, since the dislocation model computes three-dimensional (3D) surface displacement of points in a two-dimensional (2D) mesh, computation of the objective function leads to many floating-point operations proportional to the mesh size. Therefore, gradient or Hessian approximation of the objective function becomes too expensive. Furthermore, the nonlinear optimization is iterated more than 1000 times for the Monte Carlo method, and thus the Monte Carlo method takes an extreme computational burden. In the recent decade, graphics processing units (GPUs) have been gaining attention as a powerful solution to overcome the performance limitations of conventional multi-core CPUs, especially for applications demanding a massive amount of computations. Furthermore, the release of NVIDIA's computed unified device architecture (CUDA) [13], a new general-purpose parallel programming interface, also facilitated the use of GPU platforms for general-purpose computing applications. Since the CUDA architecture can reorganize the computation around data when the data can be processed independently, we can execute the independent computation in parallel on a vast number of threads to improve the overall performance significantly. Therefore, the Okada model, the dislocation model discussed in this paper, fits well for the GPU implementation because the 3D displacement computation of each point needs a massive number of floating-point operations and has no dependencies between these operations. This paper proposes a performance acceleration framework based on CUDA GPU for the earthquake source estimation with the Okada dislocation model. For this purpose, (1) we analyze the effects of various combinations from options of CUDA optimization techniques and, based on the analysis, perform the CUDA kernel optimization for the parallel implementation of the Okada model; (2) we investigate the performance and correctness of line search algorithms used for the nonlinear optimization to derive the earthquake source parameters in an efficient manner; (3) finally, we verify that the earthquake source parameters derived by our proposed approach also fit the geodetic data correctly by visualizing the residual data between the geodetic data and our modeled data. Considering the target application characteristics, three techniques, common subexpression elimination, thread block size adjustment, and constant caching, are employed for the CUDA kernel optimization.