孙志超,祁雨薇,汪东川,姜 杰,韩明利,赵人杰

1 天津绿茵景观生态建设股份有限公司,天津 300384

2 天津城建大学地质与测绘学院,天津 300384

3 天津城建大学,天津市土木建筑结构防护与加固重点实验室,天津 300384

土壤盐碱化是中国半干旱地区面临的重要资源与生态环境问题,制约着土地的可持续利用与农业的稳定高效发展[1—3]。盐碱地作为农业方面重要的潜在后备资源,对其改良后利用是响应“国家高标准农田建设”,实现“2022年建成10亿亩高标准农田”宏伟目标的重要举措。目前,盐碱地土壤改良技术已经发展成水利工程、化学生物、农艺调理相结合的综合改良技术,且已在不同地区获得较好的生态、经济及社会效益[4—5]。盐碱地土壤改良项目完成后,存在3a或5a的改良技术效果追踪,最常用、最基本的方法为测定土壤含盐量变化直接反映技术效果,其中常用的含盐量测定方法为饱合泥浆法、电导率烘干法等[6]。然而土壤含盐量测定需要现场采集大量的土壤样本并检测,费时、费力、成本高,且盐碱地土壤改良的最终目的是增加植被生态效益,即提升植被生态服务功能或增加农作物产量。相较于土壤取样检测方法,基于遥感技术监测植被的生态效益具有便捷高效、全域成像、动态可视、成本较低等优势,更适用于生态修复项目的改良技术效果追踪。

扶松林研究表明,时空尺度上单位面积粮食产量与植被总初级生产力(GPP)之间显著正相关[7],且GPP反映了植被的固碳生态效益[8],故通过分析植被GPP变化量反映工程修复的生态效益切实可行。另外,植被生态效益与植被长势显著正相关[9—11],可通过分析MODIS、Landsat、Sentinel等遥感影像的植被光谱差异,构建归一化植被指数(NDVI)、增强型植被指数(EVI)等反映植被长势[12—13],且已在华北平原[14]、青藏高原[15]、新疆[16]等地区得到应用,同样适用于小麦、玉米等农作物长势监测[17—19]。目前,空间分辨率最高且能够直接提取GPP数据的遥感数据为MOD17A2H产品,已普遍应用于流域、城市群、省域等大尺度地区的生态效益研究中[20—23]。此GPP数据的空间分辨率为500m,然而生态修复工程项目的面积一般较小,明显不适用于项目改良技术效果的区域差异分析,故必须进行空间尺度转换[24],提高GPP数据的空间分辨率。成方妍等研究表明Landsat NDVI 与净初级生产力(NPP)的线性关系显著[8],钱娅等研究表明基于Landsat数据计算的吸收性光合有效辐射(APAR)可以有效响应玉米实际 GPP 值的季节性波动[25],李振旺等研究利用多源多尺度卫星遥感数据(环境、哨兵-2、Landsat、MODIS、GOME等)识别了中国北方草原区GPP分布特征[26]。以上研究表明,基于Landsat或Sentinel数据的植被指数能够实现GPP数据的空间分辨率提高至30m或10m,属于两种遥感数据的融合方式。若在Sentinel数据、GPP数据的空间尺度转换中,增加Landsat数据,基于逐层次、逐像元的思路,依旧能够实现GPP数据的空间分辨率提高至10m,而且利用两种空间分辨率植被指数的信息融合,可以增加GPP数据的稳定性和准确性,减轻数据间分辨率差距大对尺度转换精度的影响程度,是多源数据间尺度转换的一种新思路。

本研究以内蒙古巴彦淖尔市赵贵圪旦组盐碱地为研究区,以植被净生产力反演模型(CASA)中应用的NDVI、比值植被指数(SR)为基础,利用30m空间分辨率的Landsat-8 OLI影像与10m空间分辨率的Sentinel-2号两种遥感影像,基于移动窗口法原理,首先在边长500m的窗口内将空间分辨率为500m的GPP数据修正至30m,然后在边长30m的窗口内将空间分辨率为30m的GPP数据修正至10m,完成逐层逐像元修正后,应用热点分析方法探究盐碱地治理效果的空间异质性,同时利用转移矩阵方法对土壤深度0—20cm、20—40cm、40—60cm的含盐量变化进行了分析,以期反映盐碱地改良技术带来的全局变化效果,为种植区域合理布置与重点治理区域精准定位提供参考,同时该技术适用于其他生态修复工程的生态效益评估。

1 研究区概况

研究区属于2018年千亿斤粮食增产工程“改盐增草(饲)兴牧”示范项目,位于内蒙古巴彦淖尔市乌拉特前旗乌拉山镇蓿亥村赵贵圪旦组,面积约163.87hm2。基于研究区内的排盐渠沟,将项目区划分为39个地块,如图1所示。研究区气候属于中温带大陆性季风气候,雨热同期,昼夜温差大。水道属于黄河水系,距离黄河主干道仅3km,引黄灌溉水的入渗是地下水的主要来源,且平均地下水位在0.5m左右。地貌类型为黄河冲积平原,土壤以黄河泛滥沉积物为主,土壤类型主要为重壤和壤粘土,土壤盐渍化现象严重,受土地盐渍化影响,碱蓬、碱草等盐生植被零星分布[27],种植农作物主要为食葵,食葵年均产量在750kg/hm2至3750kg/hm2范围内。

图1 研究区行政边界Fig.1 Administrative boundaries of the study area

2 数据与方法

2.1 数据来源与预处理

2.1.1土壤数据

根据《盐碱化耕地普查技术标准》对研究区进行土壤普查,以多点取样检测方法为主。土壤取样点如图2所示,以每点控制0.1km2土地面积为原则进行网格式分布布点16个,利用内径为3cm的土钻对每个取样点分三次、分三层0—20cm、20—40cm、40—60cm进行取土样,2018年9月与2020年8月分别取样48个,基于饱合泥浆法测定土壤含盐量后,利用克里金空间插值方法生成研究区土壤含盐量分布数据,数据空间分辨率为10m,见图3。由于该项目运用的土壤改良技术符合暗管排盐技术基本原理,故本研究参考《中国人民共和国土地管理行业标准(2013)》,暗管改良盐碱地技术规程第一部分:土壤调查(标准代号:TD/T1043.1—2013)中规定的土壤盐化分级指标,将研究区土壤含盐量分为5个等级[6],如表1。同时,此标准明确了该分级指标适用于滨海、半湿润、半干旱、干旱区。

表1 土壤含盐量等级Table 1 Classification of soil salt content

表2 显著性热点/冷点区域分类Table 2 Classification of significant hot/cold spots

图3 土壤深度0—60cm盐化等级空间分布图Fig.3 Spatial distribution map of salinization grade in soil depth 0—60cm

2.1.2植被数据

植被数据分为NDVI数据与GPP数据,NDVI数据来自美国Earth Date search网站(https://search. earthdata.nasa.gov/)下载的Landsat-8 OLI数据与从欧空局Open Hub网站(https://scihub.copernicus.eu/)下载的Sentinel-2号数据,Landsat-8 OLI数据日期为2018年8月3日、2020年9月7日,Sentinel-2号数据日期为2018年8月20日、2020年8月9日。GPP数据来源于从美国Earth Date search网站下载的MOD17A2H数据,日期为2018年8月30日、2020年8月29日,以上影像无云层干扰。将以上三种数据采用统一的空间坐标系,同时将空间分辨率统一重采样为10m。

2.2 方法

2.2.1土壤盐化等级变化分析方法

转移矩阵方法最初适用于反映某一区域不同地类面积在两个不同时间节点内相互转化的动态过程信息[28],本研究中基于此方法探讨2018年与2020年五类土壤盐化等级面积相互转换的动态信息,转移矩阵公式[29]为:

(1)

式中,S表示面积;i与j分别表示转移前后的土壤盐化等级;Sij表示转移前第i类盐化等级向转移后第j类盐化等级转移的面积;n表示土壤盐化等级总数,本研究中n为5。

2.2.2GPP数据空间分辨率提高方法

光合有效辐射比率(FPAR)是衡量植被冠层对光合有效辐射能量吸收能力的一个重要参数, 其估算精度直接影响陆地生态系统GPP的估算精度[30—32]。CASA模型是基于遥感技术估算FPAR的主要方法,其原理是利用NDVI与比值植被指数(SR)构建关系模型进行估算[33]。为精准反映研究区内的植被总初级生产力,本文基于Landsat-8 OLI数据与Sentinel-2数据依次构建GPP 30m与10m 空间分辨率的FPAR调节因子,实现对MODIS 500m空间分辨率GPP数据的空间分辨率提高,以期更精确反映研究区内植被总初级生产力的空间差异及其动态变化特征。技术路线图见图4。

图4 技术路线图Fig.4 Methodological workflow

基于Landsat-8 OLI数据与Sentinel-2数据分别计算NDVI与SR,计算公式如下:

(2)

(3)

式中,BNIR与BR依次代表遥感影像的近红外波段与红波段。

基于以上数据利用以下公式分别计算Landsat-8与Sentinel-2数据的FPAR[32]。

(4)

式中,FPAR是光合有效辐射,NDVImax与 NDVImin依次为研究区植被的NDVI最大值与最小值,SRmax与SRmin依次为研究区植被的SR最大值与最小值。

依据图5,利用以下公式分别构建基于Landsat-8与Sentinel-2植被指数数据的GPP修正因子,修正因子在Arcgis 10.2软件中采用Phthon语言二次开发计算所得。

图5 三种数据空间分辨率关系示意图Fig.5 Schematic diagram of the relationship between three kinds of spatial resolution data

(5)

(6)

(7)

(8)

式中,AL是基于Landsat-8影像生成的GPP修正因子,AS是基于Sentinel-2影像生成的GPP修正因子。FPARL是MOD17A2H陆地标准产品分辨率(500m×500m)窗口内对应基于Landsat-8生成修正因子值,FPARLmean是MOD17A2H陆地标准产品分辨率(500m×500m)窗口内对应基于Landsat-8生成的修正因子平均值。FPARS是Landsat-8影像分辨率(30m×30m)窗口内对应基于Sentinel-2生成的修正因子值,FPARSmean是Landsat-8影像分辨率(30m×30m)窗口内对应基于Sentinel-2生成的修正因子平均值。

利用GPP修正公式对MOD17A2H陆地标准产品中提取的GPP数据进行修正。

GPPA=GPPAS=AS×GPPAL

(9)

GPPAL=AL×GPP

(10)

式中,GPPA为修正后的GPP数据,空间分辨率为10m;GPPAS与GPPAL分明为利用Sentinel-2与Landsat-8数据修正后的GPP数据。

2.2.3GPP变化分析方法

为探讨竖井排盐工程等改良技术对农作物长势的影响效果,本研究利用公式11分析工程施工前后GPP的变化量(GPP Change,GC),公式为:

GC=GPP施工后-GPP施工前

(11)

(12)

(13)

3 结果与分析

3.1 土壤含盐量变化结果分析

从图3、表3中,发现土壤改良效果最明显的区域为12—20号地块,且不同土壤深度中不同盐化等级转入转出面积存在差异,土壤含盐量改良效果随深度增加逐渐降低,这种现象与竖井排盐工程导致地下水位下降,使得上层土壤中无机盐溶于水后,随水位逐渐下降至深层土壤有关。具体表现为:土壤深度0—20cm,强度盐土、盐土的面积为转出状态,转出面积共约64.69hm2;而土壤深度40—60cm,强度盐土、盐土的面积为转入状态,转入面积共约30.76hm2。

表3 土壤深度0—60cm盐化等级转移面积表/hm2Table 3 Table of salinization grade transfer area in soil depth 0—60cm

表4中知,土壤深度0—20cm,中度盐土面积大量增加,盐土面积大量减少,其中中度盐土面积增加的主要来源是盐土,面积约47.48hm2,其次为强度盐土,面积约14.77hm2。表5中发现,土壤深度20—40cm,中度盐土的转入量最大,面积约43.21hm2,强度盐土与盐土的转入量次之,分别约19hm2、10.48hm2。其中,中度盐土面积增加的主要来源是强度盐土,面积约20.96m2,其次为轻度盐土,面积约15.68hm2。强度盐土与盐土转移面积的主要来源分别是中度盐土、强度盐土,转移量分别约15.79hm2、8.81hm2。表6中发现,土壤深度40—60cm,含盐量较低的非盐土、轻度盐土、中度盐土均为转出状态,强度盐土与盐土均为转入状态。其中中度盐土面积转出量最大,约20.39hm2,主要转移方向为强度盐土与盐土,转移面积分别约24.35hm2、12.65hm2;轻度盐土的转出量次之,面积约8.37hm2,主要转移方向为中度盐土;强度盐土与盐土的转入面积分别约36.79hm2、21.87hm2。强度盐土的转移来源主要为中度盐土与盐土,面积依次约24.35hm2、12.44hm2;盐土的转移来源主要为中度盐土与强度盐土,面积依次约12.65hm2、9.22hm2。

表4 土壤深度0—20cm盐化等级面积转移矩阵表/hm2Table 4 Transfer matrix table of salinization grade area in soil depth 0—20cm

表5 土壤深度20—40cm盐化等级面积转移矩阵表/hm2Table 5 Transfer matrix table of salinization grade area in soil depth 20—40cm

表6 土壤深度40—60cm盐化等级面积转移矩阵表/hm2Table 6 Transfer matrix table of salinization grade area in soil depth 40—60cm

3.2 植被总初级生产力变化结果分析

3.2.1植被总初级生产力变化量分析

由图6得,修正后的GPP数据在项目区内的空间异质性明显,可以准确判断不同地块盐碱地土壤改良工程的修复效果。为进一步探讨研究区植被长势的变化规律,该部分将未种植农作物的地块剔除,统计分析种植区域的GPP变化量。基于此,利用分位数法将2018年GPP数据与GPP提升/降低区域数据分别分为5级,GPP提升区域从Ⅰ级至Ⅴ级变化强度逐渐增加,而GPP降低区域从Ⅰ级至Ⅴ级变化强度逐渐降低,分级标准见表7。

表7 数据等级表/(g C/m2)Table 7 Table of data grading

图6 GPP数据修正前后对比图Fig.6 Comparison of GPP data before and after correction

表8、表9表明,研究区内GPP提升区域面积明显大于其他区域面积,面积约多60.36hm2,说明竖井排盐工程治理效果显著。2018年GPP数据的前四类,其提升区域面积明显大于降低区域面积,而其第五类数据的提升区域面积小于降低区域面积,这说明盐碱地改良技术对GPP低值区域的改良效果优于GPP高值区域。

表8 GPP提升区域变化强度分级统计表/hm2Table 8 Hierarchical statistical table of Change intensity in the lifting area of GPP

表9 GPP降低区域变化强度分级统计表/hm2Table 9 Hierarchical statistical table of Change intensity in the lowering area of GPP

在2018年GPP数据的Ⅰ、Ⅱ等级,提升区域的Ⅳ、Ⅴ等级的面积稍大于Ⅰ、Ⅱ、Ⅲ等级,降低区域的Ⅳ、Ⅴ等级的面积稍大于Ⅰ、Ⅱ等级,这说明在GPP数据低值区域,提升区域变化强度大的面积要大于变化强度小的面积,而降低区域变化强度大的面积要小于变化强度小的面积。

在2018年GPP数据的Ⅳ、Ⅴ等级,提升区域变化强度从Ⅰ级至Ⅴ级,面积量逐渐减少,降低区域变化强度从Ⅰ级至Ⅴ级,面积量同样有逐渐降低趋势,这说明在GPP数据高值区域,提升区域变化强度大的面积要小于变化强度小的面积,而降低区域变化强度大的面积要大于变化强度小的面积。

在2018年GPP数据的Ⅲ等级,提升区域变化强度从Ⅰ级至Ⅴ级,面积程递增趋势,降低区域变化强度Ⅰ级、Ⅴ级的面积存在明显小于其他等级面积的趋势,这说明在GPP数据中值区域,提升区域变化强度大的面积要大于变化强度小的面积,而降低区域变化强度最小与最大的面积要小于变化强度中等的面积。

3.2.2植被总初级生产力变化量的空间异质性分析

本研究基于热点分析方法得到GPP提升与降低区域内变化程度的空间分布图,以期说明GPP变化强度的空间分布差异。图7表明,提升区域的冷点区域分布较零散,在4—9号、18—24号地块面积占比较大,说明此区域盐碱地改良技术对植被GPP的提升效果不明显;而提升区域的热点区域集中分布在12—17、25—27、34号地块,说明工程技术对此区域的植被GPP提高有明显改善效果;除4、28、30、36号地块不存在提升区域中的其他区域,其他地块均含有,说明工程技术的改良效应在以上区域内表现不强烈。降低区域的其他区域多分布在地块之间的边界区域,说明工程技术改良效果较差,反作用程度一般;冷点区域集中分布在6、12—14、22—23、27—28、30—31号地块内,说明这些区域内工程技术对植被GPP有明显的降低作用,更说明改良效果最差,有严重的反作用;而其热点区域多分布在3—4、18—19号地块内,说明工程技术改良效果较差,反作用程度较轻。研究结果将为次年该区域农作物种植提供参考,即重点关注降低区域内植被的生长状态。

图7 GPP提升/降低区域变化强度分布图Fig.7 Distribution map of change intensity in the lifting/lowering area of GPP

4 结论与讨论

通过分析赵贵圪旦组盐碱地改良前后土壤含盐量与农作物GPP变化情况,主要结论包括:

(1)12—20号地块为土壤改良效果的显著区域,且土壤深度0—20cm含盐量的降低效果显著,其中盐土面积减少61.69hm2。表层土壤盐分下渗,造成土壤深度20—40cm与40—60cm非盐土、轻度盐土、中度盐土的面积转出量增加,而强度盐土与盐土的面积转入量增加显著。

(2)农作物种植区域内,GPP提高区域面积占比70%以上,且GPP提高区域约为降低区域面积的2.5倍,说明盐碱地改良技术对70%以上的区域存在正向作用,且该技术对2018年GPP低值区域的改良效果优于GPP高值区域。

(3)基于移动窗口法+逐层逐像元的修正原理,修正后的GPP数据能够精准反映GPP变化的空间差异,具体表现为改良效果在农作物种植的东部区域优于西部区域,在最北部区域效果最差,且在6、12—14、22—23、27—28、30—31号地块内,GPP降低效果显著。

本研究通过构建逐层逐像元修正因子,将Landsat-8 OLI影像、Sentinel-2号影像、MODIS GPP数据融合,将GPP数据的空间分辨率从500m提高至10m。基于此,融合空间分辨率为亚米级的高分数据,空间分辨率将为亚米级,将有助于进一步精确定位土壤改良效果较差区域,有利于生态效益的空间异质性分析。另外,通过参考本文构建的GPP数据融合技术,不仅能够实现GPP数据的空间尺度转换,而且为其他类型数据融合提供思路,从而实现多源数据的空间尺度转换,故将该技术推广至其他数据集的融合需进一步通过实验验证效果。

为获取较高空间分辨率的GPP数据,最基本方法为利用较高空间分辨率的GPP计算模型因子修正原始GPP数据,仅采用一种数据修正GPP数据的研究较多[24],本研究引用Landsat、Sentinel两种数据,依次将GPP的空间分辨率提高至30m、10m,增加逐层降尺度过程,符合数据融合的优势,增加GPP数据的稳定性和准确性。然而,本研究未揭示两种数据修正结果与一种数据修正结果间的差异,将需要依据不同场景量化证明。