Temporal changes of seismic velocities may provide important information on the processes that occur inside the Earth. However, using body wave data with passive sources faces the problem of an uneven distribution of rays, which may cause artifacts with stronger amplitudes than the actual velocity changes in the Earth. We propose an algorithm for the selection of similar data sets in different time periods that minimize the artifacts related to variable data distributions. In this study, we used the data of the Japan Meteorological Agency for several years before and after the Mw 9.0 Tohoku-Oki event that occurred on 11 March 2011. We performed careful testing using different synthetic models, showing that the selected data subsets allow detecting weak velocity changes with amplitudes above 0.2%. The analysis of the experimental data revealed important features associated with the stress and deformation distributions after the megathrust event. In the upper crust, we found a large zone along the coast with significant P velocity increase likely caused by compression of crustal rocks. This zone was cut by several elongated anomalies with local velocity decrease coinciding with the limits of the maximum slip area. These anomalies possibly mark the areas of major ruptures and deformations after the Tohoku-Oki earthquake. In the coupling zone at a depth of 40 km, we observe a velocity decrease in the area of the Mw 7.7 aftershock representing strong fracturing in the focal zone. Beneath the volcanic arc, we observe significant (up to 0.5%) decrease of P velocity but less prominent S velocity changes.