Gravity inversion, assuming a two-layered structure consisting of basement and upper sediment layers, is commonly used to estimate the subsurface structure in geothermal resource exploration. In geothermal development, clarifying the fracture distribution, which can be a fluid pathway, is important for determining drilling targets and elucidating accurate geological structures that contribute to a precise reservoir numerical simulation. However, the assumption of a two-layer structure may lead to misinterpretation because geothermal areas generally have complex stratigraphy, including different rock bodies such as volcanic and sedimentary rocks. Therefore, we developed a gravity inversion procedure considering a multi-layer configuration. To reduce the non-uniqueness of the model and obtain a solution close to the actual lithology, this method introduces smoothing terms, absolute constraints based on borehole and surface geology data, and a weighting factor for adjusting the sensitivity of gravity to the target depth. Previous studies have presented that the initial model is crucial to obtaining a satisfactory solution using multi-layer inversion. They also proposed obtaining an appropriate initial solution by delineating a geological cross-section with reference to active seismic surveys and borehole data. However, seismic exploration of geothermal resources has not yet become widespread and is still under development due to the difficulty of implementation and interpretation. Therefore, we adopted a horizontal multi-layer structure as the initial model, and the boundary depths between them were obtained from a power spectrum analysis of the gravity anomalies. Synthetic tests assuming half-graben or graben structures suggested that for obtaining an appropriate solution using this methodology, the initial model must assign the proper thickness to each layer, whereas the setting of the absolute layer depths was not important. The proposed procedure was applied to actual gravity anomalies observed in the Onikobe geothermal field, in Japan. The estimated subsurface density structure was consistent with previous literature, demonstrating the effectiveness of this method.