土壤有机质、全氮具有为植物提供营养、促进植物生长发育、改善土壤的物理性质等作用，并指示土壤肥力总体状况^{[1, 2, 3]}。快速分析与预测土壤有机质、全氮的空间分布特征和动态变化趋势不仅为减少耕作成本、防止土壤肥力退化提供理论决策依据^{[1, 2]}，也对实现土壤资源的合理利用及其可持续发展具有重要意义。

采用传统方法对土壤肥力进行区域连续监测时，由于采样和分析成本高、周期长，难以获得实时和连续的土壤信息。遥感以其大区域、实时快速、可重复获取多时相信息、高效、经济等优势，为实现土壤肥力变化动态监测提供了有效途径。近年来，越来越多的研究开始利用植被指数反映表层土壤有机质、全氮含量的空间变化，并取得了一定进展^{[4, 5, 6, 7]}。这些研究基于这样的认识，即植被的生长可通过凋落物、根系等增加土壤有机质输入，促使土壤有机碳储量的增加^[8]。研究表明，植被净初级生产力(NPP)与枯落物质氮、植被的年氮累积量以及土壤氮矿化速率之间均存在较为显著的相关关系，其中NPP与氮矿化速率、年植被氮积累量的相关系数R²分别达到0.74、0.60^[9]。

尽管国内外研究表明土壤有机质、全氮与遥感因子(NDVI、NPP、植被覆盖度等)间存在一定的相关关系，并可以用卫星遥感监测土壤养分含量，但研究仍然处于初步阶段，且大多侧重于土壤肥力因子与遥感因子的直接关系，很少涉及土壤有机质、全氮的动态变化遥感监测及其区域制图研究，因而难以进一步揭示区域土壤肥力变化状况。本研究从土壤有机质、全氮变化的角度出发，通过研究广西中南部土壤有机质、全氮及其变化与多个遥感因子的耦合关系，建立表层土壤有机质、全氮含量变化的遥感监测模型，以期为广西区更准确获取土壤肥力因子变化信息提供一条新的途径。

研究区主要位于广西中南部柳州、来宾、南宁、崇左、钦州等市，地理范围为107°8′E—109°37′E，21°36′N—24°36′N，总面积大约5.4万km²。该区属于亚热带季风气候区，气候温和，雨量充沛，干湿分明。年平均温度17.8—22.4℃，年均降雨量1200—1300 mm，年日照时数1500—1700 h。全年光照充足，且光、水同季，全年夏长冬短，适宜作物四季生长。研究区以丘陵山地为主，整个地势自西北向东南倾斜，山体较多，岭谷相间，中部包含大明山，东南方向为六万大山，西南与四方岭相邻，地形复杂多变，喀斯特地貌突出，高程范围为0—1738 m，坡度为0°—60.5°。研究区林草地分布最为广泛，且分布较为集中，其中林草地面积占79.24%，耕地面积占总面积的14.76%。该区主要土壤类型为铁铝土、初育土、人为土等，三者分别占总面积的50.32%、32.23%和14.18%。

采用点对点比较法对研究区土壤养分变化进行分析研究。1981年的土壤数据来自于研究区各县第二次全国土壤普查资料(部分点于1980—1982年采集)。按照各样点的空间位置描述及其土壤类型进行空间定位，同时获取其土壤属性数据，并于2011年和2012年进行土壤样品采集(以下统称为2011年)。采用五点取样法采集，采样时尽量靠近原普查样点位置采样，取样深度为0—20cm。研究区内共采集了两个时期土地利用类型都为耕地的土壤样品117个。土壤样品测定方法与第二次土壤普查测定方法一致，有机质采用重铬酸钾容量法，全氮采用半微量凯氏法^[10]。

采用AVHRR NDVI以及MODIS NDVI两种时间序列数据。NOAA时序数据是每半月最大化合成的NDVI数据，来源于中国西部环境与生态科学数据中心，时间范围为1981—2006年，空间分辨率为8km×8km。MODIS NDVI数据是每16d最大化合成的NDVI数据，来源于中国科学院计算机网络信息中心国际科学数据服务平台，空间分辨率为250m×250m。运用ENVI等遥感图像处理软件对AVHRR和MODIS NDVI进行影像拼接、重投影、裁剪等数据预处理。

本文应用NDVI遥感数据，基于像元二分模型估算研究区植被覆盖度。像元二分模型的原理是：假设每个像元是由纯植被部分和非植被部分两部分组成，那么像元信息可表达为由纯植被成分所贡献的信息和由非植被成分所贡献的信息之和。任一像元NDVI的值为两部分植被指数的加权平均值和，权重为各部分在像元中的比例，计算公式为：

式中,NDVI为混合像元的植被指数值；NDVI_veg为纯植被像元时的植被指数；NDVI_soil为非植被像元时的植被指数，一般指纯土壤像元的植被指数；f_veg为植被覆盖度。

像元二分模型法求植被覆盖度的公式如下：

式中,NDVI_soil代表纯土壤像元的最小值；NDVI_veg代表纯植被覆盖像元的最大值^[11]。本研究采用研究区不同植被类型的NDVI的最大值代表NDVI_veg(即NDVI_veg=NDVI_max)，NDVI_soil取值为0.023。

本文采用朱文泉等改进的植被净初级生产力(NPP)遥感估算模型^[12]对研究区1981—2011年的NPP进行估算。该模型将植被覆盖分类引入光能利用率模型进行植被第一净初级生产力的估算，即：提取不同植被类型的NDVI最大值，来进行植被光合有效辐射的吸收比例(FPAR)的估算。利用气象数据(温度、降雨、太阳辐射)来实现水分胁迫因子的估算。计算公式如表 1所示。首先进行月植被净初级生产力的估算，然后计算每年总的净初级生产力。其中，气象数据来源于中国气象局国家气象中心，由研究区及其周围19个气象站点以及北海、南宁、桂林3个辐射站点的1981—2011年月值数据集，利用克里格插值(Kriging)法得到。

表 1 基于光能利用率模型的植被净初级生产力计算 Table 1 Calculation of the vegetation net primary production based on the light use efficiency model

编号 Number	公式 Formula	备注 Notes
1	NPP x,t =APAR x,t ×ε x,t	t表示时间，x表示空间位置；APAR(x,t)表示像元吸收的光合有效辐射；ε(x,t)表示像元的实际光能利用率
2	APAR x,t =SOL x,t ×FPAR x,t ×0.5	SOL (x,t)是t月份在像元x处的太阳总辐射量。
3		FPAR(x,t)为植被层对入射光合有效辐射的吸收比例
4		SR为简单比值植被指数，SR max取值为4.14—6.17，SR min取值为1.05
5	ε x,t =Tε1 x,t ×Tε2 x,t ×Wε x,t ×ε	Tε1和 Tε2表示低温和高温对光能利用率的胁迫作用；W ε为水分胁迫影响系数；ε*表示在理想条件下的最大光能利用率。
6	Tε1 x = 0.8+0.02×Topt x －0.0005× Topt x 2	TOP(x)表示某一区域一年内NDVI值达到最高时月份的平均气温
7	Tε2(x,t)= 1.1814/{1+exp[0.2×TOP(x)-10-T(x,t)]}/ {1+exp[0.3×(-TOP(x)-10+T(x,t)]}	当某一月平均温度T( x,t)比最适温度TOP(x)高10℃或低13℃时，该月 Tε2( x,t)值为最适温度时的一半
8	Wε x,t =0.5+0.5×EET x,t /PET x,t	EET( x,t)为区域实际蒸散量(mm)；PET( x,t)为区域潜在蒸散量(mm)

表选项

最大光能利用率ε^*的取值因植被类型的不同而不同。本文最大光能利用率利用Running 等人基于生态生理过程模型对十种植被类型模拟的结果^[13]，其余类型取值为0.389 g C/MJ，其取值如表 2所示。

提取年NDVI最大值，用于研究NDVI序列变化与土壤肥力因子的耦合关系，即在半月或16d合成最大值的基础上，进一步提取年NDVI最大值；提取月NDVI最大值，用于对研究区植被覆盖度和植被净初级生产力的估算。根据采样点坐标，提取采样点所对应的遥感因子的值，并计算其平均值。

为研究遥感因子变化与土壤肥力因子变化之间的相关关系，采用两种方式来表示遥感因子的变化状况，一种是直线斜率法，另一种是差值法。采用一元线性回归分析方法分析数据的变化趋势，反映植被长期变化趋势。回归直线斜率(trend)采用最小二乘法求得。计算公式为：

式中,变量j为年序号，Y为时间序列遥感数据；θ_trend分别表示NDVI、植被覆盖度、NPP的变化趋势线的斜率，计算AVHRR NDVI、MODIS NDVI、植被覆盖度、NPP的θ_trend时，n分别取值为26、12、31和29。θ_trend>0说明其变化趋势是增加的，反之则是减少。

考虑到年遥感数据可能会受气候、环境等因素的影响，为了尽可能准确地分析遥感数据变化与土壤养分变化的耦合关系，同时采用与求土壤养分变化相同的差值法计算遥感因子的变化。为了减少年遥感数据的偶然性，计算其差值时用研究时段中最后两年数据的平均值减去开始两年的平均值。差值的计算公式为：

式中,k、l表示最后两年的年份，m、n表示起始两年的年份，比如计算AVHRR NDVI变化量时，k、l和m、n分别取值为2005、2006和1981、1982；θ_chg表示遥感因子数据的变化量，θ_chg >0说明该时间段内升高，反之则减少。

采用统计分析软件，首先对两期土壤有机质、全氮进行差异性比较分析，并计算各遥感因子的变化量及变化趋势；然后分别进行遥感因子的变化量、变化趋势分别与土壤养分变化量的耦合关系分析；运用逐步线性回归方法建立土壤有机质、全氮变化的遥感监测模型。

将样本划分为建模集和验证集两个部分，对样本按照有机质变化值进行排序，由小到大每隔两个点选择一个作为验证点，共选择39个验证点组成验证集，剩余78个样点组成建模集。

建模集用于构建土壤养分变化监测模型，验证集用于评价所建模型的制图精度。利用ArcGIS软件空间分析模块，以研究区遥感因子为基础，根据建立的土壤养分变化监测模型，估算研究区土地利用类型为耕地的土壤养分变化，进一步提取出验证点上的土壤养分变化值，比较制图结果与实测结果，从而对所建模型的精度进行评价。

为减少异常样本对模型预测和验证精度的影响，本文剔除残差(Residual)大于3倍样本集标准差，以及杠杆值(leverage)大于其3倍建模样本集杠杆值均值的异常样本^[14]。

采用平均误差(ME)和均方根误差(RMSE)来评价模型建立的精度。其公式分别为：

式中，i为预测值，μi为实测值，n为验证样本个数，i为观测点。平均误差(ME)是对模型总体准确度的度量，其值越接近0，说明模型的总体偏差越小，准确度越高；均方根误差(RMSE)是对模型精度及稳定性的度量，其值越小说明模型的精度与稳定性越高。

由表 3可以看出，两期的土壤有机质、全氮含量存在极显著差异(P<0.01)，研究区有机质、全氮均值从1981年的23.11、1.26g/kg分别上升到2011年的28.55g/kg和1.47g/kg，分别增加了23.51%和16.66%。其中，含量上升的样点占总样点的72.12%和64.77%。土壤有机质、全氮的变异系数相近，都属于中等程度变异，且在不同时期二者的变异系数变化不大。土壤有机质和全氮含量的上升，主要是由于耕地长期施肥以及秸秆还田等农业措施增加了土壤有机物质的投入，使得耕地土壤有机质、氮素的投入量大于其消耗量，使有机质和氮素得到积累^{[15, 16, 17, 18]}。

研究区年最大NDVI、年最大植被覆盖度以及年总NPP的均值变化趋势见图 1。1981—2011年，AVHRR NDVI与MODIS NDVI两种时间序列数据的NDVI均值都略有增加；2000—2006年，研究区基于两种遥感数据源估算的植被覆盖度、NPP值接近；近30年来，研究区3个遥感因子均值变化趋势保持一致，都呈缓慢上升趋势，并出现4—6a的周期性波动。结果表明，研究区域整体植被生长状况有所改善。

图 1 研究区遥感因子变化趋势 Fig. 1 The trend of remote sensing parameters in the study area

图选项

本文采用两种NDVI数据源来表述研究区NDVI变化状况，并用于估算其它两个遥感因子值，基于两种数据源结果的一致性与前人研究结果基本相同^{[19, 20, 21]}，说明基于年最大NDVI、年最大植被覆盖度以及年总NPP进行遥感因子变化分析具有可靠性。

提取采样点对应的遥感因子值，并分析与有机质、全氮含量的耦合关系。为避免遥感数据的偶然性，采用相邻两年遥感因子的平均值进行相关性分析。从表 3看出，遥感因子与土壤有机质、全氮存在一定的相关关系，但相关系数不高，最大相关系数为0.232，说明利用低分辨率的NDVI遥感影像数据可能不能很好地反映出土壤有机质、全氮的含量状况。

分别计算NDVI、植被覆盖度以及NPP的θ_trend和θ_chg，并与土壤有机质、全氮含量变化进行相关性分析。因AVHRR NDVI数据与样品数据在时间跨度上基本一致，因此选择AVHRR时间序列数据来分析NDVI与有机质、全氮动态变化的相关性。

结果表明(表 5)，遥感因子变化与有机质、全氮变化的相关关系明显好于二者数量上的静态关系，其最大相关系数分别为0.31和0.32，在0.01水平上达到极显著；从有机质、全氮变化与遥感因子变化趋势和变化量的相关关系来看，NDVI、植被覆盖度的 变化趋势(θ_{trend_ NDVI}和θ_{trend_ veg})和NPP变化量

土壤有机质及全氮的逐步线性回归分析中，因子入选和淘汰出模型的显著性水平设置为0.1，其线性逐步回归方法构建的土壤有机质、全氮变化遥感监测模型结果见表 6。从表中可以看出，预测有机质变化的线性回归模型中只有θ_{trend_ veg}进入模型，预测全氮变化有θ_{chg_ NDVI}、θ_{trend_ NDVI}及θ_{chg_ NPP}三个变量进入模型。说明可以根据这些遥感因子的变化来预测表层土壤有机质、全氮含量变化在本研究区内的分布。

回归模型中的R²值表明，利用这些因子分别能够解释表层土壤有机质、全氮变化的16.9%和20.3%。由于受到多种因素的影响，以往预测模型的决定系数一般在15%—82%范围内^{[22, 23, 24]}，比如在提取遥感因子时，遥感影像的时间、空间分辨率的大小也会影响土壤属性与遥感因子之间的关系；土壤属性在区域上的空间变异的复杂性可能会使所选择的环境因子不能有效地解释其空间变异；此外，区域、尺度、以及土层深度的差异也会对利用环境因子预测土壤属性精度产生一定影响^[25]。

土壤有机质、全氮变化的回归模型(表 6)应用于研究区土地利用类型为耕地的区域，预测土壤有机质、全氮含量变化在空间上的分布及制图(图 2)，并提取验证点的土壤有机质、全氮变化的属性值，并与实际变化值比较，从而评价制图的精度(表 7)。

图 2 研究区耕地土壤有机质变化、全氮变化分布图 Fig. 2 The map of changes of organic matter and total nitrogen

图选项

对研究区耕地类型的有机质、全氮变化值进行统计，可知耕地类型中土壤有机质平均变化6.65 g/kg，全氮平均变化为0.31 g/kg。

从39个点的验证结果(表 7)可以看出，线性模型预测的有机质变化的平均误差(ME)和均方根误差(RMSE)分别为0.71、3.41 g/kg，与相关研究比较，其RMSE小于Thompson等研究结果(均方根误差4.5—4.8 g/kg)^[26]。本研究全氮变化的平均误差和均方根误差分别为0.09、0.35 g/kg，与潘剑君等利用遥感NDVI数据估算溧水水田土壤全氮含量的均方根误差(0.44—0.01 g/kg)比较接近^[5]。这说明模型的预测结果精度和可信度较高，本研究区域制图效果较好。

由遥感因子与表层土壤有机质、全氮的相关性分析结果可以看出，遥感因子变化与有机质、全氮含量变化呈显著正相关关系，且明显好于遥感因子与有机质、全氮含量数量上的静态相关性。已有研究表明，植被指数、单个波段等遥感因子与土壤有机质、全氮含量存在较好的相关关系^{[6, 27]}，但在本研究中，二者的相关系数不高，主要是由于研究区地形复杂多变、耕地破碎、遥感数据分辨率较低引起的。有研究证明，在这种类型的区域，由于遥感信息混合性较高，单一时相静态遥感数据可能不能完全反应土壤有机质、全氮含量状况，采用长时间序列遥感影像数据变化趋势，可以在一定程度上消除遥感信息中不变部分的信息，增强作物长势变化的信息。因而，利用遥感因子变化能更好的反映有机质、全氮变化^[7]。

本研究通过NDVI长时间序列遥感影像数据与地面样点数据相结合的方法，基于多元线性逐步回归构建表层土壤有机质、全氮含量变化的回归方程，间接监测大区域内表层土壤有机质、全氮含量变化，一方面充分发挥遥感数据大范围、便捷、高效的优点，弥补了传统方法的不足，极大地提高了工作效率。另一方面弥补了以往遥感监测方法受混合像元影响的不足，提高了遥感监测的可靠性。

本研究估算结果精度的影响因素很多。首先，由于研究采用两种遥感数据源，尽管估算出的遥感因子值(植被覆盖度、NPP)基本相同，但是不同分辨率的差异仍然对估算结果会造成一定的影响。其次，由于本研究区域较大，土壤属性及地形、遥感因子的空间变异也比较明显，建模样点有限，在一定程度会造成估算精度的降低。此外，土壤有机质、全氮变化受多种因素的影响，本研究选择的环境因子可能不能完全反映二者的变化状况。因此，后续研究应当从多方面考虑影响土壤有机质、全氮变化的因素，如区域差异、土壤厚度、空间分辨率以及气候因素等，以期提高遥感监测土壤有机质、全氮变化的精度。

本文以广西中南部南宁、柳州、钦州、崇左及来宾等地区的耕地肥力变化为研究对象，以1981年、2011年土壤数据为基础，以AVHRR NDVI和MODIS NDVI遥感影像为数据源，通过研究NDVI、植被覆盖度、NPP与土壤有机质、全氮及其变化的耦合关系，建立土壤有机质、全氮变化遥感监测模型，以实现表层土壤有机质、全氮变化区域监测及制图的目的。

(1)近30年来，研究区土壤有机质、全氮变化较为明显。有机质、全氮变化均呈极显著上升趋势，分别增加了5.43g/kg和0.21g/kg；不同时期二者的空间变异系数变化不明显。研究区内AVHRR NDVI与MODIS NDVI时间序列数据变化趋势基本相同，均呈现缓慢上升趋势；基于两种遥感数据源估算的植被覆盖度、NPP值接近，且均呈现缓慢上升趋势。

(2)遥感参数(NDVI、植被覆盖度以及NPP)变化与土壤有机质、全氮变化具有显著正相关关系。利用逐步线性回归方法，建立耕地土壤有机质、全氮变化的遥感监测模型，分别能解释土壤有机质变化及全氮变化的16.9%和20.3%。

(3)根据所建模型对研究区耕地土壤有机质、全氮变化进行预测制图，研究区内耕地表层土壤有机质、全氮变化均值分别为6.65g/kg和0.31g/kg。根据独立样本预测值与实测值比较结果表明，模型在总体上具有较高的准确度和稳定性，预测结果具有较高的可信度。与传统土壤养分变化监测方法相比，利用所建模型对研究区耕地土壤有机质、全氮变化进行预测制图，能够更加快捷、方便地反映二者变化的空间分布特征。