A computer-implemented method for minimizing artifacts in a seismic image of a subsurface region of interest, the image based on a data beam set derived from recorded seismic data and a modeled beam set derived at least in part from a velocity model related to a subsurface region and the seismic image, the method comprising: determining (20, 22) a cross-correlation function between the modeled beam set and data beam set over a predetermined range of time shifts; performing a Monte Carlo inversion including a back-projection and a global minimization to update the velocity model; and forming an updated seismic image of the subsurface region based on the updated velocity model, wherein the Monte Carlo inversion comprises: using the cross-correlation function to determine a probability distribution function over the predetermined range of time shifts; using the probability distribution function to determine a cumulative distribution function over the predetermined range of time shifts; generating a random number within a predetermined interval; selecting a time shift value from the predetermined range of time shifts such that the randomly generated number is equal to the cumulative distribution function at the selected time shift value; and using the selected time shift value to compute a velocity correction for the velocity model by the back-projection, wherein the global minimization is simulated annealing; and wherein temperature and time interval parameters of the Monte Carlo inversion are reduced from initial values to predetermined final values over successive iterations of the simulated annealing, thereby reducing the presence of artifacts, due to cycle skips or spurious alignments with coherent noise, in the updated seismic image formed based on the updated velocity model. determining (20, 22) a cross-correlation function between the modeled beam set and data beam set over a predetermined range of time shifts; performing a Monte Carlo inversion including a back-projection and a global minimization to update the velocity model; and forming an updated seismic image of the subsurface region based on the updated velocity model, wherein the Monte Carlo inversion comprises: using the cross-correlation function to determine a probability distribution function over the predetermined range of time shifts; using the probability distribution function to determine a cumulative distribution function over the predetermined range of time shifts; generating a random number within a predetermined interval; selecting a time shift value from the predetermined range of time shifts such that the randomly generated number is equal to the cumulative distribution function at the selected time shift value; and using the selected time shift value to compute a velocity correction for the velocity model by the back-projection, using the cross-correlation function to determine a probability distribution function over the predetermined range of time shifts; using the probability distribution function to determine a cumulative distribution function over the predetermined range of time shifts; generating a random number within a predetermined interval; selecting a time shift value from the predetermined range of time shifts such that the randomly generated number is equal to the cumulative distribution function at the selected time shift value; and using the selected time shift value to compute a velocity correction for the velocity model by the back-projection, wherein the global minimization is simulated annealing; and wherein temperature and time interval parameters of the Monte Carlo inversion are reduced from initial values to predetermined final values over successive iterations of the simulated annealing, thereby reducing the presence of artifacts, due to cycle skips or spurious alignments with coherent noise, in the updated seismic image formed based on the updated velocity model.