Determining crustal thickness is crucial for understanding regional crustal deformation and material cycling. However, in subduction accretionary zones, complex deformation and rapid changes in crustal thickness over narrow spaces often lead to poor imaging quality of Moho surfaces, thereby hindering the accurate determination of crustal thickness through seismic methods. Gravity data are frequently used to infer crustal thickness in these zones due to the dearth of alternative geophysical data and the fact that the conventional Parker- Oldenburg algorithm often necessitates low- pass filtering to ensure stability during the iterative calculation process. However, inappropriate filtering can lead to overly smooth inversion results, creating discrepancies between gravity inversion and seismic data regarding crustal thickness. In this paper, we utilize satellite gravity data and a substantial number of ocean bottom seismometer (OBS) profiles to develop a model for quantitatively estimating crustal thickness, taking into account a range of influencing factors. We apply this model to invert the crustal thicknesses of the Manila subduction zone and the Huatung basin. In comparison to the conventional Parker- Oldenburg method, the BP neural network gravity inversion method demonstrates notable advantages in regions characterized by significant crustal thickness variations. This approach aligns more closely with OBS profiles and is more consistent with the actual geological situation. Our results demonstrate that the crustal thickness of the South China Sea basin ranges from 4 to 9 km. A localized thickening to over 10 km occurs near the residual spreading center of the eastern basin, extending SW- NE towards the Manila trench. The crustal thickness of the eastern basin is 4 to 7 km, which is assumed to be oceanic crust. Furthermore, analysis of crustal thickness symmetry in different tectonic units suggests that the Gagua Ridge may be a product of subduction from the West Philippine basin to the Huatung basin.