Developing a forecast model of soil liquefaction helps assess the condition of the site and prevent the pertinent losses resulting from an earthquake. In recent years, various types of datasets have been utilized for this purpose, such as field test results (Standard Penetration Test ( \(SPT\) ), Cone Penetration Test ( \(CPT\) ), Shear Wave Velocity ( \(Vs\) ), etc.), laboratory tests (Fine Content ( \(FC\) ), Unconfined Compressive Strength ( \(qu\) ), etc.), and analytical equations (Cyclic Stress Ratio ( \(CSR\) ), Cyclic Resistance Ratio ( \(CRR\) ), etc.). The present study utilized field test data and group method of data handling ( \(GMDH\) )-based metaheuristic models, taking into account the complexity of the liquefaction mechanism. For this aim, Aquila Optimizer ( \(AO\) ), Arithmetic Optimization Algorithm ( \(AOA\) ), Genetic Algorithm ( \(GA\) ), Gravitational Search Algorithm ( \(GSA\) ), and Particle Swarm Optimization ( \(PSO\) ) have been employed to tune the hyperparameters of the \(GMDH\) and artificial neural network ( \(ANN\) ) models. The results showed that \(GSA-ANN\) and \(GA-GMDH\) have the best performance in predicting liquefaction triggering potential, while \(AOA-GMDH\) has the lowest accuracy. Also, a Shapley Additive Explanations ( \(SHAP\) ) analysis was done using \(GSA-GMDH\) , which showed that the parameters depth ( \(Z\) ), \(\sigma \) , and \(\sigma {\prime}\) have the most significant impact on soil liquefaction. All models greatly surpass the performance of a random classifier ( \(AUC\) = 0.5), therefore validating their predictive efficacy. The elevated \(AUC\) values of \(AOA-GMDH\) and \(GA-GMDH\) indicate that these hybrid methodologies are especially efficacious for the job, presumably owing to their refined parameter optimization and resilient feature selection processes.