摘要
采集新疆乌鲁木齐市 110 个风干的土壤样品，将筛选出的土壤样品分别过 2.00、0.50、0.25 mm 的筛，在室内进行反射率光谱数据的采集，对采集的光谱数据（R）与土壤养分含量数据剔除异常值，随后对光谱数据进行断点校正以及 10 种光谱反射率变换，筛选出与有效磷含量相关性最高的光谱变换为对数的一阶微分（lgR）′变换。在此基础上利用竞争性自适应加权算法（CARS）对经过（lgR）′处理后的数据进行特征波段的筛选，用偏最小二乘回归（PLSR）、BP 神经网络、随机森林建立土壤有效磷含量预测的高光谱分析模型。模型评价指标采用决定系数、均方根误差、相对分析误差、平均绝对误差。结果显示：土壤原始光谱特征在各个波段与有效磷含量相关性都较差，不同形式的光谱数据变换均能够提高光谱反射率与有效磷含量的相关性；对比不同粒径处理的模型预测精度，过筛粒径越小对有效磷含量的估测精度越高，3 种方法的最优拟合模型都是过 0.25 mm 筛的处理；对过 2.00、0.50、0.25 mm 筛的土壤样品有效磷含量分别采用 3 种方法建模，均以 PLSR 建立的模型预测精度最为突出，其模型决定系数、均方根误差、相对分析误差、平均绝对误差分别为 0.76、10.96、2.07、8.69；0.77、12.52、 2.00、9.74；0.77、10.90、2.13、8.81。因此，对过 0.25 mm 筛处理的土壤有效磷，利用 CARS 算法筛选出光谱特征波段结合 PLSR，能够较好的估测乌鲁木齐市土壤中有效磷的含量，为室内高光谱估测土壤养分含量提供理论基础。
Abstract
110 air-dried soil samples were collected from Urumqi，Xinjiang. The selected soil samples were sieved at 2.00， 0.50 and 0.25 mm，respectively. The reflectance spectral data were collected indoors，and the abnormal values were removed from the collected spectral data（R）and soil nutrient content data. Then the breakpoint correction and ten spectral reflectance transformations were carried out on the spectral data，and the log-first derivative（lgR）' with the highest correlation with the available phosphorus content was selected and transformed. On this basis，the competitive adaptive weighting algorithm （CARS）was used to screen the characteristic bands of the processed（lgR）′data，and the hyperspectral analysis model of soil available phosphorus content prediction was established by partial least square regression（PLSR），BP neural network and random forest. Coefficient of determination，root mean square error，relative analysis error and mean absolute error were used as evaluation indexes of the model. The results showed that the correlation between the original spectral characteristics and available phosphorus content was poor in each band，and the correlation between spectral reflectance and available phosphorus content could be improved by different spectral data transformation. Comparing the prediction accuracy of the models with different particle size treatments，the smaller the size of the sieve，the higher the estimation accuracy of the available phosphorus content. The optimal fitting models of the three methods were all over 0.25 mm sieve treatment. For the available phosphorus of 2.00，0.50 and 0.25 mm treatments，PLSR was used to establish the model with the most prominent prediction accuracy. The model determination coefficient，root mean square error，relative analysis error and average absolute error of the three methods were 0.76，10.96，2.07，8.69；0.77，12.52，2.00，9.74；0.77，10.90，2.13，8.81， respectively. Therefore，the characteristic bands selected by CARS algorithm and PLSR could better estimate the content of available phosphorus in soil of Urumqi.
Keywords
土壤有效磷含量是农田土地评价的指标，也是农作物施肥的基本指标，快速准确测量土壤有效磷含量是土壤信息化管理和资源评价的前提条件[1]，传统的土壤有效磷含量测定方法采用实验室化学法，此方法虽然相对准确，但通常比较费时、费工、有害或有污染。近年来，随着高光谱技术的快速发展，为土壤快速、准确检测提供了可能[2-3]。高光谱遥感波段多、分辨率高，具有定量获取土壤有效磷含量的潜力[4-6]，目前利用高光谱遥感技术对土壤有效磷含量进行快速检测，已成为土壤有效磷含量检测的一种新方法[7]。
目前的研究普遍认为，对光谱数据采用不同的预处理和模型能够很好地构建出反演模型，并提高模型的预测精度，如对原始光谱数据进行倒数、对数、一阶微分、倒数的对数、连续统去除等变换[8-14]，采用多元线性回归、偏最小二乘回归 （PLSR）、BP 神经网络，随机森林、支持向量机等建模方法[15-20]。Bogrekci 等[21]利用可见光-近红外光谱对不同类型土壤风干样本和鲜样的磷含量进行了预测，应用偏最小二乘法来预测其均方根偏差分别为 9.4% 和 12.9%；Mouazen 等[22]通过偏最小二乘法建立了土壤有效磷含量的预测模型，发现利用可见光-近红外波段预测有效磷是可行的；Dhawale 等[23]利用 PLSR 方法，将 86 个土壤样品的 3 次重复光谱的平均值与 Mehlich-3 可提取磷测定值相关联。光谱校正模型的留一交叉验证结果表明，模型预测结果与参考测量值之间存在线性关系，决定系数（R2）为 0.86，预测标准误差为 31 mg/kg；Sarathjith 等[24]在 3500～2500 nm 波段范围，利用离散小波和支持向量回归构建模型来预测有效磷，其有效磷的相对分析误差（RPD）为 2.27；周鼎浩等[25]采用 PLSR 结合留一法交叉验证，建立水稻土全磷光谱反演模型，R2 为 0.85；薛利红等[26]研究发现土地利用方式对全磷反演影响较大，模型中 PLSR 分析精度较高且稳定；高会等[1]采用相关分析方法得出土壤有效磷含量的敏感波段为 711 nm，利用一元线性回归方法，基于该波段的估算模型为最佳估算有效磷含量模型，该模型的拟合优度 R2 达 0.8221，验证 R2 达 0.9591。
土壤的光谱特征是土壤理化特性的综合反映，而土壤颗粒粒径（比表面积）对土壤光谱反射率有着显著的影响[27]，前人的研究大都针对同一土壤粒径、不同光谱变换与建模的方法开展，但是土壤的光谱反射率与土壤的颗粒大小有着密切的联系，不同粒径的土壤光谱反射率存在差异。因此，本研究以新疆乌鲁木齐市的土壤为研究对象，分析不同过筛粒径的土壤光谱特征与有效磷含量的相关关系，筛选出最优光谱变化形式，分别建立有效磷含量的 3 种估测模型，确定有效磷含量的最优估测模型，为室内快速、高效测定土壤有效磷含量提供技术支持。
1 材料与方法
1.1 研究区域概况
新疆维吾尔自治区乌鲁木齐市位于亚欧大陆腹地，地处北天山北坡，准噶尔盆地南缘。其地理坐标北纬 42°45′32″～44°08′00″，东经 86°37′33″～88°58′24″。乌鲁木齐属中温带大陆性干旱气候，春秋两季较短，冬夏两季较长，昼夜温差大。年平均降水量为 194 mm。乌鲁木齐市三面环山，地势东南高、西北低，自然坡度 12%～15%，海拔 680～920 m。本研究选取沙依巴克区、头屯河区、新市区、高新区、达坂城区、米东区、水磨沟区、乌鲁木齐县均属于乌鲁木齐市。研究区主要以玉米、水稻、向日葵、红薯为主，土壤质地主要以砂壤土和黏土为主。
1.2 土壤采集与预处理
土壤样品采集时间为 2022 年，采集地点位于新疆维吾尔自治区乌鲁木齐市（图1），共采集 110 份土壤样本。样本的采集用网格布点法，样本采集前，按照 666.67 hm2 一个样点均匀布设，采样深度 0～20 cm，将采集好的土壤样品混合均匀放入密封袋中，贴好标签，标注采样时间和地点等信息，单个土样采集重量约为 500 g。土样带回室内后，去除杂草、砾石及表面杂质后自然风干。将风干后的土样进行研磨，分别过 2.00、0.50、0.25 mm 筛，混匀后分成两份，一份用于光谱数据的测定，另一份用于实验室土壤有效磷含量的测定。实验室检测土壤速效养分含量时，有效磷的检测使用碳酸氢钠浸提剂：称取 42.0 g 碳酸氢钠溶于约 950 mL 水中，用氢氧化钠溶液调节 pH 至 8.5，用水稀释至 1 L，贮存于聚乙烯瓶或玻璃瓶中备用，如贮存期超过 20 d，使用时必须检查并校准 pH。
图1研究区地理位置图
1.3 光谱数据的采集
光谱的采集使用美国 ASD 公司的 FieldSpec ProHR NG 型光谱仪进行土样光谱数据的测定，测定波长范围为 350～2500 nm，光谱分辨率在 350～1000 nm 为 3 nm，在 1000～2500 nm 为 6 nm，样品采样间隔在 350～1000 nm 为 1.4 nm，在 1001～2500 nm 为 1.1 nm，重采样间隔为 1 nm。光谱测定前采用标准白板对光谱仪进行校准和调整。为了提高光谱测定的准确性，光谱采集前需要对仪器进行 30 min 左右的预热，采集时将土样放于直径 18 cm、深 2.5 cm 的黑色盛样皿，将 ASD 光纤探头垂直置于土壤样品表面，采用 25°视场角探头，在光谱测定过程中转动器皿 3 次，每次转动 90°，各个土壤样本在每个角度分别采集 10 条光谱反射率数据，单个样本共采集 30 条反射率数据。
1.4 光谱数据的预处理
1.4.1 断点校正
由于受到仪器本身影响，光谱反射率曲线在 350～399 和 2401～2500 nm 这2个谱段反射率出现剧烈变换的锯齿现象，噪声影响较大，为了提高数据的准确性，故去掉这两段光谱。本研究采用 400～2400 nm 范围内的光谱数据进行土壤养分信息的分析。光谱探测仪器在 300～1000 和 1000～2500 nm 波长使用探测元件不同，采集的光谱曲线在 1000 nm 两侧有个跳跃点，为了消除这一误差，对光谱曲线进行断点校正，断点校正利用 FieldSpec ProHR NG 光谱仪自带的软件 ViewSpecPro 进行处理，校正前后的结果如图2所示。
图2断点校正前后对比
1.4.2 异常值剔除
土壤样品在采集和室内养分含量测定的过程中，容易受到外界环境的影响，从而导致测定的土壤样品含量值出现异常。为了保证所建立的模型有较高的精度，需要在模型建立前对土壤样品含量进行异常值剔除，本研究采用箱线图与蒙特卡罗离群点检测 2 种方法进行异常值的剔除。
1.4.2.1 箱线图离群点检测
箱线图又称盒子图和箱形图，箱线图主要由 6 个统计量：最小值、第一四分位数、中位数、第三四分位数、最大值和异常值来描述一组数据的分布情况。第一四分位数是该样本中所有数据从小到大排列后第 25% 的数字为 Q1，第三四分位数等于该样本中所有数值由小到大排列后第 75% 的数字为 Q3，四分位范围为Q3-Q1=∆Q，则在区间（Q1-1.5ΔQ，Q3+1.5ΔQ） 范围外的点就被称为离群值（如图3a）。采用箱线图法剔除有效磷异常值样本数 5 个，剩余样本 105 个。
1.4.2.2 蒙特卡洛离群点检测
蒙特卡洛离群点检测是一种基于统计学原理的离群点检测方法。它通过模拟数据的随机性来判断某个数据点是否为离群点。蒙特卡洛离群点检测的基本思想是如果某个数据点是离群点，那么它与其他数据点之间的差异应该显著大于与随机数据点之间的差异。因此，通过计算数据点与随机数据点之间的差异，并将其与差异的分布进行比较，可以判断数据点是否为离群点。图3b中平均值高或标准差高的样本可能是异常值，应在建模前将其剔除，本研究利用蒙特卡洛离群点检测剔除异常值 8 个，剩余样本 97 个。
图3有效磷异常值剔除
注：图中数据为样品编号。
1.5 光谱反射率变换
为了突出反应光谱的重要信息，寻找用于建模的最好的光谱变换形式，本研究采用的光谱反射率变换形式有倒数（1/R）、对数（lgR）、倒数的对数[lg（1/R）]、倒数的一阶微分[（1/R）']、对数的一阶微分[（lgR）']、倒数的对数一阶微分｛[lg （1/R）]'｝、多元散射校正（MSC）、连续统去除 （CR）、一阶微分（R'）、二阶微分（R＂）10 种光谱变化形式，光谱变换采用 Matlab2019 完成。
1.6 相关性分析
对土壤有效磷含量和原始光谱反射率以及经过倒数、对数、倒数的对数、倒数的一阶微分、对数的一阶微分、倒数的对数一阶微分、多元散射校正、连续统去除、一阶微分、二阶微分 10 种光谱变化后的光谱数据进行相关性分析。利用皮尔逊相关系数计算不同波长的光谱反射率与土壤样品元素含量之间的相关系数。公式如下：
上述公式中 Q 代表某一波段上土壤元素含量与光谱反射率的相关系数，n 代表样本数，xi 为 i 样本在某一波段的上的反射率值，为光谱平均值，yi 为 i 样本的土壤养分含量，为元素含量的平均值。 Q 绝对值越接近于 1，相关性越强。
1.7 特征波段筛选
本研究采用竞争性自适应重加权算法对各养分元素的光谱波段进行特征降维，得到相互独立的特征波段变量用于建立估测模型。作为一种新的变量选择算法，竞争自适应加权算法（CARS）是由 Li 等[28]在 2009 年提出的，该方法与偏最小二乘法相结合，变量选择主要包括 4 个步骤：
（1）蒙特卡罗采样：每次釆样都需从建模样本集中随机抽出特定比率的样本建立 PLSR 模型。
（2）对于值过小的波段，采用指数衰减函数强制消除。
（3）经过 N 次采样后，可以得到 N 个变量子集，对比分析采样过程中计算出的变量子集的交互验证均方根误差（RMSE），当计算得出变量子集的误差相对较小时，进行统计并作为最优变量子集。
（4）根据以上步骤对土壤养分元素波段进行特征选择提取，基于 Matlab2019 编写 CARS 算法计算程序，计算过程如图4所示。
图4为土壤有效磷特征波段的选取过程图，图中反应了样品采样次数与波段数量变化的关系，以及样品采样次数与交叉验证 RMSE 的关系。从图4a中可以看出，由于样本的采集次数增加，使得指数衰减函数算法保留的波段越来越少，且减少的速度越来越缓慢；从图4b可以看出，随着样品采样次数的增加，RMSE 也随着变化，变化趋势整体表现为先减少后增加；当运行到第 29 次时 RMSE 最小，说明此时建模效果最好，在此之前的 28 次采样计算，说明在不断剔除共线性较强的波段，进而在不断提高模型的精度，31 次以后再剔除波段将会降低模型的精度。由此说明经过 29 次运行模型所筛选出的波段，含有最少的赘余信息。此时所筛选出的 29 个光谱波段就是有效磷的特征波段，所选波段：401，405，417，451，455，458，471，477， 500，501，510，514，515，531，961，974，1001， 1042，1069，1894，1895，2211，2212，2213， 2282，2325，2347，2353，2373。
图4特征波段的筛选
1.8 模型的建立与验证
模型精度采用 R2、RMSE、RPD、平均绝对误差（MAE）4 个参数衡量。R2 反映模型建立和验证的稳定性，R2 越接近于 1，说明模型的稳定性越好、拟合程度越高。RMSE 用来检验模型的预报能力，RMSE 越小模型估算能力越好。RPD 是样本标准差与 RMSE 的比值，判断模型的预测能力，当 RPD<1.4 时，模型无法对样品进行预测；当1.4 ≤ RPD<2 时，认为模型效果一般，可以用来对样品进行粗略评估；当 RPD ≥ 2 时，模型具有极好的预测能力[29-30]。MAE 表示模型预测值与样本真实值之间距离的平均值，MAE 值越小说明样本预测值与真实值之间的误差越小。
式中：是预测值；y、yi 均是实测值；是实测值的均值；n 是样本数量，用 i=1，2，…，n 表示。
表1为有效磷建模集与验证集的划分。
2 结果与分析
2.1 不同粒径土壤原始光谱特征与有效磷含量相关性
图5表示不同粒径土壤样品光谱特征与有效磷含量的相关性，土壤原始光谱反射率在各个波段与有效磷含量的相关系数均不高，从其变化趋势来看，不同粒径处理的土壤样品与有效磷的相关性变化相似，波长在 400～646 nm 内表现为正相关，在 647～1191 nm 内相关性为负相关，随后又变为正相关。不同粒径土壤样品在 756～780 nm 之间出现谷值，相关系数在-0.055～-0.056 之间；而后在 1402～1413 nm 处与土壤有效磷含量的相关性出现峰值，相关系数为 0.082～0.086；随着波长的增加，相关性增高，在 1904～1921 nm 内形成峰值，相关系数约为 0.184～0.195；随后在 2100 nm 附近形成谷值，相关系数在 0.15 左右；从相关性来看，不同粒径组成的土壤原始光谱与有效磷含量相关性均不高，其中 2.00 mm 粒径处理相关性最低， 0.25 mm 粒径处理的相关性最高，试验结果表明，土壤粒径较小的处理在相关性上优于粒径较大的处理。
图5不同粒径土壤原始光谱特征与有效磷含量相关性
2.2 土壤有效磷含量与光谱反射率的相关性分析
通过 10 种变换对比发现，所有的微分变换形式后与土壤养分含量的相关性曲线，较前面几种变换后的曲线，曲线的波动更为明显、波峰与波谷更多。有效磷的所有变换中对数的一阶微分效果最好，相关系数值介于-0.476～0.492 之间 （图6）。
图6不同光谱变换的有效磷相关性分析
2.3 不同粒径的土壤有效磷估测模型的构建
经过 10 种光谱变换筛选出与有效磷含量相关性最高的一种光谱变换，应用 CARS 竞争性自适应加权算法选择光谱的特征波段作为模型的自变量，以土壤的有效磷含量作为因变量，建立土壤有效磷含量的 3 种估测模型。
分别对 2.00、0.50、0.25 mm 3 种粒径的土壤有效磷，建立 PLSR 模型、BP 神经网络模型、随机森林模型（表2）。对比不同粒径的土壤有效磷所建立的估测模型发现，基于 PLSR 所建立的 3 种粒径的估测模型，模型的 RPD 值分别为 2.07、2.00、2.13，均达到 2.00 以上，能够较好的估测土壤有效磷含量；对比 MAE 发现，3 种粒径的土壤有效磷建立的 PLSR 模型，预测值与实测值的误差由大到小依次为 MAE0.50>MAE0.25>MAE2.00；对比基于 BP 神经网络所建立的 3 种粒径的土壤有效磷含量的估测模型，模型的 RPD 值分别为 1.92、1.98、1.86，均未达到 2.00 以上，只能粗略的估测土壤有效磷含量；对比 MAE 发现，3 种粒径的土壤有效磷建立的模型，预测值与实测值的误差由大到小依次为 MAE0.25>MAE2.00>MAE0.50；对比基于随机森林所建立的 3 种粒径的土壤有效磷含量的估测模型，模型的 RPD 值分别为 1.90、1.90、1.94，均未达到 2.00 以上，只能粗略地估测土壤有效磷含量；对比 MAE 发现，3 种粒径的土壤有效磷建立的模型，预测值与实测值的误差由大到小依次 MAE0.50>MAE2.00>MAE0.25； 综合比较了 3 种粒径的土壤有效磷建立的 PLSR 模型、BP 神经网络模型、随机森林回归模型发现，基于 PLSR 所建立的估测模型能够较好地估测土壤中的有效磷含量，且 3 种粒径中 0.25 mm 粒径的有效磷所建立的估测模型精度最高；而基于 BP 神经网络与随机森林所建立的有效磷含量的估测模型只能对有效磷含量进行粗略的估测。为了更直观地分析 PLSR 模型预测的效果，将模型的预测值和实测值绘制成散点图，如图7所示。图中的虚线表示预测值的趋势线，实线代表 y=x；模型的精度越高，预测值分布越均匀，整体越靠近 y=x 这条直线。
表2不同粒径土壤有效磷含量的模型构建与验证
注：PLSR、BP、RF 分别表示偏最小二乘回归、BP 神经网络、随机森林。
图72.00 mm 土壤有效磷含量的偏最小二乘回归估测
图80.50 mm 土壤有效磷含量的偏最小二乘回归估测
图90.25 mm 土壤有效磷含量的偏最小二乘回归估测
通过对比 2.00、0.50、0.25 mm 3 种粒径的土壤有效磷含量的估测模型可以看出，0.25 mm 粒径的土壤有效磷含量建立的估测模型，模型的趋势线更靠近 y=x 这条实线，且建模集与验证集的模型 R2 分别为 0.83、0.78，均大于 2.00、0.50 mm 所建立模型的 R2；模型的 RMSE 也均低于 2.00、0.50 mm 粒径所建模型的 RMSE。
3 讨论
本研究在对土壤有效磷原始光谱数据进行 10 种光谱变换，利用皮尔逊相关性分析，找出与土壤有效磷含量相关性最高的一种光谱变换，经过 CARS 算法筛选出有效磷的特征波段后利用 PLSR、 BP 神经网络与随机森林 3 种方法，分别构建了有效磷含量的高光谱估测模型，对比不同粒径的模型估测精度后发现，过 0.25 mm 筛处理的模型精度最高，其原因或许是土壤颗粒越小，彼此的结合越紧密，土壤表面也就越平滑，反射率就越大[31]。因此，基于以上研究可推断出土壤过筛越细，土壤颗粒越小，利用土壤光谱特征对土壤有效磷含量的估测能力越强。
3.1 土壤光谱反射率受不同因素的影响
分析土壤光谱反射率及其与土壤属性（如有效磷含量）的关系时，一个关键的考虑因素是土壤的多样性及其地理特异性。利用乌鲁木齐市的土壤样本，通过高光谱技术在室内进行的一系列预处理 （包括风干、研磨、过筛）和随后的光谱采集与分析，揭示了该地区土壤有效磷含量与光谱反射率之间的关系。然而，这一基于乌鲁木齐市土壤样本建立的模型在其他地区的适用性仍需进一步验证，因为不同区域的土壤特性（如母质、有机质含量、质地和氧化物组成）可能存在显著差异，这些都对光谱反射率有着重要影响。此外，本研究由于样本数量的限制，仅考察了 2.00、0.50、0.25 mm 3 种粒径土壤有效磷的光谱反射率，未来研究可以扩展到更多不同粒径的土壤样本，以评估粒径变化对模型精度的影响。通过比较不同粒径土壤的光谱反射率及所建立模型的预测能力，可以更全面地理解土壤粒径对土壤属性遥感监测的影响，从而为提高模型的普遍适用性和预测精度提供科学依据。
3.2 相关性受光谱数学变换的影响
土壤原始光谱特征与有效磷含量的相关性波长在 400～646 nm 内表现为正相关，波长在 647～1191 nm 内的相关性为负相关，随后又变为正相关。然而原始有效磷光谱反射率与有效磷养分含量相关性不高，大量的研究证明，通过对原始光谱数值进行一阶微分、二阶微分、连续统去除等变换后，能够显著提高光谱特征与土壤特性参数的相关性[32-34]。因此，文中通过 10 种光谱变换后的数据与有效磷养分含量进行相关性分析得出，基于对数的一阶微分变换后的光谱数据与有效磷含量的相关性最高，最大相关系数为 0.492。然而，对土壤有效磷原始光谱数据进行变换的方法并不止这 10 种，未来可以探究更多变换形式下土壤有效磷含量相关性更高的光谱变换形式，建立估测精度更好的模型。
4 结论
对土壤原始光谱反射率进行不同的光谱变换，能够显著提升光谱反射率与土壤有效磷含量的相关性，其中以对数一阶微分变换的效果最好，采用对数的一阶微分建立有效磷含量的 PLSR 模型，在不同粒径处理中都能够较为准确地估测土壤有效磷的含量，模型的 RPD 均在 2.00 以上；而采用对数的一阶微分建立的 BP 神经网络与随机森林模型，在不同粒径处理中只能粗略的估测有效磷含量，模型的 RPD 均小于 2.00。
对比不同粒径的土壤有效磷含量建立的估测模型发现，土壤过筛粒径越细，所建立的模型精度越高；过 0.25 mm 筛的处理，采用对数的一阶微分变换构建的 PLSR 模型，模型 RPD 为 2.13；BP 神经网络与随机森林回归建立的模型，模型 RPD 分别为 1.86 和 1.94；3 种模型中只有 PLSR 建立的模型能够对土壤有效磷含量进行精确估测，而 BP 神经网路与随机森林只能对土壤有效磷含量进行粗略估测，PLSR 建立的模型最优。