Physically based multispectral image simulation consists of sensor system modeling, bottom-of-atmosphere (BOA) image generation, and top-of-atmosphere (TOA) image calculation. TOA radiance images are usually generated using a lookup table (LUT) for computational efficiency, which is calculated by means of atmospheric radiative transfer codes with different combination of input variables, including viewing zenith, solar zenith, and relative azimuth angles; visibility; columnar water vapor; and ground elevation. In this paper, a new strategy is proposed for TOA radiance image simulation, where transmitted surface radiance and atmospheric radiance at the TOA are calculated, respectively, to improve accuracy as well as efficiency. The transmitted surface radiance image is obtained from pixel-by-pixel calculation of BOA radiance and path transmittance. In calculating the atmospheric radiance of TOA, two LUTs are built for the emitted and the scattered radiance from each atmospheric layer, respectively. The effects of visibility and columnar water vapor on the atmospheric radiance are characterized by means of an equivalent path transmittance, which is related to the scene geometry as well as the thickness of atmospheric layer. In this way, when a new scene is simulated, except for three variables, i.e., viewing and solar zenith angles and atmospheric layer number, other parameters are set as constants in building the LUTs, enabling more combinations of input variables without adding excessive computational burden. Multispectral images in different bands with moderate spatial resolution are simulated and compared with the moderate-resolution imaging spectroradiometer (MODIS) images to demonstrate the accuracy and the usefulness of the proposed strategy.