SummaryAlthough researchers have applied many methods to history matching, such as Monte Carlo methods, ensemble-based methods, and optimization algorithms, history matching fractured reservoirs is still challenging. The key challenges are effectively representing the fracture network and coping with large amounts of reservoir-model parameters. With increasing numbers of fractures, the dimension becomes larger, resulting in heavy computational work in the inversion of fractures. This paper proposes a new characterization method for the multiscale fracture network, and a powerful dimensionality-reduction method by means of an autoencoder for model parameters. The characterization method of the fracture network is dependent on the length, orientation, and position of fractures, including large-scale and small-scale fractures. To significantly reduce the dimension of parameters, the deep sparse autoencoder (DSAE) transforms the input to the low-dimensional latent variables through encoding and decoding. Integrated with the greedy layer-wise algorithm, we set up a DSAE and then take the latent variables as optimization variables. The performance of the DSAE with fewer activating nodes is excellent because it reduces the redundant information of the input and avoids overfitting. Then, we adopt the ensemble smoother (ES) with multiple data assimilation (ES-MDA) to solve this minimization problem. We test our proposed method in three synthetic reservoir history-matching problems, compared with the no-dimensionality-reduction method and the principal-component analysis (PCA). The numerical results show that the characterization method integrated with the DSAE could simplify the fracture network, preserve the distribution of fractures during the update, and improve the quality of history matching naturally fractured reservoirs.