The quadratic Wasserstein metric has shown its power in comparing probability densities. It is successfully applied in waveform inversion by generating objective functions robust to cycle skipping and insensitive to data noise. As an alternative approach that converts seismic signals to probability densities, the squaring scaling method has good convexity and thus is worth exploring. In this work, we apply the quadratic Wasserstein metric with squaring scaling to regional seismic tomography. However, there may be interference between different seismic phases in a broad time window. The squaring scaling distorts the signal by magnifying the unbalance of the mass of different seismic phases and also breaks the linear superposition property. As a result, illegal mass transportation between different seismic phases will occur when comparing signals using the quadratic Wasserstein metric. Furthermore, it gives inaccurate Fréchet derivative, which in turn affects the inversion results. By combining the prior seismic knowledge of clear seismic phase separation and carefully designing the normalization method, we overcome the above problems. Therefore, we develop a robust and efficient inversion method based on optimal transport theory to reveal subsurface velocity structures. Several numerical experiments are conducted to verify our method.