A complete set of 13 morphological parameters in accordance with standard ISO 4287 was applied to quantifying a series of random profiles. These profiles were imported into a discrete numerical model to perform 480 interface shear tests on coarse-grained soils. The relevant morphological parameters were selected using Spearman's rank correlation coefficient for model selection. An optimal metaheuristic model was developed using a genetic algorithm and was further compared with the existing predicted formulas. The 2D discrete element method (DEM) results indicate that the highest correlation with shear strength was obtained for the hybrid parameter Pdq which represents not only the amplitude information but also the surface slope information on a random surface. The optimal model with one significant input variable (Pdq) was effectively selected through the Bayesian nonparametric general regression analysis. For irregular interface shearing widely existing in most geotechnical engineering, Pdq is more efficient and accurate to quantify the surface morphology or estimate the interface shear strength compared with relative roughness.