We have developed an original method for migration velocity analysis for vertical seismic profile (VSP) data. The method is based on VSP imaging with the primary reflected data and surface-related multiples. Velocity updates are carried out by matching the images obtained by these two types of waves. To do it automatically, we have introduced an objective function of the crosscorrelation type and developed a gradient-based maximization algorithm. On synthetic walk-away VSP data sets, we found that maximization of the crosscorrelation of the images provided by the two different types of waves allowed the retrieval of the interval velocities below the borehole receivers — a problem that was considered to be challenging with a VSP acquisition. Also, the possibility of the lateral velocity gradient reconstruction was developed on a simple synthetic model. On the real walk-away VSP data set from the North Sea, we found that the method could also be used for estimation of effective anisotropy parameters below the borehole receivers. In all of our examples, the velocity model refinement based on the presented method led to significant improvements of the seismic images.