研究区分别位于我国黑龙江省建三江地区和海伦地区(图1), 两个地区均属于东北黑土区, 但在地形地貌、 土壤类型、 农作物类型等方面均存在差异。 建三江地区处于乌苏里江、 松花江、 黑龙江冲积而成的三江平原腹地, 地势平坦, 平均海拔50 m左右, 相对高差10 m左右, 农作物类型主要为水稻。 土壤类型有暗棕壤、 白浆土、 黑土、 草甸土、 沼泽土、 泥炭土、 水稻土。 研究区面积约2 600 km²。 海伦地处松嫩平原东北端, 小兴安岭西麓, 位于东北黑土区的中心区域, 地形为丘陵、 漫岗, 属中温带大陆性气候, 农作物种植类型主要为水稻、 玉米和大豆。 主要土壤类型为黑土、 暗棕壤和草甸土等。 研究区面积约15 00 km²。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区位置分布图Fig.1 Location map of the study area

航空高光谱数据获取采用CASI/SASI成像光谱仪(加拿大ITRES公司), 选择Y-12飞机为搭载平台。 CASI传感器谱段范围为380~1 050 nm, 空间分辨率1.5 m, 波段数72个, 光谱带宽10 nm; SASI传感器谱段范围为950~2 450 nm, 空间分辨率为3.75 m, 波段数100个, 光谱带宽15 nm。 飞行相对高度3 000 m。 建三江航空数据获取时间为2017年4月, 海伦数据获取时间为2018年5月, 地表均处于裸土期。

数据预处理主要包括辐射校正、 几何校正、 大气纠正和光谱重建。 航空成像光谱测量系统配备了用于辐射校正和几何校正的软件, 其中辐射校正软件为RCX(Radiometric Calibration Xpress) 9.3.5.1, 几何校正软件为Geocor 3.0。 处理步骤包括: ①辐射校正; ②传感器姿态数据处理; ③GPS定位数据处理; ④姿态数据与定位数据时间同步与集成; ⑤精细的几何校正。 在此基础上采用大气辐射传输模型对航空高光谱遥感数据进行光谱重建, 软件采用ENVI5.3版本。 首先在ENVI的FLAASH模块中利用大气辐射传输模型进行校正, 大气水汽含量主要根据820, 940和1 135 nm附近的水气吸收谱带的吸收面积与大气水汽含量的统计关系来反演, 在校正过程中及时对水汽含量反演结果进行评估, 反演生成的水汽图像应呈云雾状, 地物的轮廓基本分辨不清。 然后利用地面黑、 白两种定标布的光谱测量数据, 在ENVI的empirical line模块中进行地-空回归校正, 计算地面光谱与大气校正光谱(基于辐射传输模型)之间的乘性与加性系数, 进一步消除因辐射定标、 波段间相对定标、 波段配准、 大气参数选取等误差因素而造成的光谱误差, 最终获得地物光谱反射率。 反射率数据去除了水汽吸收波段(1 370~1 415和1 835~1 925 nm)及噪声较大波段(380~499 nm)。

为了进一步增强光谱信息, 采用多种方法进行光谱反射率的数学变换, 包括光谱一阶导数、 倒数对数(log(1/R))、 连续统去除(continuum remove, CR)、 多元散射校正(multivariate scattering correction, MSC)和变量标准化(standard normalized variable, SNV)。

在航空数据获取时同步进行地面土壤采样, 综合土壤类型、 地形地貌及农作物类型等进行均匀布点。 采样点布置于农田中, 长期处于翻耕状态, 为了保证样品的代表性, 采样深度为地表0~20 cm层位的耕作层土壤, 由4个子样等量均匀混合成1件样品后采用四分法留取样品。 建三江地区采样数量为225个, 海伦地区采样数量为121个。 指标分析方法如下: Cr的全量, 样品经盐酸、 硝酸、 氢氟酸、 高氯酸分解后用盐酸提取, 采用OPTIMA 2100DV等离子体发射光谱仪测量; 有机质(SOM)采用硫酸、 重铬酸钾消解法滴定; 全氮, 样品采用石墨消解仪消解后用凯氏定氮仪测量; P, K₂O, SiO₂, Al₂O₃, Fe₂O₃, CaO, MgO, Na₂O采用X射线荧光光谱法测量; pH值采用PHS-25酸度计测定。 数据统计分析在SPSS 22软件中完成。

建模采用偏最小二乘法。 建模数据准备: 以地面采样点位置为中心, 将对应位置航空影像中3× 3像元范围的光谱反射率进行平均后作为该点的光谱数据。 建三江研究区共整理225组样本数据, 其中145组作为建模样本, 80组作为验证样本。 海伦研究区共整理121组样本数据, 其中80组样本为建模样本, 剩余41组样本为验证样本。 建模样本和验证样本均按照空间分布均匀挑选。 模型评价指标包括决定系数(R²)和均方根误差(RMSE)。

土壤理化参数SOM, N, P, K₂O, SiO₂, Al₂O₃, Fe₂O₃, CaO, MgO, Na₂O和pH采用光谱数据直接进行建模。 Cr的反演模型自变量包括光谱和理化参数数据, 以探讨间接反演模型的效果。

2.1.1 土壤指标含量描述统计

分别统计了海伦和建三江地面样品中Cr、 养分、 pH和相关主量成分的含量(表1)。 参照规范(GB15618— 2018)要求, 两个地区所有样品的Cr含量均没有超过土壤污染风险筛选值, 表明目前两个地区土壤Cr污染的风险较小。

表1

Table 1

表1(Table 1)

表1 土壤指标含量统计 Table 1 Statistics of soil index content 指标 海伦 建三江 最小值 最大值 均值 变异系数 最小值 最大值 均值 变异系数 Cr/(mg· kg-1) 43.9 88.1 62.76 0.112 37.5 80.8 59.852 0.126 SOM/% 2.138 7.93 5.131 0.193 0.52 9.26 4.099 0.327 N/% 0.107 0.357 0.235 0.155 0.056 0.497 0.221 0.312 P/(mg· kg-1) 0.51 1.355 0.926 0.126 0.322 5.481 0.802 0.251 K2O/% 2.26 2.69 2.502 0.030 2.13 2.88 2.509 0.051 pH 5.15 8.45 6.366 0.091 4.7 7.91 6.126 0.080 SiO2/% 58.26 70.2 64.349 0.022 57.63 72.93 66.001 0.042 Al2O3/% 11.73 15.05 14.002 0.031 11.85 16.47 13.475 0.084 Fe2O3/% 3.27 5.73 4.645 0.062 2.53 5.77 4.154 0.155 MgO/% 0.84 1.8 1.303 0.099 0.59 1.36 0.888 0.167 CaO/% 0.89 4.19 1.637 0.246 0.86 3.65 1.209 0.222 Na2O/% 1.37 1.89 1.641 0.077 1.31 2.34 1.842 0.105

表1 土壤指标含量统计 Table 1 Statistics of soil index content

土壤养分方面, 参照规范(DZT0295— 2016), 两个地区土壤有机质、 全氮、 全磷、 全钾整体上均处于较丰富级和丰富级, 显示出土壤养分条件较好。 海伦地区的有机质、 全氮和全磷的均值高于建三江, 其中有机质更为明显, 这与海伦地区黑土资源更为丰富有关。

土壤其他主量成分方面, 海伦地区土壤Fe₂O₃, MgO和CaO含量整体高于建三江, Na₂O略低于建三江, SiO₂和Fe₂O₃含量两区基本持平。

根据变异程度分级标准^[5], 海伦地区土壤Cr, SOM, N, P和CaO属于中等变异, 其余指标为弱变异; 建三江地区Cr, SOM, N, P, Fe₂O₃, MgO, CaO和Na₂O属于中等变异, 其余为弱变异。 除了Cr, pH和CaO外, 其余指标在建三江地区的变异系数均明显高于海伦地区, 表明在建三江地区这些指标的空间分布更为不均匀。

2.1.2 Cr与土壤理化指标的相关性

分别计算建三江和海伦研究区土壤Cr与理化指标的相关系数(表2), 可以发现两个地区相关系数具有较强的规律性, 表现为: Cr与Al₂O₃, Fe₂O₃, MgO, K₂O和pH均呈极显著的正相关, 与SiO₂, Na₂O和SOM均呈极显著的负相关, 与CaO没有显著相关性。 在个别参数上的相关性略有差异, 表现为: 建三江地区Cr与P呈显著的负相关, 海伦地区Cr与P没有显著相关性; 海伦地区Cr与N呈极显著的负相关, 建三江地区则没有显著的相关性。 总体上两个地区Cr与土壤理化指标的相关性规律具有较高的一致性, 虽然与N和P的相关系数显著性有差异, 但相关系数绝对值均不高。

表2

Table 2

表2(Table 2)

表2 Cr与土壤理化指标的相关系数 Table 2 Correlation coefficient between Cr and soil physical and chemical indexes Cr 建三江 海伦 SiO2 -0.691* * -0.494* * Al2O3 0.865* * 0.742* * Fe2O3 0.775* * 0.570* * CaO -0.024 -0.046 MgO 0.810* * 0.674* * K2O 0.233* * 0.497* * Na2O -0.724* * -0.566* * P -0.162* -0.025 SOM -0.200* * -0.396* * pH 0.196* * 0.287* * N 0.011 -0.160* * 注: * * 代表0.01水平的相关性, * 代表0.05水平的相关性 Note: correlations are significant at 0.01 level (* * ) and 0.05 level (* )

表2 Cr与土壤理化指标的相关系数 Table 2 Correlation coefficient between Cr and soil physical and chemical indexes

2.1.3 Cr的赋存规律分析

通过Cr的相关性进一步分析其赋存规律。 pH是影响土壤Cr赋存状态的重要影响因素, 其对Cr(Ⅲ )和Cr(Ⅵ )的影响机制不同。 土壤中的Cr(Ⅲ )主要以阳离子形式存在, Cr(Ⅲ )进入土壤后容易形成铬和铁氢氧化物的混合物, 吸附在土壤表面形成稳定的沉淀物。 土壤对Cr(Ⅲ )的吸附主要是因为形成Cr(OH)₃的沉淀作用, 其次是阳离子交换吸附作用。 陈英旭等^[6]实验结果表明, pH值4~6属于Cr(Ⅲ )的吸附-沉淀区间, 在此pH值区间内Cr(Ⅲ )开始沉淀, 土壤对Cr(Ⅲ )的吸附能力随pH升高而增强, pH大于6时Cr(Ⅲ )进入稳定沉淀区间。 而土壤对Cr(Ⅵ )的吸附量整体上随pH值增高而减少。 建三江和海伦的土壤pH整体处于中性和偏酸性, 根据上述pH对Cr(Ⅲ )和Cr(Ⅵ )的影响作用, pH与土壤Cr含量整体应呈正相关, 这与相关系数计算结果是较为吻合的。

土壤金属氧化物是吸附Cr的重要载体。 当pH> 4时, 水合铁氧化物表面的负电荷开始增加, 对Cr(Ⅲ )阳离子的吸附能力随着pH的升高而增强, 此外在pH 4~6区间, Cr(Ⅲ )还易与Fe³⁺和Al³⁺发生吸附共沉淀。 土壤游离氧化铁是吸附Cr(Ⅵ )的主要成分, 此外铁、 铝氧化物对Cr(Ⅵ )的吸附能力明显高于粘土矿物^[7]。 Vega^[8]的研究显示土壤水铝矿、 针铁矿和云母对Cr有较强的吸附性。 建三江地区Cr与Al₂O₃, Fe₂O₃和MgO的相关系数分别达到了0.865, 0.775和0.810, 海伦地区对应的相关系数分别为0.742, 0.570和0.674, 显示了在两个研究区上述金属氧化物对Cr均具有较高的吸附性, 其中建三江地区相关性更高。

此外有机质对土壤Cr的形态和含量也有重要影响。 有机质对Cr(Ⅵ )有较强的还原作用, 有机质中的腐殖质能够增加土壤的表面积和表面活性, 从而增加土壤的吸附能力^[9], 促进碳酸盐结合态Cr向有机结合态Cr的转化^[10]。 土壤有机酸与Cr离子会产生络合作用, 从而对Cr在土壤中的吸附解吸产生影响^[7], 陈英旭实验^[6]显示柠檬酸和富啡酸降低土壤对Cr(Ⅲ )的吸附和沉淀作用, 而胡敏酸的相关作用不明显。 建三江地区Cr与有机质的相关系数为-0.2, 海伦地区为-0.396, 虽然相关系数绝对值并不高, 但两个地区Cr与有机质的显著性均达到了极显著水平, 表明在有机质的综合作用下, 两个地区Cr与有机质整体上均表现为负相关。

综上, 基于Cr在土壤中的赋存规律以及相关性分析, 认为两个研究区Cr的赋存规律有较高的一致性, Cr均与Al₂O₃, Fe₂O₃, MgO, K₂O和pH呈极显著的正相关关系, 与SiO₂, Na₂O和SOM呈极显著的负相关关系, 为Cr的建模提供了相关依据。

2.2.1 Cr与不同光谱预处理的相关性

计算Cr与光谱反射率(R)、 反射率变量标准化(SNV)、 反射率多元散射校正(MSC)、 连续统去除(CR)、 反射率一阶导数(R')、 反射率倒数对数(log(1/R))共六种光谱预处理的相关系数。

从图2可看到, 建三江地区反射率R的相关系数曲线较为平缓, 随着波段增大曲线缓慢下降并保持平稳, 整体呈负相关关系, 在2 210 nm处有微弱的峰值(相关系数-0.275); SNV和MSC曲线较为吻合, 在1 800 nm之前呈正相关关系, 之后转为负相关关系, MSC曲线较为明显的峰值为0.424(660 nm), 0.486(1 520 nm), -0.495(2 210 nm), SNV曲线较明显的峰值为0.401(670 nm), 0.326(1 535 nm), -0.492(2 210 nm); CR曲线在1 370 nm之前呈持续下降的趋势, 在2 000 nm之后整体呈负相关关系, 较明显的峰值有0.405(540 nm), -0.590(2 195 nm), -0.479(2 345 nm); R'曲线在1 200 nm之前呈缓慢下降的趋势, 1 415 nm之后波动较大, 较明显的峰值有-0.494(2 105 nm), 0.374(2 390 nm); log(1/R)曲线整体形态与原始反射率曲线呈镜像分布, 相关性也基本等同。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 建三江地区土壤Cr与不同光谱预处理的相关系数曲线Fig.2 Correlation coefficient curves of Cr and different spectral pretreatments in Jiansanjiang area

从图3可见, 海伦地区反射率R的相关系数曲线整体随波段增大而缓慢下降, 形态较为平缓, 相关系数最大为0.165(500 nm); SNV和MSC曲线整体上同样较为吻合, SNV曲线较明显的峰值为0.071(920 nm), -0.089(1 520 nm), -0.086(2 195 nm), MSC曲线较明显的峰值为0.066(940 nm), -0.080(2195 nm); CR曲线在1 100 nm之前波动较小, 之后波动幅度变大, 较明显的峰值有-0.213(1 310 nm), 0.220(1 700 nm), -0.206(2 165 nm), -0.205(2 270 nm), -0.200(2 345 nm); R'曲线在1 000 nm之前呈缓慢下降趋势, 波动幅度较小, 之后波动较为剧烈, 较明显的峰值有-0.270(1 280 nm), -0.183(2 045 nm), 0.225(2 360 nm); log(1/R)曲线整体形态与原始反射率曲线呈镜像分布, 相关性略低于反射率曲线。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 海伦地区土壤Cr与不同光谱预处理的相关系数曲线Fig.3 Correlation coefficient curves of Cr and different spectral pretreatments in Hailun area

2.2.2 不同变换形式对Cr与光谱相关性的影响

对比上述六种光谱预处理的相关性计算结果可以看到, 反射率的相关性曲线整体较为平缓, 没有明显的峰值出现, 这是由于土壤是复杂的混合体, 光谱的混合程度较高。

MSC和SNV对应的曲线相关特征得到增强, 出现了不同程度的峰值。 MSC是在标准光谱的参考下对每条光谱的基线平移和偏移进行修正, 消除散射影响, 增强与成分含量相关的光谱吸收信息; SNV是各波段间反射率数据的标准化, 在一定程度上也可消除光谱中线性平移的影响, 与多元散射校正有类似之处。 因此通过MSC和SNV变换能在一定程度上增强土壤光谱吸收特征。

CR曲线整体上在1 000 nm之前相对平缓, 之后波动幅度较大, 这是由于在1 000 nm之前土壤光谱主要受有机质和铁锰氧化物的影响, 影响结果是光谱整体反射率的影响, 诊断特征相对较弱; 1 000 nm之后光谱主要受黏土矿物、 水等影响, 诊断特征相对明显。 因此经CR变换后在1 000 nm之后, 尤其是2 000 nm之后黏土矿物吸收区间出现了明显的峰值。

R'曲线整体形态特征和CR曲线有相似之处, 大致在1 000 nm之前较平缓, 之后波动加剧, 因为两种变换方式都是为了增强局部特征。 但是R'曲线的局部波动幅度明显更高, 与其计算方式有关, R'变换在增强局部吸收特征的同时也增强了噪声信息。

此次log(1/R)变换没有明显增强Cr与光谱的相关性。

2.2.3 不同地区Cr的光谱特征对比

在整体相关性水平上, 每种光谱形式下建三江地区Cr与光谱的相关性均高于海伦地区。 建三江地区正相关值最高为0.486(MSC变换), 负相关值最小为-0.590(CR变换), 海伦地区相应的值分别为0.220(CR变换)、 -0.270(R'变换)。 建三江地区Cr与光谱相关性更高的原因之一认为是该区Cr与土壤矿物、 有机质等组分的相关性更强, 间接影响到了光谱的相关性。

虽然两区土壤Cr与光谱的相关程度存在差异, 但是特征波段的位置存在相似性。 对比图2和图3可以看到: 两区反射率曲线均没有明显的相关峰值; SNV和MSC曲线在2 200 nm附近均存在明显的负相关峰值; CR曲线在2 345 nm处均存在明显的负相关峰值; R'曲线在2 360 nm附近存在明显的正相关峰值。 两区特征波段位置的相似性在一定程度上表明了土壤Cr赋存的控制因素存在一致性。

2.2.4 Cr的特征波段筛选

为了增强Cr反演模型的精度及可迁移性, 需要综合两个研究区的光谱特征筛选特征波段。 根据上述相关性分析结果筛选每种光谱预处理对应的特征波段, 结果见表3。

表3

Table 3

表3(Table 3)

表3 Cr的特征波段筛选 Table 3 Screening result of characteristic band of Cr 光谱预处理 波段/nm R 全波段 SNV 1 520, 2 195, 2 210, 2 225 MSC 2 195, 2 210, 2 225 CR 2 345 R' 2 360, 1 275, 2 390, 2 405 log(1/R) 全波段

表3 Cr的特征波段筛选 Table 3 Screening result of characteristic band of Cr

2.3.1 不同光谱预处理建模结果分析

利用不同光谱预处理进行建模及验证, 从结果(表4)可以看到, 建三江和海伦研究区中SNV对应的建模R²和验证R²均为最高。 整体上对于两个研究区, SNV和MSC对应的建模R²均明显高于其他四种光谱预处理, RMSE也明显低于其他四种光谱预处理, 验证结果的趋势整体上与建模结果较为吻合。 综合上文中不同光谱变换形式对Cr与光谱间相关性的影响结果, 认为本研究中SNV和MSC变换能够较好地提升Cr与光谱的相关性, 使建模结果得到了稳定的提升。 CR和R'变换在某些波段范围内能大幅提升微弱的光谱特征, 但同时也引入了噪声干扰, 建模结果中只有建三江地区的CR对应的模型R²略高于原始反射率, 其余均有一定的降低。 log(1/R)对应的模型精度基本与原始反射率持平。

表4

Table 4

表4(Table 4)

表4 不同光谱形式的建模精度 Table 4 Modeling accuracy of different spectral pretreatments 光谱形式 建三江 海伦 建模 验证 建模 验证 RMSE R2 RMSE R2 RMSE R2 RMSE R2 R 5.36 0.536 5.56 0.483 5.83 0.446 6.03 0.386 SNV 4.76 0.643 5.11 0.571 5.34 0.537 5.54 0.471 MSC 4.83 0.631 5.23 0.556 5.41 0.521 5.61 0.456 CR 5.31 0.544 5.58 0.487 5.88 0.433 6.05 0.377 R' 5.46 0.515 5.65 0.445 5.86 0.441 6.08 0.381 log(1/R) 5.41 0.533 5.59 0.479 5.81 0.449 6.06 0.383

表4 不同光谱形式的建模精度 Table 4 Modeling accuracy of different spectral pretreatments

2.3.2 间接反演模型

间接反演模型的自变量集成了光谱数据和土壤组分数据, 为了后续模型迁移性研究的需要, 将两个研究区的模型自变量设置为完全一致。 根据上述Cr的光谱特征分析与光谱建模结果, 光谱数据采用SNV变换后的特征波段数据(表3); 根据Cr在土壤中的赋存规律分析, 参与建模的土壤组分数据设定为Al₂O₃, Fe₂O₃, MgO, K₂O, pH, SiO₂, Na₂O和SOM。 建三江研究区建立的模型为: Cr=-0.123b_{1 520}-0.098b_{2 195}-0.096b_{2 210}-0.094b_{2 225}-0.705SiO₂+3.932Al₂O₃+0.715Fe₂O₃+0.358MgO+0.253K₂O-0.312Na₂O-0.601SOM-1.498pH+61.813, 建模结果显示R²为0.751, RMSE为3.741; 模型验证结果显示R²为0.687, RMSE为3.873。 海伦研究区建立的模型为: Cr=-0.053b_{1 520}+0.031b_{2 195}+0.033b_{2 210}+0.035b_{2 225}-0.808SiO₂+8.235Al₂O₃+3.972Fe₂O₃+1.945MgO+1.357K₂O-1.634Na₂O-0.399SOM+0.024pH-21.11, 建模结果显示R²为0.676, RMSE为4.133; 模型验证结果显示R²为0.643, RMSE为4.316。

通过对比两个研究区Cr的间接反演模型与纯光谱模型的模型精度, 可以发现间接反演模型明显提升了Cr的模型精度。 以纯光谱模型中精度最高的SNV结果(表4)为比较对象, 建三江研究区的间接反演模型将建模R²由0.643提升到了0.751, 验证R²由0.571提升到了0.687; 海伦研究区的间接反演模型将建模R²由0.537提升到了0.676, 验证R²由0.471提升到了0.643。 间接反演模型包括了更多Cr的赋存规律信息, 因此模型精度有明显提升。 在航空高光谱反演制图时, 首先利用光谱信息进行土壤相关组分的含量反演, 然后将相关组分的反演数据纳入间接反演模型进行Cr的含量反演。 此次建三江研究区Cr间接反演模型中相关组分的建模R²均在0.783以上, 最高为0.887(Fe₂O₃); 海伦研究区相关组分的建模R²均在0.736以上, 最高为0.836(MgO), 均达到了较好的反演效果, 能够为Cr的间接反演模型提供数据支持。

2.3.3 模型迁移性分析

分别将两个研究区各自建立的SNV光谱模型和间接模型进行互相迁移, 以探讨直接与间接反演模型在迁移性上的差异及迁移的可行性。 从图4(a)可看到, 海伦研究区的直接光谱模型应用到建三江研究区后, 实测值和预测值之间几乎没有相关性, R²仅为0.019 4, 并且预测值出现了若干负值, 可以认为模型不具有迁移的价值。 而图4(b)的结果与图4(a)类似, 同样不具有可迁移性。 出现此类现象的原因是不同地区间土壤光谱信息的影响因素有较多差异, 而包括Cr在内的重金属元素在可见光-近红外区间又没有明显的诊断光谱特征, 从而影响了不同地区间模型的通用性。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 SNV光谱模型迁移结果Fig.4 Migration results of SNV spectral model

图5显示, 将海伦研究区Cr的间接反演模型应用到建三江研究区后, 实测与预测值的回归R²达到0.597 5, 而建三江的间接反演模型应用到海伦后, 回归R²为0.577 3。 可以看到间接模型迁移后预测值和实测值的线性关系较为明显, 相比纯光谱模型, 间接模型在两个研究区间的可迁移性得到了明显提升。 这是由于两个地区中Cr与模型中土壤相关组分的相关性是一致的, 同时相关组分表现的光谱特征具有一致性, 提升了模型迁移后的稳健性。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 间接反演模型迁移结果Fig.5 Migration results of indirect inversion model