The　global　spatial　distributions　of　direct　solar　radiation　(DSR)　on　celestial　surfaces　over　rugged　topography　serve　as　a　fingerprint　of　their　surface　energy　and　climatic　environment.　The　widespread　availability　of　digital　elevation　models　(DEMs)　presents　an　inspiring　opportunity　to　simulate　and　further　understand　the　planetary　solar　radiation　environment.　In　response,　we　develop　and　assess　a　DEM-based　theoretical　framework　for　accurately　simulating　DSR　on　celestial　surfaces.　Building　upon　this　framework,　we　draw　the　global　spatial　distribution　map　of　DSR　for　the　Moon,　Mars,　Ceres,　and　Mercury.　The　proposed　method　demonstrates　a　satisfactory　capacity　to　restore　the　visual　effect　of　planetary　aerial　images.　Validation　against　meteorological　data　from　Earth＇s　stations　reveals　a　Root　Mean　Squared　Error　of　22.43　MJ/m2　and　a　bias　of　-7.62　MJ/m2.　Subsequently,　we　performed　an　exhaustive　spatial　analysis,　encompassing　both　qualitative　and　quantitative　assessments　of　the　DSR　distribution,　terrain　influences,　and　the　underlying　driving　mechanisms.　Specifically,　DSR　distribution　exhibits　latitudinal　anisotropy　in　intensity　and　terrain　effects　across　all　celestial　bodies,　while　Mercury　displays　longitudinal　anisotropy,　governed　by　its　axial　tilt　and　rotational　dynamics.　We　found　that　the　area　proportion　of　the　permanently　shadowed　regions　(PSR)　is　potentially　correlated　with　the　axial　tilt　angle,　while　the　PSR　was　underrated　for　Ceres　in　previous　studies.　Using　geographical　detectors　and　interpretable　machine　learning,　we　quantitatively　assessed　the　contribution　of　various　factors　to　the　spatial　distribution　and　intensity　of　DSR.　Astronomical　factors　(the　sunshine　duration　or　mean　solar　zenith　angle　at　noon),　exert　the　greatest　influence　on　the　spatial　distribution　of　DSR,　followed　by　geographical　position　and　terrain　factors.　However,　terrain　factors　(terrain　openness)　contribute　the　most　to　DSR　intensity,　followed　by　astronomical　factors　and　geographical　position　factors.