长河流域位于山西省晋城市泽州县境内(35° 30'— 35° 31'N, 112° 36'— 112° 39'E), 包括川底、 大东沟、 下村三个乡镇47个村庄, 属于温带大陆性气候, 年均温11.2 ℃, 最高温40.2 ℃, 最低温-18.7 ℃, 年均降水量为616.8 mm。 该区处于黄土高原东南端, 境内以山地、 丘陵为主, 地形起伏多变, 地势西高东低, 海拔730~1 185 m。 研究区处于黄土向棕壤过渡地带, 土壤成土母质主要为黄土性和黄土状母质, 土壤类型主要为褐土, 是黄土高原南部主要分布土类, 易粘化, 板结效应明显^[9]。 区域内主要粮食作物是玉米、 小麦、 农地利用方式主要为玉米-小麦旱作方式。 长河贯穿南北全境, 属于沁河支流。 流域内有12个正在开采的煤矿几乎覆盖流域西部, 煤矿的生产规模均在45万t· a^-1以上。 受长期采煤影响, 流域内地表沉陷严重, 土壤碳库受到极大扰动。

研究区土壤类型为褐土, 质地类别是壤土, 土壤粘粒含量为12.90%~47.70%; 土壤pH范围6.49~8.51, 均值为7.97; 土壤有机质含量为7.80~20.60 g· kg^-1, 均值为13.18 g· kg^-1。 于2015年7月, 以当地矿区为中心, 在考虑地形图和实地勘探基础上, 紧邻矿区的四周采用“ Z” 形路线采样, 使用网格法进行布点, 网格规格为1 km× 1 km, 在网格中选取相对平整的土壤进行样本采集, 在5 m半径内, 随机选择3个不同位置采集表层土壤组成1个混合样, 采集深度为0~20 cm, 若遇到房屋、 道路、 河流等不适合采样的情况则选择临近同类土壤采取, 共采集117个表层土壤样品, 每份土样约1 kg。 使用Garmin GPS72记录采样点坐标, 以确定其于遥感影像中的对应位置。 为了减少混合像元的影响, 提高端元丰度, 在采样后又重新在研究区进行调研, 根据地理坐标对样本点进行筛选, 在保证样本点覆盖所有地貌类型的前提下, 剔除影像单个像元范围内(30 m× 30 m)出现建设用地的样本点, 避免建设用地平滑表面导致的过高反射率对同像元内耕地反射率造成噪声干扰, 使得提取的褐土带耕地反射率真实反映耕地的光谱特性。 筛选后采样点共63个, 如图1所示。 带回室内, 自然风干, 研磨后剔除植物根系、 小石块、 动物残体等杂质, 过2 mm筛子, 然后采用重铬酸钾(K₂Cr₂O₇)— 外加热法测定样本SOC含量^[10]。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区土壤采样点分布Fig.1 Distribution of soil sampling points in study area

为满足经典统计分析建模的要求, 对SOC数据进行正态性检验, 结果服从正态分布。

使用的Landsat8 OLI影像于https: //earthexplorer.usgs.gov/下载。 Landsat8卫星OLI传感器相较于Landsat7的ETM+传感器增加了Coastal波段, 波长更短, 利于气溶胶性质的检测。 考虑到7月份遥感影像中植被、 云量对于SOM反射光谱造成的信息干扰, 选用2013年— 2015年春季、 秋季的影像, 轨道号125-035, 云量低于2%, 植被较少, 地表裸露, 可以较好的还原SOM原有的光谱信息, 有利于提高SOM的反演精度。 使用ENVI5.3对影像进行包括裁剪、 几何校正(校正误差小于1个像元)的预处理。

整体技术路线流程见图2。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 整体技术路线图Fig.2 Flow chart of overall technology

1.4.1 大气校正

使用ENVI5.3对预处理后的影像进行大气校正, 同时利用RSD(remote sensing desktop)针对Landsat8 OLI的6SV大气校正模块, 调用成像当天辅助气象数据, 对影像进行大气校正。

ENVI中的FLAASH大气校正模型是基于MORTRAN4+内核对定标后的影像进行大气校正。 它假设太阳波谱范围内(不包括热辐射)地表为标准的平面朗伯体^[11], 模型中辅助气象数据的输入局限于几种固定的模式, 对于地形复杂区域大气校正的精度受限。

RSD引进了USGS基于6SV2.1模型为Landsat8设计的大气校正程序Landsat 8 Surface Reflectance Code(LaSRC)。 首先引进了LaSRC的核心校正代码并使之适应Windows窗口环境, 改造和重写了属性文件解析、 辅助文件导入以及全部I/O部分代码, 使之在Windows环境顺利运行。 6SV模型是6S模型的矢量版本, 在精确地模拟卫星和平面观测高度、 合理考虑异向性和非均一地表、 非朗伯体等因素的基础上, 可以对连续有序的散射进行模拟, 计算辐射极化, 进而求解辐射传输方程, 能够更加有效地消除瑞利散射和气溶胶的影响。

辅助气象数据由USGS通过MODIS逐日为LaSRC计算而得, 与FLAASH相比, 大气压、 水蒸气、 臭氧等参数有更高的时间、 空间分辨率。 此外Levy和Vermote等发现红、 蓝波段比中红外、 蓝波段反演气溶胶精度更高, 以此结合Landsat8 OLI Coastal波段可以进行更高精度的大气校正。 6SV模型传感器接收到大气反射率为

R* (θS, θV, ΦV)=Tg(θS, θV){Ra(θS, θV, ΦV)+T(θS)1-RPS[RCe-τμV+RPTd(θV)]}(1)

式中, R^* (θ _S, θ _V, Φ _V)为传感器所接受的大气顶层表观反射率, R_a(θ _S, θ _V, Φ _V)是瑞利散射和气溶胶散射引起的路径辐射, θ _S为太阳天顶角, θ _V为观测天顶角, Φ _V为方位角; S是大气球面反射率; T(θ _S)为下行辐射总透射率, T_d(θ _V)为上行辐射总透射率, e-τμV是直接投射到传感器的上行辐射, μ _V=cosθ _V, 为卫星天顶角余弦; τ 是大气光学厚度; R_P为非均一目标反射率, R_C表示邻近效应^[12]。

1.4.2 光谱特征分析及特征波段选择

依据泽州县长河流域土地利用现状图和实地调查资料, 对FLAASH, 6SV影像进行SOC的光谱曲线分析, 掌握矿区褐土SOC的光谱特征。 利用波段数据进行PLSR, 根据回归系数进行波段选择。 最终在光谱特征分析和回归系数等信息的基础上, 选择对SOC反映较好的特征光谱波段。

1.4.3 模型构建与精度比较

选取50个样点构成建模样本集, 13个样本点构成验证样本集, 在所选特征波段基础上, 利用两种大气校正影像各波段反射率及其变换形式, 分别建立MLR, BP神经网络和PLSR的SOC反演模型。 经过建模比较, FLAASH大气校正影像MLR, BP神经网络模型选用平方根( R), PLSR模型选用log(1/R), 6SV大气校正影像选用一阶微分(R')的变换形式。 MLR在SPSS 22中实现, BP神经网络在MATLAB R2015b中实现, PLSR在SIMCA-P 13.0中实现。 用验证集对模型精度进行检验, MLR模型精度评价指标为决定系数R², 统计量F及显著性水平Sig。 BP模型精度评价指标为决定系数R²、 均方根误差RMSE。 以决定系数R²、 累积交叉有效性系数Q²评价PLSR模型的精度。 大量模拟实验后确定BP模型参数设置如表1所示。

表1

Table 1

表1(Table 1)

表1 BP模型主要参数设置 Table 1 Main parameters of BP model BP数据源 BP神经网络 隐含层节点数 学习速率 FLAASH影像 3 0.01 FLAASH'影像 3 0.01 6SV影像 4 0.01 6SV'影像 4 0.01 Note: ' indicates mathematical transformation of the original data 注: '表示对原始反射率进行数学变换

表1 BP模型主要参数设置 Table 1 Main parameters of BP model

为突出SOC变动特征, 以土壤有机质含量等级变化揭示有机碳含量变化特征。 由表2可知, 研究区土壤有机质含量为7.80~20.60 g· kg^-1, 均值为13.18 g· kg^-1。 根据第二次土壤普查养分分级标准, 以三级(20~30 g· kg^-1)、 四级(10~20 g· kg^-1)、 五级(6~10 g· kg^-1)3级别划分土壤有机质含量。 含量水平主要集中于四、 五级两级, 总体水平较低。 建模集和验证集的褐土样本统计特征与总集接近, 说明两个子集可以很好代表总体样本。 根据野外采样点实测数据, 绘制并分析2015年三种影像上各级土壤有机质含量的反射率, 结果如图3所示。

表2

Table 2

表2(Table 2)

表2 土壤有机质含量统计特征 Table 2 Statistical characteristics of SOM content 样本集 样本数 最小值/ (g· kg-1) 最大值/ (g· kg-1) 均值/ (g· kg-1) 总集 63 7.80 20.60 13.18 建模集 50 7.80 20.60 13.26 验证集 13 9.00 16.80 12.86

表2 土壤有机质含量统计特征 Table 2 Statistical characteristics of SOM content

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 褐土光谱曲线Fig.3 Spectral curves of cinnamon soil

总体而言: (1)图3(a)和(b)两幅影像与图3(c)影像相比, 可见光波段反射率均有一定程度的降低, 而红外波段反射率均得到一定程度的提高, 对光学窗口内的波段而言, 辐射亮度失真的关键在于散射, 受大气氧分子、 水分子吸收作用的影响并不明显, 但在更长的红外窗口波段, 大气的主要影响是吸收, 而非散射。 瑞利散射和无选择性散射直接作用于可见光波段, 一部分未到达地面的辐射经大气散射后直接进入传感器, 传感器接收到无效信号, 导致图像质量下降。 而水汽、 CO₂和O₃等对于红外波段的强吸收, 使得传感器接收到该波段的反射能量衰减, 真实反射光谱信息丢失严重。 通过大气校正摒除了大气散射、 吸收的影响后, 无效辐射能量对光谱细节干扰明显减弱。 在可见光波段, 土壤反射率降低; 在近红外和短波红外波段, 土壤反射率增高。 说明FLAASH大气校正和6SV大气校正均可减弱大气吸收和散射的影响; (2)即使不同采样点土壤有机质含量不同, 但是大气校正以后它们的土壤光谱曲线形态十分相似, 且总体呈递增趋势。

图3(a)为FLAASH影像各级土壤有机质含量的光谱反射率, 其反射率在可见光波段较低, 五级有机质的反射率相较于三、 四级有机质低, 而四级有机质反射率相较于三级有机质含量来说要高, 呈现出土壤有机质含量与反射率不相对应的状态, 各级土壤有机质可见光波段反射率可分性差, 这是因为FLLASH大气校正假设在太阳波谱范围内(不包括热辐射)地表为标准的平面朗伯体, 并没有考虑地表非均一的状态, 对于地形起伏较大、 地物类型多样的环境校正效果并不准确, 并且FLAASH大气校正对于气溶胶、 大气模型等参数相较6S大气校正分为几种固定模式, 并没有影像当天详细的气象数据, 在微观区域上进行大气校正精度不足。 四级土壤有机质红外波段与可见光波段差值仅为0.08, 而波动最大的三级土壤有机质红外波段与可见光波段差值不超过0.19, 说明FLAASH影像包含的信息量较少, 区分不同等级土壤有机质含量的能力较低。

图3(b)为6SV影像各级土壤有机质含量的光谱反射率, 其反射率在可见光波段较低, 三级土壤有机含量可见光波段反射率最低, 五级土壤有机质含量可见光波段反射率最高, 土壤有机质含量与反射率各自对应, 反射率可分性强。 这是因为6SV模型不仅考虑了地表为非均一朗伯平面体的情况, 并且结合从MODIS影像中计算得到的当天的气象辅助数据, 对于大气环境、 气溶胶等参数有精确的记录, 可以计算极化辐射, 能大幅度提高地表反射率的提取精度。 土壤有机质含量与土壤可见光(1~4 Band: 0.43~0.45, 0.45~0.51, 0.53~0.59, 0.59~0.67 μ m)反射率总体上呈负相关关系, 反射率随土壤有机质的增加而减少, 沙晋明等的研究也反映了这个规律。 一方面因为有机质含量低, 湿度低、 质地粗的土壤对可见光的吸收率低, 反射率高, 而研究区土壤为褐土, 湿度低, 质地粗, 故而会出现低含量土壤有机质反射率高于高含量土壤有机质的现象; 另一方面由于采矿活动中穿孔、 凿岩、 爆破、 装运、 破碎等工序的影响, 耕地表面矿物粉尘覆盖增多, 且褐土本身铁元素会发生水解, 氧化和迁移比较明显, 氧化铁含量较未扰动土壤增加, 会导致可见光波段土壤反射率的明显下降, 这也在一定程度上说明了6SV影像可以精准反映褐土光谱特性。 图3(b)中三级土壤有机质含量因为大于2%, 所引起的土壤反射率的下降可以掩盖其他成分的光谱特征。 导致近红外波段受到铁氧化合物的吸收作用影响减弱, 反射率得到提高。 由图3(b)可知各级土壤有机质可见光波段与红外波段差值均达到0.18, 且波动最大的三级土壤有机质红外波段与可见光波段差值达到0.25, 说明6SV影像包含的信息量较多, 6SV影像区分不同等级土壤有机质能力较高。

以PLSR分别对FLAASH影像、 6SV影像7个波段的反射率与土壤有机质含量之间的关系进行分析, 并绘制回归系数曲线, 如图4所示。 FLAASH影像的曲线于1— 3波段平缓上升, 3— 5波段上升速度加快, 5— 7波段迅速下降, 最高峰位于5波段; 6SV影像曲线于1— 2波段短暂上升后, 在2— 4波段连续下降, 4— 5波段又转为上升, 5— 7波段间下降上升波动强烈, 最高峰同样在5波段。 两幅影像曲线在4— 7波段均存在波峰或波谷, 与图3光谱曲线特征相近, 可见4— 7波段对于反映土壤有机质含量作用明显。 此外将FLAASH影像与6SV影像回归系数取绝对值对比, 后者明显大于前者, 即6SV大气校正所得地表反射率对于褐土带矿区土壤有机质信息细节捕捉更优, 除5波段外, 两影像波段回归系数绝对值差值较小, 说明波段间存在相关性。 因此以6SV影像为主要参考, 兼顾FLAASH影像, 综合考虑曲线变化特征波段、 回归系数与邻近波段的相关性, 最终确定土壤有机质的特征波段4— 7。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 偏最小二乘回归系数Fig.4 Partial least squares regression coefficients

对于FLAASH影像和6SV影像, 以选择的特征波段反射率及相应数学变换形式为自变量, 以土壤有机质含量为因变量, 分别通过MLR, BP神经网络和PLSR建立反演模型, 结果见表3。

表3

Table 3

表3(Table 3)

表3 SOC反演模型与精度 Table 3 SOC retrieval models and accuracies 模型名称 评价 指标 FLAASH 影像 FLAASH' 影像 6SV 影像 6SV' 影像 多元线性 回归 R2 0.063 0.067 0.765 0.765 F 1.818 1.879 40.879* * 40.820* * BP神经 网络 R2 0.385 0.389 0.764 0.767 RMSE 0.133 0.116 0.190 0.189 偏最小二乘 回归 R2 0.060 0.141 0.775 0.778 Q2 0.003 0.093 0.706 0.764 Note: * * means significant at p< 0.01 level 注: * * 表示p< 0.01水平显著

表3 SOC反演模型与精度 Table 3 SOC retrieval models and accuracies

在未进行数学变换的情况下, FLAASH的MLR模型的R²仅为0.063, 稍优于PLSR模型, 最佳模型为BP模型, R²为0.385, RMSE为0.133。 6SV影像的MLR模型R²为0.765, 优于BP模型, PLSR模型R²高于MLR模型, 为最优模型。

在对土壤反射率数据进行相应的数学变换之后, FLLASH影像所有模型的精度均得到提高; 6SV影像MLR模型F检验值低于变换之前, 模型精度略有下降, 除此之外其他模型精度均有相应改善。 这在一定程度上说明对土壤反射率进行适当的数学变换可以降低影像预处理过程中没有消除掉的噪声。 FLAASH影像的PLSR模型R²为0.141, 优于MLR模型, 最佳模型仍为BP模型, R²为0.389, RMSE为0.116。 6SV影像的MLR模型的R²为0.765, 为三种模型中最低。 BP模型略优于MLR模型, 最优模型仍为PLSR模型, R²高达0.778, Q²为0.764, 接近1, 说明模型预测效果与逼近效果相近, 建模效果良好。

由结果可知, 6SV影像的建模效果远远优于FLAASH影像。 这表明6SV大气校正因为考虑到了非均一朗伯体的问题而较FLAASH在微观区域的大气校正上精度有所提高。 此外由于气溶胶、 大气环境等辅助性数据的输入并非仅仅局限于几种固定的模式进行套用, 而是通过USGS官网所提供的相关辅助数据中进行特定区域特定时间的气象信息读取, 在空间和时间分辨率上较普通的FLAASH大气校正有了明显的提高, 所以大气校正后地表反射率的精度更为接近真实的地表反射率, 在反演相应的土壤理化性质上更有优势。

再者, 由上述比较可知对土壤反射率进行适当的数学变换一定程度上可以提高模型的精度, 这是因为数学变换可以降低影像预处理过程中没有消除掉的噪声, 使得土壤有机质与地表反射率之间的关系更为显著。

从具体模型上来看, MLR在三种模型中建模效果最差, 这可能是由于不能完全消除特征波段中存在的线性相关问题, 导致函数逼近过程中一部分重合信息被多次利用, 降低了模型的精度。 BP模型优于MLR模型而低于PLSR模型, 主要因为BP神经网络对于样本量的要求低于经典的统计分析模型, 所以在样本量较少的情况下BP模型精度要高于MLR模型, 但是BP模型样本依赖性严重且一定程度上会有过拟合现象的出现, 所以在预测精度上要低于PLSR模型。 而由于PLSR模型对数据进行降维后可以更有效地把握数据的特征信息, 更容易辨识系统信息与噪声, 甚至一些非随机噪声, 避免了因子内部高度线性相关、 因子结构不确定性和模型不能识别等问题, 在提取出自变量最大变异信息的同时, 还具有一个响应矩阵对因变量进行最优模拟, 有效克服了典型相关分析、 多元线性相关分析和主成分分析的不足, 并继承其优点^{[13, 14]}, 因此它的预测效果要优于一般的模型。

由图5基于6SV影像的三种模型验证散点图可知, PLSR模型检验样本拟合效果最好, R²为0.702, BP神经网络次之, R²为0.696, MLR模型验证效果最差, R²为0.675。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 验证散点图Fig.5 Scatter plot of validation

利用6SV影像和样本数据建立PLSR模型, 对研究区2013年— 2015年土壤有机质含量进行预测, 并乘以Van Bemmelen因数, 得到有机碳反演结果, 如图6、 表4所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 2013年— 2015年SOC反演结果Fig.6 SOC retrieval results from 2013 to 2015

表4

Table 4

表4(Table 4)

表4 土壤有机碳反演结果统计表 Table 4 Statistical table of SOC results 土壤有机碳含量 /(g· kg-1) 占图幅百分比/% 2013 2014 2015 < 5 0.14 0.04 0.09 5~10 99.82 95.75 84.04 10~15 0.04 4.21 15.87

表4 土壤有机碳反演结果统计表 Table 4 Statistical table of SOC results

褐土矿区有机质含量中部高、 两侧低, 有机碳含量最低区域集中在西北部。 西北部地区海拔高, 丘陵密布, 重力作用加速地表物质的迁移, 土层薄弱, 抗蚀抗冲性差, 为矿区最早开发的地区, 扰动时间长, 扰动强度大, 地表沉陷、 水土流失严重, 地表裂缝分布多, SOC含量降低; 中部地区开发晚, 扰动强度较低, 故而含量较高。 由表可知, 自2013年矿区进行复垦以来, SOC含量< 5 g· kg^-1的土地面积逐年减少, 而SOC含量5~15 g· kg^-1的土地面积逐年增多, 截至2015年, SOC含量< 5 g· kg^-1的土地由0.14%降低至0.09%, 而SOC含量10~15 g· kg^-1的土地由0.04%增加至15.87%。 结果表明, 长河流域褐土矿区复垦稳固甚至增加了SOC的含量, 复垦对于矿区生态稳定和环境恢复具有积极作用。