This paper introduces a parameter estimation method for a general class of statistical models. The method exclusively relies on the possibility to conduct simulations for the construction of interpolation-based metamodels of informative empirical characteristics and some subjectively chosen correlation structure of the underlying spatial random process. In the absence of likelihood functions for such statistical models, which is often the case in stochastic geometric modelling, the idea is to follow a quasi-likelihood (QL) approach to construct an optimal estimating function surrogate based on a set of interpolated summary statistics. Solving these estimating equations one can account for both the random errors due to simulations and the uncertainty about the meta-models. Thus, putting the QL approach to parameter estimation into a stochastic simulation setting the proposed method essentially consists of finding roots to a sequence of approximating quasiscore functions. As a simple demonstrating example, the proposed method is applied to a special parameter estimation problem of a planar Boolean model with discs. Here, the quasi-score function has a half-analytical, numerically tractable representation and allows for the comparison of the model parameter estimates found by the simulation-based method and obtained from solving the exact quasi-score equations.