The numerical approximation of Gaussian integrals is of great relevance for a wide variety of statistical calculations. In recent decades, numerous efficient cubature rules of high polynomial exactness for the case of high-dimensional Gaussian integrals have been proposed. In addition to the question of efficiency, however, the question of stability, which depends on the influence of negative weights, is also of crucial importance. For certain degrees of polynomial exactness no stable rules are known. Therefore, a new approach based on metaheuristic optimization is presented which allows the computation of cubature rules which significantly exceed comparable rules known from literature with respect to stability and efficiency. Especially with regard to iterative numerical calculations, the use of unstable cubature rules over time can lead to a distortion and thus a significant deterioration of the results. To illustrate this fact,
the cubature-based Kalman filter is utilized. Based on the results of simulation studies, the necessity for the use of stable cubature rules is motivated.