In geophysical applications, inverse problems arise in order to characterize the earth structure using seismic data. We focus in least squares methods for the joint inversion of different datasets from the same region. An improvement in the numerical results is expected when using the constraints given by each independent dataset versus one dataset. We are conducting a study of some regularization techniques and Primal-Dual Interior-Point methods for incorporating restrictions in the parameterized subsurface layered velocities. We present a numerical experimentation using synthetic data from teleseismic P-wave receiver functions and surface wave dispersion velocities, and conclude which strategy characterize better a region with lower CPU time and higher accuracy.