Shear wave splitting (SWS) analysis has been widely employed for fracture characterization in both global seismology and seismic exploration. Two key SWS attributes—fast shear wave polarization and the time delay between fast and slow shear waves—can be inverted from four-component seismic data (two horizontal sources and two horizontal receivers). These SWS attributes enable the characterization of subsurface fracture parameters, such as fracture strike and density. In this study, a nine-component vertical seismic profile (VSP) survey was acquired in the Sanhu Depression of the eastern Qaidam Basin, northwestern China. Preliminary analysis of the shear-wave source VSP data reveals two distinct SWS signatures at different depths, corresponding to two separate fractured layers. However, characterizing multiple fractured layers presents significant challenges, as the SWS attributes of deeper fractured layers are strongly influenced by those of overlying fractured formations. Existing approaches for predicting multi-layer fracture parameters are predominantly data-driven and are largely limited to qualitative analysis. To address these challenges, we propose a robust, rock-physics-model-guided method that enables the quantitative estimation of both the fracture strike and fracture density in multiple fractured layers. First, the parameters of the shallow fractured layer are directly estimated from the SWS attributes. Then, synthetic VSP records of the deeper fractured layer are modeled by incorporating Hudson's theory and the reflectivity method. The fracture parameters of the deep fractured layer are inverted by minimizing the difference between the SWS attributes of synthetic records and those of the actual seismic data. A hierarchical search strategy (coarse-scale + refined-scale) is employed to accelerate convergence toward the optimal solution. This investigation provides a practical tool for quantitative characterization of subsurface formations with multiple fractured layers.