土壤有机质(Soil Organic Matter，SOM)是土壤的重要组成部分，在土壤多种功能中均扮演着核心角色，是土壤质量评估、土壤功能评价等诸多应用中的关键参数，同时也是土壤资源利用与管理的关键目标[1]，因此有效获取SOM的空间分布信息至关重要。近年来，遥感技术的快速发展和有效应用为区域SOM预测提供了新的途径和有利条件。SOM是土壤的重要染色物质[2]，其在可见光波段乃至短波红外波段都能够显著吸收入射光，SOM含量与土壤光谱反射率之间通常呈显著负相关，因此利用土壤光谱信息能够有效反映SOM含量的变化[3]。

通常利用室内测量的土壤光谱能够高精度预测SOM含量。但利用卫星影像反演预测土壤属性时，土壤自身条件(如土壤湿度和粗糙度)和土壤表面条件(如作物和作物残茬覆盖)等会显著影响土壤反射率，进而严重降低土壤属性的预测精度[4-5]。有研究表明，利用土壤光谱指数可以一定程度上抑制土壤水分、粗糙度和大气等因素的影响[6]。与直接利用卫星影像波段反射率相比，基于线性或非线性波段组合形成的光谱指数所构建的土壤属性预测模型精度通常较好。例如，杨邵文[7]在建立北京市密云区耕地SOM预测模型时，选用了Landsat 8卫星数据计算的差值光谱指数和比值光谱指数。赵瑞等[8]基于GF-5号卫星高光谱影像，通过构建多种光谱指数(包括土壤沙化指数、土壤退化指数、归一化亮度指数和土壤盐分指数)反演SOM，实现了对内蒙古与黑龙江交界区域的SOM预测。上述SOM遥感反演研究中，涉及到多种类型的光谱指数，表明光谱指数方法具有提升区域SOM预测精度的潜力，但这些研究构建的光谱指数没有考虑时相信息，均是利用单一影像不同波段间组合构建的光谱指数。

为进一步提高SOM反演精度，许多研究者引入了时相信息，通过将两景不同时期遥感影像波段进行组合，构建双时相光谱指数[9-10]。其中，刘焕军等[10]认为，不同SOM含量的土壤在相同外部条件(诸如土壤水分)影响下其反射光谱敏感程度不同，通过选择不同时期影像构建的相同波段组合的双时相光谱指数，具有提高有机质遥感反演的潜力；其利用同一年内多景裸土期MODIS影像，构建了包括差值指数等的3种双时相光谱指数进行SOM反演，获得了较高的SOM预测精度。Dou等[6]基于相似的理论假设，认为选择更长时期的遥感影像组合构建双时相光谱指数同样具有提高SOM预测的潜力；研究中，其选择与土壤采集时间相距超过10年的多年MODIS影像，构建了包括差值指数、比值指数等4种双时相光谱指数反演东北黑土区SOM，取得了较为理想的SOM预测结果。上述研究利用双时相光谱指数方法，针对东北黑土区SOM反演均取得了较好的预测结果。但与东北黑土区(一年一熟)相比，我国黄淮海地区(一年两熟)裸土期时间窗口更短，同时还需要考虑云等因素对影像质量的影响，因此在黄淮海地区进行土壤SOM反演可能需要更长时期的遥感影像。另外，在遥感影像获取日期跨度较大的情形下(比如5年间隔两景不同影像)，影响不同影像间土壤光谱差异的因素可能更为复杂，除采用相同波段外，还可以考虑采用不同影像不同波段构建双时相光谱指数。因此，本研究拟针对黄淮海典型县域，进一步探索双时相光谱指数在农田土壤SOM遥感反演中的潜力。

当运用双时相光谱指数建立SOM预测模型时，可能会出现光谱指数和土壤样品数量接近的情形。例如，在刘焕军等[10]的研究中，共计使用了61个土壤样品，双时相光谱指数高达49个；在Dou等[6]的研究中，共有56个双时相光谱指数，土壤样品为147个，其中仅一半样品用于建模。这些研究中的模型均属于小样本高维模型，在相关研究中其通过相关分析和逐步多元回归法来筛选变量，降低变量间共线性。此外，有研究表明，筛选有效的建模变量和选择更为先进的机器学习方法可以提高SOM的预测精度[11-12]。Wang等[13]提出可以将最小绝对收缩和选择算子(Least Absolute Shrinkage and Selection Operator，LASSO)算法用于变量选择，然后结合支持向量机(Support Vector Machine，SVM)算法建立目标对象的预测模型。该方法已被成功用于SOM含量估算、森林生物量估算和数字土壤制图等多个领域[11, 14-15]，并取得较为理想的结果。

因此，本研究以具有短裸土期特点的黄淮海平原典型县域封丘县为研究区，利用Landsat 8影像，通过构建多种类型的双时相光谱指数，并结合LASSO筛选变量方法与SVM算法建立SOM预测模型，探讨了基于两景影像波段组合构建的双时相光谱指数在SOM预测中的表现，以期提高SOM预测精度，并实现区域SOM空间分布预测。

研究区选择在河南省封丘县，该县地处黄淮海平原中部，其东、南两面毗邻黄河，地理坐标为114°14' E ~ 114°43'E，34°53' N ~ 35°14'N，面积约1 220 km2，区域内主要土壤类型为潮土[16]。该地区属于典型的季风气候，年平均气温为15℃，年平均降水量在500 ~ 600 mm。封丘县是重要的粮食生产县，区域内普遍采用小麦–玉米一年两熟的种植制度[17]，夏玉米一般在6月小麦收获后种植，冬小麦在同年10月玉米收获后种植。

通常，研究区内小麦收获后采用秸秆还田的保护性耕作，此时土壤表面留有较多的小麦秸秆和残茬，不利于土壤属性信息遥感监测，而在玉米收获之后通常会进行土壤翻耕。因此，研究区内10—11月内有一段土壤相对裸露、受作物和作物残茬影响相对较小的时期，是土壤属性遥感反演的潜在时间窗口。

土壤样品于2014年6月小麦收获后采集。采样前，在室内按照3 km×3 km的规则网格在研究区布设土样采集点位。野外土壤样品采集时，根据交通通达性等因素进行一定的点位调整，研究共计采集表层(0 ~ 20 cm)土壤样品117件。采样同时使用亚米级高精度GPS设备记录采样点坐标，其点位分布如图 1所示。

采样时，应尽可能地远离周边道路、防护林等非耕地区域，以避免采样点周边道路和田间防护林所导致的混合像元问题。土壤样品采集采用五点采样法，点位间距离设定为5 ~ 10 m，5个不同位置采集的土样混合后作为该采样点的代表性土样。采集的土样运至实验室，经风干、研磨、过筛处理后，采用重铬酸钾容量法测定SOM含量。

选择Landsat 8卫星平台OLI传感器影像作为影像数据源，其空间分辨率为30 m。选择前7个光学范围内的波段用于SOM含量反演，波段1至波段7的光谱范围分别为：0.44 ~ 0.45、0.45 ~ 0.51、0.53 ~ 0.59、0.64 ~ 0.67、0.85 ~ 0.88、1.57 ~ 1.65、2.11 ~ 2.29 µm。获取的影像为Landsat 8地表反射率产品，已经经过辐射定标、大气处理和几何校正处理。

研究区的潜在裸土窗口期为10—11月，该时间窗口内受土壤表面作物和作物残茬的影响相对较小，能够较好地反映土壤光谱信息。此外，由于SOM含量随时间变化相对较小，并且为了获得更多有效的影像，首选初选了采样时间近5年内(2013—2017年间)10—11月之间云量小于20% 的影像，经过筛选能够覆盖全部土壤采样点的Landsat 8影像共有7景(条带号：123，行编号：036)，其时间和云量信息如表 1所示。对所有影像人工检查后发现，部分影像仍然受云的较大影响，因此本研究进一步选取云量小于1% 的高质量影像来预测SOM，其获取日期分别为2014年10月6日和2017年10月30日。为了尽可能去除云的影响，研究中还针对选取的两景影像利用Landsat 8影像数据中的QA_PIXEL波段(质量评价波段)，对云和云阴影进行掩膜处理，其中2014年10月6日获取的影像有1个采样点被云覆盖，2017年10月30日没有采样点被云覆盖，因此，以去除该采样点后共计116个采样点用于数据分析。

已有研究表明，与直接使用影像波段反射率相比，光谱指数能够提高土壤遥感反演精度[6-7]。本研究中，采用比值光谱指数(Ratio Spectral Index，RI)、差值光谱指数(Difference Spectral Index，DI)、归一化光谱指数(Normalized Spectral Index，NDI)以及优化光谱指数(Optimized Spectral Index，OSI)[18]进行SOM预测。构建的光谱指数包括单时相光谱指数和双时相光谱指数。单时相光谱指数指的是同一景影像不同波段之间构建的指数，双时相光谱指数是指两景不同影像波段之间构建的指数。计算公式如表 2所示。

研究中共有2014年10月6日和2017年10月30日获取的两景有效遥感影像，每一景影像可以得到单时相光谱指数168个，其中RI、DI、NDI以及OSI指数各42个。两景遥感影像共计得到双时相光谱指数196个，各类光谱指数均为49个。

变量选择是模型构建的一个重要方面，特别是当自变量个数较多时，通过变量筛选能够有效降低变量个数，提高模型构建效率及稳定性。本研究选择LASSO算法进行变量筛选，并与未经变量筛选的模型结果进行对比分析。LASSO算法在化学计量学相关研究中已有应用，其在最小化模型残差平方和的基础上增加了惩罚函数，将模型回归系数的绝对值约束在一定的阈值，能够有效改善模型精度和稳定性[19]。其通过优化目标函数将相关性较小的变量系数压缩至0，通过删除相关系数为0的变量，达到筛选变量的目的，可以有效防止过拟合问题。LASSO变量筛选通过使用R语言中的“glmnet”包实现。

为探究双时相光谱指数在SOM反演中的表现，本研究分别针对影像波段、单时相光谱指数以及双时相光谱指数构建模型。针对影像波段建模以及针对单时相光谱指数建模时，分别使用单景影像以及两景影像相关波段和指数构建模型；针对双时相光谱指数构建模型时，直接使用基于两景影像构建的指数建立模型。同时为了对比分析变量筛选在模型构建中的作用，除直接使用影像波段或光谱指数构建模型外，还构建了经LASSO算法变量筛选后的模型。

本研究中采用的建模算法为SVM，SVM算法是一种经典的机器学习算法，在许多应用场景中均有良好的表现[20]。该算法使用核函数将数据投影到高维空间，以找到用于对不同组进行分类或建立最佳回归模型的超平面[21]。SVM模型构建采用R语言中的“caret”包实现，选用线性核函数，并采用留一交叉验证对模型参数进行优化。模型精度评价采用留一交叉验证的结果，根据决定系数(R2)和均方根误差(RMSE) 指标评估SOM反演精度，计算公式如下：

式中：Qi为观测值；Pi为实测值；n为样本量。

研究区SOM含量范围为8.60 ~ 28.40 g/kg，中值为16.50 g/kg，均值为16.40 g/kg，标准差为3.01 g/kg，变异系数为18.36%，属于中等程度的变异[18]。SOM的峰度系数为–0.65，偏度系数为0.41，其含量分布符合正态分布(图 2)，因此建模时未对测量结果进行转换处理。同时，为去除异常值对建模的影响，选择格拉布斯检验(Grubbs Test)方法对SOM结果进行异常值检测，在P < 0.001的显著性水平下未发现异常值的存在[22]，因此本研究直接使用去除云影响之后的所有116个样本进行建模分析。

分别以单期影像和两期影像的波段反射率构建模型，SOM预测的交叉验证结果如表 3所示。单期影像中，基于2014年影像的SOM预测精度较低，R2仅为0.05，RMSE为2.88 g/kg；与之相比，基于2017年影像的SOM预测精度相对较高，R2=0.30，RMSE=2.47 g/kg。两期影像同时使用时，与基于2017年单期影像相比，SOM预测精度并未改善(R2=0.27，RMSE=2.52 g/kg)。经筛选变量后，基于单期影像和两期影像的SOM预测精度均有了一定程度的提升，其中两期影像同时使用时SOM预测R2提高到0.32，RMSE降低为2.44 g/kg，但基于2014年影像的SOM预测精度仍然较低，R2=0.09，RMSE=2.82 g/kg。

针对单时相光谱指数，分别利用单期影像和两期影像构建模型，SOM预测的交叉验证结果如表 4所示。对于2014年影像，所有光谱指数的SOM预测精度均相对较低，R2=0.17，RMSE=3.02 g/kg，经变量筛选后模型预测R2提升至0.21，RMSE降低为2.63 g/kg。基于2017年影像的SOM预测精度优于基于2014年影像的预测结果，且变量筛选后模型的预测精度也有所提升，变量筛选前、后的预测精度分别为：R2= 0.35，RMSE=2.38 g/kg和R2=0.37，RMSE=2.34 g/kg。当两景影像同时使用时，特别是经变量筛选后，模型的SOM预测精度显著高于基于单景影像的预测结果，R2=0.45，RMSE=2.20 g/kg。

双时相光谱指数经变量筛选前后的SOM预测模型评价结果如表 5所示。与使用波段或者单时相光谱指数的预测结果相比，有机质的预测精度有了较大程度的提升，未经变量筛选时，SOM的预测R2可达0.45，RMSE为2.19 g/kg，经变量筛选后SOM的预测精度更高，R2为0.53，RMSE降低至2.01 g/kg。双时相光谱指数模型的有机质预测散点图如图 3所示，可以看出，整体上有机质预测值与实测值较为接近，特别是经变量筛选后模型的SOM预测值更接近1∶1线，表现出更优的预测精度。

高精度的SOM空间分布信息对土壤资源高效管理具有重要意义。当前许多研究基于土壤–景观关系，利用数字土壤制图技术进行土壤空间分布图件的制作，在相关研究中土壤成土环境因子特别是地形相关变量被广泛采用[23]。然而，地形成土环境变量在平缓地区土壤制图中往往面临较大困难。研究区位于黄淮海平原，其海拔范围在65 ~ 72.5 m [24]，变异相对较小，因此难以利用地形相关变量进行土壤属性空间分布制图。在与本研究区类似的平缓地区，遥感可以有效获取土壤信息，直接反演土壤属性空间变异。双时相光谱指数的SOM预测结果展示出了较高的预测精度(表 5)，具有制作SOM空间分布图的能力。

利用基于变量筛选后的双时相光谱指数构建SOM预测模型对封丘县内耕地土壤进行SOM预测，结果如图 4所示。结果表明，封丘县内SOM含量分布差异较大，SOM含量较高的区域主要分布在封丘县南部与黄河平行一线，以及东北部部分区域，而紧靠黄河沿线的SOM含量较低。该预测结果空间上基本与实测值分布一致，统计后发现封丘县内SOM预测值的平均结果与SOM实测值的平均值相同，均为16.40 g/kg，表明基于双时相光谱指数的SOM空间预测结果较好，相关方法可以用作SOM制图。

利用光学遥感影像进行土壤属性定量遥感反演时容易受地表植被和土壤水分等因素的影响。本研究中选择的遥感影像在潜在裸土期时间窗口内(10月)，但直接利用影像波段进行SOM预测的精度较低，最高精度模型R2=0.32，RMSE=2.44 g/kg，这可能是土壤表面覆盖物或者土壤自身状态的影响所致。2014年10月6日获取的遥感影像，其SOM预测精度较低，R2仅为0.05，RMSE高达2.88 g/kg，而2017年10月30日获取的遥感影像预测精度相对较高，两景影像精度差异的原因很可能是由于10月6日获取的影像受到了玉米秸秆的影响，该影像获取时可能部分区域秸秆尚未收获或者土壤未开始翻耕，而10月30日时土壤均已翻耕完毕且冬小麦尚未影响到土壤光谱。两景影像的植被指数差异印证了这一推测，10月6日影像的平均归一化光谱指数为0.13，而10月30日的对应值为0.11。因此，对于黄淮海平原一年两熟区，尽管区域内普遍实行抢收抢种的农田管理模式，但小农户种植管理比重大，因此，某一时刻的土壤表面状况也可能存在较大差异，遥感影像的时相选择显得至关重要。

已有研究发现通过构建光谱指数，能够一定程度上消除外部因素对土壤遥感精度的影响，提高目标土壤属性反演精度[6, 9-10]。本研究中，与直接使用波段反射率构建的模型相比，基于单时相光谱指数的SOM预测精度明显更高，当两景影像构建的指数共同使用并经变量筛选后，SOM的预测精度达到一个相对较高的水平(R2 = 0.45，RMSE=2.20 g/kg)。但与徐夕博[25]在山东省潍坊市北部平原地区利用Landsat 8影像构建单时相光谱指数模型预测SOM的结果相比，本研究中的预测精度相对较低，其模型最高精度R2=0.65，RMSE=2.50 g/kg。邱壑等[26]在福建省三明市基于Hyperion高光谱遥感数据进行SOM预测时，同样发现与基于波段反射率建模相比，基于单时相光谱指数建模可以有效提高反演精度，其最高模型精度R2=0.55，RMSE=8.21 g/kg，而基于影像波段建模最高模型精度R2=0.17，RMSE=11.14 g/kg。该结果与本研究中的预测精度较为接近，这可能是由于尽管作者使用了高光谱卫星影像，但其研究区与黄淮海平原相比，耕地土壤更易受到作物覆盖和土壤自身状态的影响。本研究发现，与仅利用一景影像的单时相光谱指数相比，两景影像的单时相光谱指数共同使用时SOM预测精度得到了提升(表 4)，该结果与Wang等[27]的研究类似，其在黑龙江省绥化市基于Sentinel-2影像预测SOM的研究中发现，同时使用多景影像的单时相光谱指数构建的模型精度(R2=0.49，RMSE= 3.53 g/kg)优于单景单时相光谱指数的结果(R2=0.31，RMSE=3.84 g/kg)。单景影像仅能记录地表的瞬时状态，而多景影像可以反映地表的动态变化过程，因此，更多时相信息的加入可以改善SOM预测精度。

对比影像波段模型和单时相光谱指数模型，尽管使用了相同的两景影像，双时相光谱指数SOM的反演精度更高，最高模型精度R2=0.53，RMSE=2.01 g/kg。因此，遥感影像时相信息的利用方式也可能对土壤属性的预测精度产生影响。实际上，双时相光谱指数利用的是目标土壤属性对外部因素(如土壤水分、土壤粗糙度和秸秆覆盖度等)的响应差异来间接反映土壤属性的差异，而单时相光谱指数利用的是土壤属性在不同波段之间的差异。刘焕军等[10]在我国东北黑土区，选择了同一年内多景裸土期内的MODIS影像构建双时相光谱指数，SOM预测R2可达0.72，RMSE为4.42 g/kg。在黑土区田块尺度，潘越[28]利用同一年份的Landsat 8和Sentinel-2遥感影像，构建双时相光谱指数对SOM预测同样获得了较高的精度(R2=0.80，RMSE=8.98 g/kg)。黑土区为一年一熟区且区域内耕地多由农场统一管理，因此，一定区域内的农田耕作方式比较一致。但与我国东北黑土区不同，本研究区所在的黄淮海平原耕地普遍采用一年两熟的种植模式，两季作物之间的裸土窗口期短，同一年内获取多景裸土期高质量影像难度较大。当选择构建双时相光谱指数的遥感影像时，影像的获取时间一方面受到影像可获取性、云等干扰因素以及裸土期等因素的影响，另一方面还可能会受到土壤属性随时间变化的影响。但对于SOM，其短期内含量变化较小，影像选择时间可以不局限于采样当年。针对SOM预测，本研究中影像选择时间设置为5年，但去除受云等干扰因素影响的影像后仅有两景高质量影像，基于这两景非同年获取的影像构建的双时相光谱指数同样得到了较高的SOM预测结果。实际上，即便在裸土期较长的东北地区，使用更长时间间隔的影像构建双时相光谱指数有可能得到更好的预测结果。例如，Dou等[6]在针对黑土区SOM遥感反演预测中，基于2002年、2009年和2014年3年的MODIS影像构建的双时相光谱指数，使用2002年影像时SOM预测结果最优(R2=0.76，RMSE=7.66 g/kg)。因此，针对SOM等随时间变化较慢的土壤属性时，可以选择长时序的遥感影像进行反演，但对于随时间变化较快的土壤属性，遥感影像的时相选择是否对预测结果产生影响还需要进一步研究。

利用遥感影像进行土壤属性预测反演制图时，受制于成本限制，通常建模样品数量较少。本研究中，实际利用116个土壤样品构建模型，与光谱指数的数量相比，建模样本相对较少，属于小样本高维模型。这类模型构建时需要进行变量筛选，以减少变量共线性的影响，降低模型过拟合概率，提高模型稳定性和精度[29]。本研究利用LASSO算法筛选变量，然后利用SVM机器学习方法构建模型，取得了较高的SOM反演精度，所构建模型预测的SOM空间分布格局与真实分布较为吻合。在今后研究中，可以尝试使用其他筛选变量算法，例如遗传算法和竞争性自适应重加权算法，进一步分析变量筛选对模型精度的影响。此外，本研究针对土壤SOM预测反演仅使用了Landsat 8影像，在今后研究中可以尝试更高时空分辨率的影像(如Sentinel-2影像)，进一步研究双时相光谱指数的应用潜力，实现区域SOM高精度反演和高空间分辨率预测制图的结合。

1) 与直接使用影像波段反射率或者单时相光谱指数相比，利用双时相光谱指数可以更好利用影像时相信息的优势，SOM预测精度更高。针对影像波段和光谱指数，经LASSO变量筛选之后，SOM预测模型精度更优。

2) 在与研究区类似的地形平缓区域，利用双时相光谱指数构建的遥感预测模型，可以有效用于区域SOM预测与制图。对于SOM等随时间变化相对缓慢的土壤属性，遥感影像可以选择不同于采样年份内的(准)裸土期影像。