###### Abstract  
Under the influence of complex environmental conditions, the spatial heterogeneity of soil organic matter (SOM) is inevitable, and the relationship between SOM and visible and near-infrared (VNIR) spectra has the potential to be nonlinear. However, conventional VNIR-based methods for soil organic matter estimation cannot simultaneously consider the potential nonlinear relationship between the explanatory variables and predictors and the spatial heterogeneity of the relationship. Thus, the regional application of existing VNIR spectra-based SOM estimation methods is limited. This study combines the proposed partial least squares-based multivariate adaptive regression spline (PLS-MARS) method and a regional multi-variable associate rule mining and Rank-Kennard-Stone method (MVARC-R-KS) to construct a nonlinear prediction model to realize local optimality considering spatial heterogeneity. First, the MVARC-R-KS method is utilized to select representative samples and alleviate the sample global underrepresentation caused by spatial heterogeneity. Second, the PLS-MARS method is proposed to construct a nonlinear VNIR spectra-based estimation model with local optimization based on selected representative samples. PLS-MARS combined with the MVARC-R-KS method is illustrated and validated through a case study of Jianghan Plain in Hubei Province, China. Results showed that the proposed method far outweighs some available methods in terms of accuracy and robustness, suggesting the reliability of the proposed prediction model. s 2021 11 566 10 October 2020 4 January 2021 8 January 2021 8 January 2021 8 January 2  
multiple linear regression (MLR), partial least squares (PLS) regression, and geographic weighted regression (GWR) are typical linear prediction methods. These methods have been extensively used in many previous studies because they are easily accessible in most software packages [4]. A common characteristic of these prediction methods is that they assume, either explicitly or implicitly, a linear relationship among the analyzed data sets. A nonlinear relationship, however, is prevalent in practice between SOM content and VNIR spectra. Thus, a series of nonlinear prediction methods has been proposed to estimate nonlinear relationships, including, but not limited to, local weighted regression (LWR) [8; 9], artificial neural network [10; 11], and support vector machine (SVM) methods [12; 13]. Although these methods can efficiently determine nonlinear relationships in certain situations, they lead to models based on global optimization and disregard the spatial heterogeneity of the nonlinear relationship. This underlying principle will cause the accuracy of existing methods to easily change with different samples and have difficulty for application to other regions. Therefore, nonlinear estimation methods considering spatial heterogeneity are needed for SOM content prediction. The multivariate adaptive regression spline (MARS) method is a nonlinear prediction method that considers the effects of local difference. The MARS method is the expansion of splines, thereby providing immense flexibility by automatically determining the splines (i.e., the number of basis functions (BFs) and locations of knots) without human interference. This method has been successfully applied to various fields, such as geotechnical engineering and soil liquefaction assessment [14; 15; 16]. The MARS method has proven to be advantageous due to its adaptation to nonlinear relationships and to its adaptive model construction process, while considering the natural local difference in the training data sets. These characteristics of the MARS method indicate its considerable potential to represent the relationships between VNIR spectra and SOM content. Although tracked records of the application of the MARS model to soil property prediction are available in literature [11; 17], these records are not well received.