Finite‐fault inversion is a widely used method to study earthquake rupture processes. Previous studies proposed many approaches to implement finite‐fault inversion, including time‐domain, frequency‐domain, and wavelet‐domain methods. Among these methods, wavelet‐domain methods can analyze both the time and the frequency information of waveforms at different scales using the wavelet transform. Traditional wavelet‐domain methods based on the simulated annealing (SA) approach implement finite‐fault inversion with lower‐ and higher‐frequency signals simultaneously to recover both larger‐scale and smaller‐scale characteristics of the rupture process. In this article, we propose a new two‐stage SA strategy to recover larger‐scale characteristics first and then smaller‐scale features based on the multiscale wavelet‐domain method. During the first stage, we use only the wavelet scales corresponding to lower‐frequency waveforms to construct the objective function to obtain a relatively reliable larger‐scale slip model. Subsequently, we initiate the second stage by gradually adding the wavelet scales corresponding to higher‐frequency parts into the objective function to resolve smaller‐scale slip features. Our synthetic tests show that this new multiscale SA strategy works better at converging to the global minimum solution than the original SA approach. We also apply this new strategy to the 2015 Gorkha, Nepal, earthquake. The results show that our recovered slip distribution of this event is associated well with high‐frequency radiation subevents from previously published backprojection studies. These large high‐frequency subevents appear to be located at the beginning, transition, or termination of the main slip patches obtained in this study.