A new inversion scheme is presented to obtain three-dimensional images of P-wave velocity(Vp) and P–S-wave velocity ratio(Vp/Vs) using P- and S-phase pairs, i.e., the same source–receiver pairs for the P- and S-wave arrival-time data. The S-wave velocity(Vs) was separately inverted using the S-phase arrival times. The earthquake hypocenters were simultaneously relocated in the joint inversion. The method considers the Vp/Vs anomaly as a model parameter in the inversion. The proposed method thus provides a more robust calculation of the Vp/Vs anomaly than the conventional method of dividing Vp by Vs. The method also takes into account the ray path difference between P- and S-waves, and hence yields a less biased Vp–Vs ratio than the method of inverting S–P-wave data for Vp and Vp/Vs anomalies under the assumption of identical P and S ray paths. The proposed method was used to image the crust and upper mantle in northeastern(NE) Japan taking advantage of a large number of high-quality arrival times of P- and S-wave source–receiver pairs. The inverted structures suggest that the subducting slab of the Pacific plate is an inclined zone of high-Vp and Vs anomalies with low Vp/Vs perturbation. The mantle wedge is characterized by low-Vp, low-Vs, and high-Vp/Vs anomalies at shallow depths beneath active volcanoes. These features are also observed at greater depths in the back-arc region. Although these features have been previously reported, the Vp/Vs anomaly pattern obtained in this study shows much less scatter and is much better correlated with the seismic velocity perturbation patterns than previous studies. The proposed method can be used, in conjunction with velocity anomaly patterns, to quantify thermal processes associated with plate subduction.