We have developed a method for three-dimensional simultaneous velocity and density inversion using traveltimes of local earthquakes and gravity data. The purpose of this method is to constrain the velocity inversion and increase the spatial resolution of shallow velocity structures by introducing additional gravity data. The gravity data contributes to the P-and S-wave velocity models by imposing constraints between seismic velocities and density. The constraint curve is constructed so as to fit the data for porous rock samples, and deviations from the curve are taken into account in the inversion. The constraint is imposed at only the first layer, because density structure is well resolved at shallower parts and it is difficult to determine uniquely at greater depths. Synthetic inversion tests indicate that gravity data can improve the resolution of the velocity models for this layer. The method is applied to investigate the subsurface structure of Izu-Oshima volcano, Japan and velocity structures with high spatial resolution are obtained. The additional gravity data contribute primarily to improvement of the S-wave velocity model. At 0.25 km depth, a high velocity anomaly due to caldera-filling lava flows is observed. At 1.25 and 2.5 km depths, high velocity intrusive bodies are detected. A NW-SE trending high velocity belt at 1.25 km depth is interpreted as being caused by repeated intrusion of dikes.