2020, Vol. 29
PDF
2. 吉林建筑大学, 吉林 长春 130118;
3. 二十一世纪空间技术应用股份有限公司, 北京 100096;
4. 大庆钻探地质录井二公司, 吉林 松原 138000
2. Jilin Jianzhu University, Changchun 130118, China;
3. Twenty First Century Aerospace Technology Co., Ltd., Beijing 100096, China;
4. The Second Logging Branch, Daqing Drilling Engineering Corporation, Songyuan 138000, Jilin Province, China
土壤中有机碳（Corg）包括土壤中植物、动物、微生物及其分泌物，是衡量土壤肥力的重要指标，是全球碳循环的重要指示剂[1-2].土壤中有机碳在可见光、近红外都具有独特的光谱特性，可以根据遥感光谱反射率及其数学变化来模拟与有机碳含量之间的关系.传统的采用地球化学方法估算土壤有机碳库的方法效率低，大面积使用不现实，因此利用遥感手段开展大面积Corg含量反演研究，可以快速、高效地获取Corg含量值并探索其空间分布规律.当前对Corg的空间分布预测主要利用地统计法和遥感反演法[3-5]，还有基于点面扩展的土壤学方法[6].然而由于土壤的空间变异显著，使用地统计法需要大量样本点的数据支撑[3-4]，否则无法保证样点的代表性.部分学者开展综合数据分析[7-11]，另一些学者尝试将遥感数据作为辅助变量与地统计法结合[12-15]，遥感数据在研究区内为地统计法空间插值Corg提供了更多的空间变异信息.同时，在建模方法上进行了多种尝试[14-15].刘焕军、Hill、程朋根等证明了高光谱、多光谱遥感数据均可用于建立Corg定量预测模型[16-21]，这为利用遥感反演Corg提供了理论依据.
黑龙江海伦地区是东北典型黑土区，耕作方式、气候等因素决定了该地区裸土时间较长，适合于开展大面积土壤中有机碳的遥感反演研究.本研究以海伦地区裸土期耕地为研究对象，利用Landsat8 LOI卫星获取便利、光谱范围宽、覆盖面积广等诸多特点，结合该地区土壤地球化学采样分析结果，建立多波段光谱反射率及其构建的光谱指数与土壤有机碳之间的多元逐步回归预测模型，探讨典型黑土区有机碳遥感快速、高效反演方法，为估算东北黑土区有机碳库提供支持.
1 研究区概况及采样点分布海伦市是绥化市的3个县级市之一，为黑土地保护利用试点区，位于黑龙江省中部绥化市北部，南接绥化市北林区，北邻北安市，南距哈尔滨市214 km.区内土壤种类主要为典型黑土，是黑土亚类中最接近黑土概念的亚类.黑土的母质为黄褐色、深黄褐色亚黏土，黑土表层为松散的灰黑色和黑色腐殖质层，向下渐变过渡到淋溶层，其下为淀积层，淀积层发育良好.
基于2015年该地区地球化学采样分析结果，在区内旱田抽取了164个样点作为黑土有机碳反演的建模数据，采样点分布如图 1所示.其中工作区东北部主要分布林地，没有设计采样点.
|
图 1 研究区采样点分布图 Fig.1 Distribution of sampling sites in the study area 1-采样点（sampling site）；2-地名（place name）；3-公路（highway）；4-行政区界线（administrative division） |
本次研究主要使用的是Landsat8 OLI的Level 1T级别的数据，是经过辐射校正和几何校正的数据产品.出于研究的需要，为了使遥感数据能够具有更高的精度，对影像还需进行必要的处理. OLI数据的预处理主要包括辐射定标和大气校正.
2.1 辐射定标所谓辐射定标，简单来说就是将传感器图像的数字量化值（DN）转化为辐射亮度值或反射率或表面温度等物理量的处理过程.辐射定标是大气校正的准备工作.
在ENVI5.1软件中，Landsat8 OLI遥感影像数据的辐射定标具体步骤是：读取该景影像的头文件，在选择定标数据类型后，选择辐射亮度（Radiance）为该影像的的定标类型；选择BIL为定标的输出存储顺序以便进行下一步的大气校正过程；在将输出的数据类型设置为Float后，可以自动获取辐射亮度的转换系数为0.1；设置好所有参数后，输出定标文件.
2.2 大气校正遥感传感器所接收的辐射能量是经大气层发生散射、吸收的衰减能量，使其反映到遥感影像的光谱发生了变化.同时，不同的太阳、地物、传感器之间的几何位置关系，也使影像受到不同程度的大气影响.利用FLASSH模块对影像数据进行大气校正，消除大气散射的影响.随后将相邻两景影像中的一景影像作参考，将另一景影像中的相同点的辐射亮度值调至相同或相近，进行相对辐射校正.本次试验优先选择春天4~5月份数据，利用小波变换方法对2景影像进行光谱差异校正.
3 黑土有机碳反演模型建立 3.1 逐步多元线性回归模型多元线性回归按规则对输入的自变量进行筛选，找出最优组合建立回归方程.它要求自变量与因变量之间为线性关系，并且自变量之间不存在严重的共线性，因变量符合正态分布.其筛选自变量的方法有强行进入法、消去法、向前选择法、向后剔除法和逐步回归法5种.逐步回归法是向前选择和向后剔除法的结合，根据全部自变量对因变量的贡献值大小来选择优先进入方程的自变量，即通过偏F检验来筛选变量，F统计量包括了F进和F出，其计算方法为：
当模型已经引入了x1，…，xi，…，xk这k个自变量，在引入新的第k+1个自变量时，
| $ {F_{进}} = \frac{{{\mathop{\rm SSE}\nolimits} \left( {{x_1}, \cdots , {x_k}, {x_{k + 1}}} \right) - {\mathop{\rm SSE}\nolimits} \left( {{x_1}, {x_2}, \cdots , {x_k}} \right)}}{{{\mathop{\rm SSE}\nolimits} \left( {{x_1}, {x_2}, \cdots , {x_k}} \right)/(n - k - 1)}} $ | (1) |
在剔除第i个自变量xi时，
| $ {F_{出}} = \frac{{{\mathop{\rm SSE}\nolimits} \left( {{x_1}, \cdots , {x_k}} \right) - {\mathop{\rm SSE}\nolimits} \left( {{x_1}, \cdots , {x_{i - 1}}, {x_{i + 1}}, \cdots , {x_k}} \right)}}{{{\mathop{\rm SSE}\nolimits} \left( {{x_1}, {x_2}, \cdots , {x_k}} \right)/(n - k - 1)}} $ | (2) |
式中SSE为因变量的预测误差平方和，n为样本个数.
在选择自变量进入模型时，按0.05的水平进行显著性检验.从模型中剔除自变量时，按0.1的检验概率值进行.为了避免变量的进出循环，一般取偏F检验拒绝域的临界值为F进＞F出.在自变量引入后，对这些自变量进行偏F检验.若不通过，则把它从回归模型中剔除；若都通过；则进行下一步新自变量的引进.有可能在新自变量引入时，原被剔除的自变量解释作用又变得显著，通过检验，则还把它重新选入.直到所有模型外部自变量都不能通过偏F值检验时，回归过程结束.因此，与其他建模方法相比，逐步回归法不仅能通过变量选择来减小共线性影响，防止过度拟合，还能考虑到自变量之间存在的互相传递的相关关系，从而提高模型的精度.
3.2 样本奇异点去除及数据集划分样本的收集有一定的随机性，同时由于样本的测量受到各种因素的影响，会存在一些误差甚至错误.在实际检测中，检测的样本数据中总是会出现与总体数据偏离，或者不符合正常规律的数据，称之为奇异点.在建模过程中即使混入很少的奇异点，也可能对模型预测和分析产生严重干扰.因此在建立模型前，需要对样本进行筛选，去除奇异点，才能使建立的模型获得较高的精度.对于光谱定量反演，样本数据异常包括光谱异常和待测属性值异常两部分.
常见的剔除样本奇异点的方法包括：主成分的得分图、拉依达法则、马氏距离、杠杆值、预测残差等，依据的原理是判断某统计量是否超过正态分布或χ2分布等形式的临界值.首先对黑土样本的属性测量值（有机碳含量）剔除奇异点，采用3倍标准差准则来确定.计算所有的样本点有机碳含量（y1，…，yn）的平均值y以及标准差σ，当偏差|vi|=|yi-x|＞3σ时，认为是异常点，予以剔除.然后采用主成分分析对光谱数据进行奇异点筛选.将光谱数据进行主成分分析后，前两个主成分对所有光谱数据的累积贡献率为96.8%，将位于置信区间外的样本点剔除.
将筛选完毕的所有样本点按有机碳含量从低到高排列，进行训练集和验证集的划分.划分主要遵循两方面的规则：首先两个数据集需要有最大程度的相似性，在组成性质上基本相同，保证建立的模型对验证集有适用性；训练样本的物理和化学性质范围需要涵盖验证样本的范围，且分布比较均匀.
将164个黑土样本按有机碳含量从低到高排列，从第一个样本开始，每4个样本抽取一个放入验证集，剩下的样本为训练集，建模样本和预测样本的比例为3 : 1.从训练集和验证集的有机碳含量描述统计（表 1）可以看出，训练集的有机碳含量范围包括了验证集的范围，而且验证集的主成分分值位于训练集的分值范围内，说明验证集有一定代表性.
|
表 1 黑土样本集划分后的描述性统计 Table 1 Descriptive statistics of black soil samples by set division |
分别提取120个建模样本和44个验证样本在Landsat8 OLI数据中各波段反射率值，用于建模.
3.4 模型建立用多元逐步回归建立黑土有机碳含量的估测模型时，若用全谱段的反射率作为自变量，数据量大而且自变量之间的共线性会非常严重，土壤样本也远远小于自变量个数，从而导致模型无法求解或“过拟合”，所得结果也并不可靠.因此根据土壤有机碳与光谱反射率的相关性分析结果，采用相关性较大的波段来建立模型。需要考虑的是，在挑选波段时，并非只关注相关程度最高的波段，而是综合考虑在不同光谱区间的，相关系数较高且具有一定代表性的波段来进入模型.
分别对黑土样本的光谱反射率（R）、反射率的对数（lgR）、反射率的一阶微分（R′）、反射率的倒数（1/R）建立多元逐步回归线性模型. 表 2为各种变换形式的入选波段及其相关系数.
|
|
表 2 土壤光谱反射率及变换形式与有机碳含量的特征波段及其相关系数 Table 2 Characteristic bands and correlation coefficients of organic carbon content for soil spectral reflectance and variation forms |
根据上述有机碳与各波谱形式相关分析结果，采用反射率倒数形式与土壤有机碳建立有机碳含量估测模型，预测结果精度如表 3所示.
|
|
表 3 多元逐步回归模型评价指标分析结果 Table 3 Analysis results of evaluation indexes for multiple stepwise regression model |
将1/B1，1/B2，1/B3，1/B4，1/B5，1/B6，1/B7作为自变量，而将有机碳作为因变量进行逐步回归分析.经过模型识别，最终余下1/B1，1/B4共2项在模型中，R方值为0.180，意味着1/B1，1/B4可以解释Corg的18.0%变化原因.而且模型通过F检验（F=12.833，p=0.000＜0.05），说明模型有效.模型公式为：Corg=1.644+0.056×1/B1-0.038×1/B4.另外，针对模型的多重共线性进行检验发现，模型中VIF值全部小于5，意味着不存在着共线性问题；并且D-W值在数字2附近，说明模型不存在自相关性，样本数据之间并没有关联关系，模型较好.
最终具体分析可知：1/B1的回归系数值为0.056（t=3.975，p=0.000＜0.01），意味着1/B1会对Corg产生显著的正向影响关系；1/B4的回归系数值为-0.038（t=-2.040，p=0.044＜0.05），意味着1/B4会对Corg产生显著的负向影响关系.总结分析可知：1/B1会对Corg产生显著的正向影响关系.另外，1/B4会对Corg产生显著的负向影响关系.据此做出基于土壤光谱反射率倒数形式与土壤Corg含量的散点图（图 2）.
|
图 2 土壤光谱反射率倒数变换与Corg含量散点图 Fig.2 Scatter diagram of soil spectral reflectance reciprocal transformation vs. Corg content |
以ArcGIS软件为解译平台，建立水体、冲沟及建筑用地的解译标志，进行人机交互解译，提取水体、冲沟以及建筑用地信息（如图 3），目的是去除一般与黑土无关的一些土地利用类型.
|
图 3 水体、冲沟以及建筑用地掩膜信息 Fig.3 Mask information of water bodies, gully and construction land 1-冲沟（gully）；2-城镇用地（urban land）；3-工矿用地（industrial and mining land）；4-水库（reservoir）；5-池塘（pond）；6-河流（river）；7-湖泊（lake） |
结合反演模型以及遥感影像解译结果，对研究区内“纯”黑土进行实地反演，对黑土土壤有机碳含量进行分类.从反演分类结果（图 4）可以看出，研究区Corg含量总体水平较高，最高值大于4×10-3，Corg有由东北未开垦区向西南耕地区递减的趋势，说明西南区域开垦耕地后对土壤中Corg含量影响较大.
|
图 4 研究区土壤有机碳含量分布图 Fig.4 Distribution of soil organic carbon contents in the study area 1-Corg＜1×10-3；2-Corg=1×10-3~1.5×10-3；3-Corg=1.5×10-3~2×10-3；4- Corg=2×10-3~2.5×10-3；5- Corg＞2.5×10-3 |
本研究基于地面土壤地球化学分析结果和Landsat8 LOI遥感影像，建立了海伦地区的土壤有机碳遥感预测模型，得出如下结论：
1）绥化典型黑土区土壤有机碳含量高，空间变异性显著，且耕作方式、气候等因素决定了裸土时间长，因而该区适于区域土壤有机碳含量遥感反演.
2）Landsat8数据质量高、时相全，各波段反射率与有机碳含量均显著相关，其中与1波段倒数变换形式相关性最优，可在黑土区推广该数据的使用.研究区黑土有机碳光谱预测最优模型以倒数为自变量模型，决定系数R2=0.180，RMSE=0.557，预测精度较好、稳定性高，可以用于揭示黑土典型区Corg含量的空间分布特征.同时认为在不对土壤进行地面光谱测试的情况下，直接采用化学分析数据与遥感卫星相关联的方法预测模型拟合度有限.
3）研究区Corg由东北向西南呈递减趋势，Corg遥感反演结果可以为该区土壤退化监测、耕地质量评价、土壤有机碳库估算等提供数据支持，为该区土壤资源保护、土地可持续利用提供决策支持.
4）本次研究追求适用于大面积黑土有机碳反演的方法，所以在数据处理与建模方法方面有待进一步优化.接下来会考虑非线性的模拟方法，选择神经网络、随机森林、支持向量机等方法进行试验.另外要考虑多源数据的融合适用及排除土壤中含水量的影响.
| [1] |
Marchant B P, Villanneau E J, Arrouays D, et al. Quantifying and mapping topsoil inorganic carbon concentrations and stocks:Approaches tested in France[J]. Soil Use and Management, 2015, 31(1): 29-38. DOI:10.1111/sum.12158 |
| [2] |
Six J, Paustian K. Aggregate-associated soil organic matter as an ecosystem property and a measurement tool[J]. Soil Biology and Biochemistry, 2014, 68: A4-A9. DOI:10.1016/j.soilbio.2013.06.014 |
| [3] |
Webster R, Oliver M A. Sample adequately to estimate variograms of soil properties[J]. European Journal of Soil Science, 1992, 43(1): 177-192. |
| [4] |
张新乐, 窦欣, 谢雅慧, 等. 引入时相信息的耕地土壤有机质遥感反演模型[J]. 农业工程学报, 2018, 34(4): 143-150. |
| [5] |
吴才武, 张月丛, 夏建新. 基于地统计与遥感反演相结合的有机质预测制图研究[J]. 土壤学报, 2016, 53(6): 1568-1575. |
| [6] |
陆访仪, 赵永存, 黄标, 等. 海伦市耕层土壤有机质含量空间预测方法研究[J]. 土壤通报, 2012, 43(3): 662-667. |
| [7] |
于芳, 李海明. 基于地统计学和GIS的湖北省土壤有机质空间变异性研究[J]. 江苏农业科学, 2019, 47(15): 282-286. |
| [8] |
韩丹, 程先富, 谢金红, 等. 大别山区江子河流域土壤有机质的空间变异及其影响因素[J]. 土壤学报, 2012, 49(2): 408-408. |
| [9] |
赵明松, 张甘霖, 王德彩, 等. 徐淮黄泛平原土壤有机质空间变异特征及主控因素分析[J]. 土壤学报, 2013, 50(1): 1-11. |
| [10] |
张哲寰, 赵君, 戴慧敏, 等. 黑龙江省讷河市土壤-作物系统Se元素地球化学特征[J]. 地质与资源, 2020, 29(1): 38-43. |
| [11] |
陶培峰, 王建华, 李志忠, 等. 基于高光谱的土壤养分含量反演模型研究[J]. 地质与资源, 2020, 29(1): 68-75. |
| [12] |
Wu C F, Wu J P, Luo Y M, et al. Spatial prediction of soil organic matter content using cokriging with remotely sensed data[J]. Soil Science Society of America Journal, 2009, 73(4): 1202-1208. DOI:10.2136/sssaj2008.0045 |
| [13] |
Mirzaee S, Ghorbani-Dashtaki S, Mohammadi J, et al. Spatial variability of soil organic matter using remote sensing data[J]. CATENA, 2016, 145: 118-127. DOI:10.1016/j.catena.2016.05.023 |
| [14] |
王茵茵, 齐雁冰, 陈洋, 等. 基于多分辨率遥感数据与随机森林算法的土壤有机质预测研究[J]. 土壤学报, 2016, 53(2): 342-354. |
| [15] |
杨佳佳, 白磊, 吴嵩.黑龙江典型黑土区土壤侵蚀遥感监测技术研究[J]. 2019, 28(2): 193-199, 183.
|
| [16] |
丁美青, 肖红光, 陈松岭, 等. 基于BP神经网络的土地开发整理区土壤有机质含量遥感定量反演[J]. 湘潭大学自然科学学报, 2012, 34(2): 103-106. |
| [17] |
程朋根, 吴剑, 李大军, 等. 土壤有机质高光谱遥感和地统计定量预测[J]. 农业工程学报, 2009, 25(3): 142-147. |
| [18] |
刘焕军, 赵春江, 王纪华, 等. 黑土典型区土壤有机质遥感反演[J]. 农业工程学报, 2011, 27(8): 211-215. |
| [19] |
刘焕军, 吴炳方, 赵春江, 等. 光谱分辨率对黑土有机质预测模型的影响[J]. 光谱学与光谱分析, 2012, 32(3): 739-742. |
| [20] |
Hill J, Schütt B. Mapping complex patterns of erosion and stability in dry mediterranean ecosystems[J]. Remote Sensing of Environment, 2000, 74(3): 557-569. DOI:10.1016/S0034-4257(00)00146-2 |
| [21] |
周萍, 王润生, 阎柏琨, 等. 高光谱遥感土壤有机质信息提取研究[J]. 地理科学进展, 2008, 27(5): 27-34. |
| [22] |
李媛媛, 李薇, 刘远, 等. 基于高光谱遥感土壤有机质含量预测研究[J]. 土壤通报, 2014, 45(6): 1313-1318. |