Spatially distribution information on wheat yield forecasting at the large regional scale is important for market forecast and agricultural sustainable development. Assimilating remote sensing information into the crop growth model has demonstrative to be the effective approach for crop yield prediction. However, it remains a challenge to determine the crop growth model of input parameters and initial conditions at the spatial regional scale. In the paper, we proposed a Markov Chain Monte Carlo (MCMC) and 4DVAR hierarchical data assimilation scheme, which enables the winter wheat yield forecasting at the 500 m grid size ahead of harvest time in Henan province. This approach applies data assimilation algorithms at two spatial scales. At the county scale, the MCMC algorithm was used to recalibrate the uncertain and sensitive parameters of the WOFOST model using the county-level statistical yield. Then, we assimilated time-series MODIS reflectance into WOFOST-PROSAIL coupled model using the 4DVAR method for each 500 m pixel across the entire Henan province of China. The results show that the simulated yield was strong positive correlated with the statistical yield at county-level scale with R2 = 0.81 and RMSE = 877 kg/hm ² , which demonstrated the potential usage of the MCMC-4DVAR based large area yield estimation with remote sensed data and yield statistics.