The medium-resolution spectral imager (MERSI) onboard China’s new generation of polar-orbiting meteorological satellite series FengYun-3 (FY-3) is providing global observations of the Earth’s environment. One of the important characteristics of the MERSI sensor is that it has four bands at 250-m spatial resolution. However, very few high-level land surface and atmospheric products have been generated from MERSI data. This article presents an effective optimization method to estimate the continuous temporal distribution of multiple land surface and atmospheric variables at a 250-m spatial resolution from time series FY-3B MERSI top-of-atmosphere (TOA) observations. The MODIS cloud detection method was applied to the MERSI TOA reflectance to determine the clear observations. Then, a physics-based BRDF correction was adopt to reduce the topographic effects using the slope angle and aspect angle of the slope. Finally, a coupled land surface-atmosphere radiative transfer (RT) model and the shuffled complex evolution (SCE) optimization algorithm were used to estimate a suit of variables, including the leaf area index (LAI), the aerosol optical depth (AOD), the cloud optical thickness (COT), the cloud effective particle radius (CER), the land surface reflectance, the shortwave and visible albedo, the incident shortwave radiation (ISR), the incident photosynthetically active radiation (PAR), the fraction of absorbed PAR (FAPAR), the surface broadband emissivity (BBE), and the TOA shortwave albedo. The method was applied to the data from the 18 SURFRAD, the America FluxNet (AmeriFlux), and the Images and Beijing Normal University net (BNUnet) sites. All the estimated variables were also compared with Moderate Resolution Imaging Spectroradiometer (MODIS) and Global Land Surface Satellite (GLASS) products. The estimated LAI, PAR, FAPAR, shortwave albedo, and ISR were also validated using ground measurements from the selected sites. The results show that our method can simultaneously estimate multiple temporally continuous variables in both clear-and cloudy-sky conditions, with an accuracy comparable to those of MODIS and GLASS products. The estimated LAI, PAR, FAPAR, shortwave albedo, and ISR are consistent with the field measurements, with coefficients of determination ( <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$R^{2}$ </tex-math></inline-formula> ) of 0.692, 0.783, 0.788, 0.559, and 0.833, and root mean square errors (RMSEs) of 0.427, 56.681 W/m <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> , 0.092, 0.084, and 115.305 W/m <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> , respectively. The proposed algorithm can effectively estimate the instantaneous atmospheric and land surface variables directly from satellite observations without cumbersome atmospheric corrections. Moreover, it can potentially be applied to multiple satellite observations.