In this study, a novel approach of the landslide numerical risk factor (LNRF) bivariate model was used in ensemble with linear multivariate regression (LMR) and boosted regression tree (BRT) models, coupled with radar remote sensing data and geographic information system (GIS), for landslide susceptibility mapping (LSM) in the Gorganroud watershed, Iran. Fifteen topographic, hydrological, geological and environmental conditioning factors and a landslide inventory (70%, or 298 landslides) were used in mapping. Phased array-type L-band synthetic aperture radar data were used to extract topographic parameters. Coefficients of tolerance and variance inflation factor were used to determine the coherence among conditioning factors. Data for the landslide inventory map were obtained from various resources, such as Iranian Landslide Working Party (ILWP), Forestry, Rangeland and Watershed Organisation (FRWO), extensive field surveys, interpretation of aerial photos and satellite images, and radar data. Of the total data, 30% were used to validate LSMs, using area under the curve (AUC), frequency ratio (FR) and seed cell area index (SCAI). Normalised difference vegetation index, land use/ land cover and slope degree in BRT model elevation, rainfall and distance from stream were found to be important factors and were given the highest weightage in modelling. Validation results using AUC showed that the ensemble LNRF-BRT and LNRFLMR models (AUC = 0.912 (91.2%) and 0.907 (90.7%), respectively) had high predictive accuracy than the LNRF model alone (AUC = 0.855 (85.5%)). The FR and SCAI analyses showed that all models divided the parameter classes with high precision. Overall, our novel approach of combining multivariate and machine learning methods with bivariate models, radar remote sensing data and GIS proved to be a powerful tool for landslide susceptibility mapping.