A four-dimensional ensemble-variational hybrid data assimilation (DA) system based on the dimension-reduced projection (DRP) technique was developed for global numerical weather predictions (NWPs). Instead of the adjoint technique, an ensemble approach is utilized in this technique to calculate the gradient of the cost function of the standard four-dimensional variational (4DVar) DA. The flow-dependence of the background error covariance (BEC) is realized in the variational configuration by dynamically updating the initial perturbation samples during the assimilation cycle. A limited number of leading eigenvectors of the correlation function of localization are selected to filter out the spurious correlations in the BEC matrix (B-matrix). A linear combination of the ensemble analysis increment sample with the random perturbation sample satisfying balance constraints is used as the inflation technique to prevent the BEC from underestimation and to achieve the hybrid of the flow-dependent and static B-matrices when updating the initial perturbations. In order to evaluate the new system, single-point observation experiments (SOEs) and observing system simulation experiments (OSSEs) were conducted with sounding and cloud-derived wind data. The flow-dependent characteristic was verified by the SOEs that utilized localized ensemble covariance. In the OSSEs, DRP-4DVar produced better analysis than 4DVar. Moreover, the ensemble mean forecast with DRP-4DVar analyses as the initial conditions reduced errors in geopotential height, 24-h accumulated precipitation and other variables relative to the 4DVar-based forecast. Significant improvements of analysis and forecast achieved in the data-sparse Southern Hemisphere by DRP-4DVar indicate a remarkable advantage of the ensemble covariance in sufficient use of sparse observations.