Naturally fractured vuggy carbonate reservoirs play an important role in the world's oil and gas production. Pore-scale modeling is treated as an important approach for unveiling the flow mechanism in carbonate reservoirs. The realistic regeneration of fractured vuggy reservoirs needs to preserve the geological structures and complex pores morphologies. Nevertheless, regeneration of such reservoirs is a great challenge because of the co-existence of fractures and vugs over a wide range of scales. In this paper, the reference fracture-vug network will be divided into two systems: matrix-fracture and matrix-vug. The two systems are reconstructed in two ways: (I), the matrix-vug and matrix-fracture systems are generated by improved simulated annealing method (ISAM) and multiple point statistics algorithm (MPS), respectively; (II), both of the matrix-vug and matrix-fracture systems are generated by MPS. When applying for ISAM, the multipoint histogram mainly depends on the N-point template. The template is a unit configuration, which is considered as the characteristic of the geometric features. Based on this, we introduce the correlation coefficient matrix to determine an optimal template, and the accuracy of template selection algorithm is also discussed. The optimal template represents essential geological features of the training image, and its size is much smaller than the template size selected by entropy. Therefore, it can improve the computing efficiency while preserving the key properties of the network. Then, a superimposing method is proposed to create the fracture-vug network by combining the matrix-vug with matrix-fracture models. The superimposing method using ISAM (for matrix-vug model generation) and MPS (for matrix-fracture model generation) is called SMIM, while superimposing method using MPS (for both of matrix-vug and matrix-fracture models generation) is called SMM. Finally, the spatial connectivity, local percolation probability and patterns reproduction of the generated models are used to compare with those of the reference model. The results show that the 7-point template selected by the correlation coefficient matrix is enough to get the spatial features of the matrix-vug system. The generated fracture-vug network using the selected template is able to retrieve the patterns of training image (reference model) more accurately than other used templates. The fracture-vug network generated by superimposing method produces better results than the other methods (all named as “direct method”). The network simulated with SMIM algorithm matches the reference model better than SMM algorithm. The probable reason maybe exist in the different generation methods used for matrix-vug systems. The local percolation probability (LPP), treated as a general geometric characterization method, is used to characterize the connectivity of measured units at a given local porosity. The results of the probability analysis show that SMIM is better in depicting the expected patterns of geological heterogeneities and in producing higher degree of connectivity of pore space than the SMM. The pattern reproduction of SMIM model is assessed by scanning the generated model, and the ratio of pattern reproduction number to that of reference model (i.e. the similarity) is around 97.28%. These results may guide the further application of the proposed method to field-scale modeling of naturally fractured vuggy reservoirs.