A two-level radiosity algorithm was developed here for the computation of global radiation transfer at complex land cover called HRAD, which rapidly simulates the surface leaving radiances in three dimensional (3D) landscapes. The algorithm is an integration of a ray casting module and an adapted hierarchical radiosity method with heterogeneous volume clusters. The time cost of our algorithm was nearly linear to the polygon number. The algorithm was first validated by comparing the simulated radiances with that of Radiosity-Graphics combined Model (RGM) at small scenes. The single scattering results were highly correlated with that of RGM (R2 > 0.99 and RMSE <; 0.005). The multiple scattering results were also perfect for visible bands (R2 > 0.99 and RMSE <; 0.01), and acceptable for NIR (R2 > 0.98 and RMSE = 0.015). By using the Gauss scene in RAMI-III website, we also validated our algorithm.