The distribution of fractures is highly uncertain in naturally fractured reservoirs (NFRs) and may be predicted by using the assisted-history-matching (AHM) that calibrates the reservoir model according to some high-quality static data combined with dynamic production data. A general AHM approach for NFRs is to construct a discrete fracture network (DFN) model and estimate model parameters given the observations. However, the large number of fractures prediction required in the AHM process could pose a high-dimensional optimization problem. This difficulty is particularly challenging when the fractures form a complex multi-scale fracture network. We present in this paper an integrated AHM approach of NFRs to tackle these challenges. Two essential ingredients of the method are (1) a 2D fractal-DFN model constructed as the geological simulation model to describe the complex fracture network, and (2) a mixture of multi-scale parameters, built according to the fractal-DNF model, as an inversion parameter model to alleviate the high-dimensional optimization burden caused by complex fracture networks. A reservoir with a multi-scale fracture network is set up to test the performance of the proposed method. Numerical results demonstrate that by use of the proposed method, the fractures well recognized by assimilating production data.