In this study, the authors first incorporate the multigrid implementation strategy into the non-linear least squares-based three-dimensional variational data assimilation system (NLS-3DVar), which is applied for temperature data fusion. A merged temperature dataset at 1° resolution and 6-hour interval is produced based on in situ observations at 2400 observational sites over China and NCEP (National Centers for Environmental Prediction) final global tropospheric analyses. Another set of independent validation data (from January to December except April and May in 2014) is used to evaluate the merged dataset. The dataset of NLS-3DVar is compared with the gridded data at 1° resolution produced by the widely used Cressman interpolation method. NLS-3DVar product always has lower RMSE (Root Mean Square Errors) of 1.961℃ d^-1 and higher correlation coefficient of 0.924 compared to the dataset produced by Cressman interpolation. The precision of merged temperature product of NLS-3DVar is higher in most stations and independent of validation data, especially at those stations in Xinjiang, Gansu, Yunnan, Shanxi, and so on. The performances of NLS-3DVar based on both the single grid and multigrid strategies are also compared. Both RMSE and correlation coefficient have little differences. Although multigrid-based NLS-3DVar uses the sparse process, the precision is almost the same as that of single-grid based NLS-3DVar. However, its computational costs are greatly reduced due to the sparse process. Compared with the STMAS algorithm (Space-Time Multiscale Analysis System), multigrid-based NLS-3DVar performs better regarding the precision of product with almost the same computational efficiency.