As the main approach to studying the rock mechanical properties and revealing the evolution of meso-cracks, numerical simulation methods are widely adopted and discussed. The existing numerical simulation parameter calibration methods are mainly trial-and-error methods and orthogonal test methods, both of which fail to fully consider the influence of mesoscopic parameter interactions, resulting in inaccurate simulation results and inconsistency with laboratory tests in terms of macroscopic failure modes. Therefore, the numerical analysis method combining orthogonal-response surface method is adopted. Firstly, the parallel bond model (PBM) mesoscopic parameters with significant influence are selected by orthogonal test. Secondly, the response surface method (RSM) is used to study the influence of the interactions between mesoscopic parameters on the macroscopic parameters of the model sample. Finally, a set of PBM parameter calibration process is proposed in combination with the macroscopic failure morphology of the rock. The results show that effective modulus (E*), stiffness ratio (kn/ks) and their interactions have significant effects on elastic modulus (E); kn/ks, contact friction coefficient (μ) and minimum particle radius (Rmin) have a significant effect on Poisson’s ratio (ν); the influence of cohesion (c), normal bond strength (σc) and their interactions on uniaxial compressive strength (UCS) is significant. The error between the simulated results by RSM and experimental results of the references of the mesoscopic parameters is less than 7%, the characteristics of the stress-strain curves of the two are consistent, and the macroscopic failure modes are the same, which proves that the proposed calibration process of PBM mesoscopic parameters is accurate and reliable.