孫志超,祁雨薇,汪東川,姜 杰,韓明利,趙人杰

1 天津綠茵景觀生態(tài)建設股份有限公司,天津 300384

2 天津城建大學地質(zhì)與測繪學院,天津 300384

3 天津城建大學,天津市土木建筑結構防護與加固重點實驗室,天津 300384

土壤鹽堿化是中國半干旱地區(qū)面臨的重要資源與生態(tài)環(huán)境問題,制約著土地的可持續(xù)利用與農(nóng)業(yè)的穩(wěn)定高效發(fā)展[1—3]。鹽堿地作為農(nóng)業(yè)方面重要的潛在后備資源,對其改良后利用是響應“國家高標準農(nóng)田建設”,實現(xiàn)“2022年建成10億畝高標準農(nóng)田”宏偉目標的重要舉措。目前,鹽堿地土壤改良技術已經(jīng)發(fā)展成水利工程、化學生物、農(nóng)藝調(diào)理相結合的綜合改良技術,且已在不同地區(qū)獲得較好的生態(tài)、經(jīng)濟及社會效益[4—5]。鹽堿地土壤改良項目完成后,存在3a或5a的改良技術效果追蹤,最常用、最基本的方法為測定土壤含鹽量變化直接反映技術效果,其中常用的含鹽量測定方法為飽合泥漿法、電導率烘干法等[6]。然而土壤含鹽量測定需要現(xiàn)場采集大量的土壤樣本并檢測,費時、費力、成本高,且鹽堿地土壤改良的最終目的是增加植被生態(tài)效益,即提升植被生態(tài)服務功能或增加農(nóng)作物產(chǎn)量。相較于土壤取樣檢測方法,基于遙感技術監(jiān)測植被的生態(tài)效益具有便捷高效、全域成像、動態(tài)可視、成本較低等優(yōu)勢,更適用于生態(tài)修復項目的改良技術效果追蹤。

扶松林研究表明,時空尺度上單位面積糧食產(chǎn)量與植被總初級生產(chǎn)力(GPP)之間顯著正相關[7],且GPP反映了植被的固碳生態(tài)效益[8],故通過分析植被GPP變化量反映工程修復的生態(tài)效益切實可行。另外,植被生態(tài)效益與植被長勢顯著正相關[9—11],可通過分析MODIS、Landsat、Sentinel等遙感影像的植被光譜差異,構建歸一化植被指數(shù)(NDVI)、增強型植被指數(shù)(EVI)等反映植被長勢[12—13],且已在華北平原[14]、青藏高原[15]、新疆[16]等地區(qū)得到應用,同樣適用于小麥、玉米等農(nóng)作物長勢監(jiān)測[17—19]。目前,空間分辨率最高且能夠直接提取GPP數(shù)據(jù)的遙感數(shù)據(jù)為MOD17A2H產(chǎn)品,已普遍應用于流域、城市群、省域等大尺度地區(qū)的生態(tài)效益研究中[20—23]。此GPP數(shù)據(jù)的空間分辨率為500m,然而生態(tài)修復工程項目的面積一般較小,明顯不適用于項目改良技術效果的區(qū)域差異分析,故必須進行空間尺度轉換[24],提高GPP數(shù)據(jù)的空間分辨率。成方妍等研究表明Landsat NDVI 與凈初級生產(chǎn)力(NPP)的線性關系顯著[8],錢婭等研究表明基于Landsat數(shù)據(jù)計算的吸收性光合有效輻射(APAR)可以有效響應玉米實際 GPP 值的季節(jié)性波動[25],李振旺等研究利用多源多尺度衛(wèi)星遙感數(shù)據(jù)(環(huán)境、哨兵-2、Landsat、MODIS、GOME等)識別了中國北方草原區(qū)GPP分布特征[26]。以上研究表明,基于Landsat或Sentinel數(shù)據(jù)的植被指數(shù)能夠?qū)崿F(xiàn)GPP數(shù)據(jù)的空間分辨率提高至30m或10m,屬于兩種遙感數(shù)據(jù)的融合方式。若在Sentinel數(shù)據(jù)、GPP數(shù)據(jù)的空間尺度轉換中,增加Landsat數(shù)據(jù),基于逐層次、逐像元的思路,依舊能夠?qū)崿F(xiàn)GPP數(shù)據(jù)的空間分辨率提高至10m,而且利用兩種空間分辨率植被指數(shù)的信息融合,可以增加GPP數(shù)據(jù)的穩(wěn)定性和準確性,減輕數(shù)據(jù)間分辨率差距大對尺度轉換精度的影響程度,是多源數(shù)據(jù)間尺度轉換的一種新思路。

本研究以內(nèi)蒙古巴彥淖爾市趙貴圪旦組鹽堿地為研究區(qū),以植被凈生產(chǎn)力反演模型(CASA)中應用的NDVI、比值植被指數(shù)(SR)為基礎,利用30m空間分辨率的Landsat-8 OLI影像與10m空間分辨率的Sentinel-2號兩種遙感影像,基于移動窗口法原理,首先在邊長500m的窗口內(nèi)將空間分辨率為500m的GPP數(shù)據(jù)修正至30m,然后在邊長30m的窗口內(nèi)將空間分辨率為30m的GPP數(shù)據(jù)修正至10m,完成逐層逐像元修正后,應用熱點分析方法探究鹽堿地治理效果的空間異質(zhì)性,同時利用轉移矩陣方法對土壤深度0—20cm、20—40cm、40—60cm的含鹽量變化進行了分析,以期反映鹽堿地改良技術帶來的全局變化效果,為種植區(qū)域合理布置與重點治理區(qū)域精準定位提供參考,同時該技術適用于其他生態(tài)修復工程的生態(tài)效益評估。

1 研究區(qū)概況

研究區(qū)屬于2018年千億斤糧食增產(chǎn)工程“改鹽增草(飼)興牧”示范項目,位于內(nèi)蒙古巴彥淖爾市烏拉特前旗烏拉山鎮(zhèn)蓿亥村趙貴圪旦組,面積約163.87hm2?；谘芯繀^(qū)內(nèi)的排鹽渠溝,將項目區(qū)劃分為39個地塊,如圖1所示。研究區(qū)氣候?qū)儆谥袦貛Т箨懶约撅L氣候,雨熱同期,晝夜溫差大。水道屬于黃河水系,距離黃河主干道僅3km,引黃灌溉水的入滲是地下水的主要來源,且平均地下水位在0.5m左右。地貌類型為黃河沖積平原,土壤以黃河泛濫沉積物為主,土壤類型主要為重壤和壤粘土,土壤鹽漬化現(xiàn)象嚴重,受土地鹽漬化影響,堿蓬、堿草等鹽生植被零星分布[27],種植農(nóng)作物主要為食葵,食葵年均產(chǎn)量在750kg/hm2至3750kg/hm2范圍內(nèi)。

圖1 研究區(qū)行政邊界Fig.1 Administrative boundaries of the study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源與預處理

2.1.1土壤數(shù)據(jù)

根據(jù)《鹽堿化耕地普查技術標準》對研究區(qū)進行土壤普查,以多點取樣檢測方法為主。土壤取樣點如圖2所示,以每點控制0.1km2土地面積為原則進行網(wǎng)格式分布布點16個,利用內(nèi)徑為3cm的土鉆對每個取樣點分三次、分三層0—20cm、20—40cm、40—60cm進行取土樣,2018年9月與2020年8月分別取樣48個,基于飽合泥漿法測定土壤含鹽量后,利用克里金空間插值方法生成研究區(qū)土壤含鹽量分布數(shù)據(jù),數(shù)據(jù)空間分辨率為10m,見圖3。由于該項目運用的土壤改良技術符合暗管排鹽技術基本原理,故本研究參考《中國人民共和國土地管理行業(yè)標準(2013)》,暗管改良鹽堿地技術規(guī)程第一部分:土壤調(diào)查(標準代號:TD/T1043.1—2013)中規(guī)定的土壤鹽化分級指標,將研究區(qū)土壤含鹽量分為5個等級[6],如表1。同時,此標準明確了該分級指標適用于濱海、半濕潤、半干旱、干旱區(qū)。

表1 土壤含鹽量等級Table 1 Classification of soil salt content

表2 顯著性熱點/冷點區(qū)域分類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植被數(shù)據(jù)

植被數(shù)據(jù)分為NDVI數(shù)據(jù)與GPP數(shù)據(jù),NDVI數(shù)據(jù)來自美國Earth Date search網(wǎng)站(https://search. earthdata.nasa.gov/)下載的Landsat-8 OLI數(shù)據(jù)與從歐空局Open Hub網(wǎng)站(https://scihub.copernicus.eu/)下載的Sentinel-2號數(shù)據(jù),Landsat-8 OLI數(shù)據(jù)日期為2018年8月3日、2020年9月7日,Sentinel-2號數(shù)據(jù)日期為2018年8月20日、2020年8月9日。GPP數(shù)據(jù)來源于從美國Earth Date search網(wǎng)站下載的MOD17A2H數(shù)據(jù),日期為2018年8月30日、2020年8月29日,以上影像無云層干擾。將以上三種數(shù)據(jù)采用統(tǒng)一的空間坐標系,同時將空間分辨率統(tǒng)一重采樣為10m。

2.2 方法

2.2.1土壤鹽化等級變化分析方法

轉移矩陣方法最初適用于反映某一區(qū)域不同地類面積在兩個不同時間節(jié)點內(nèi)相互轉化的動態(tài)過程信息[28],本研究中基于此方法探討2018年與2020年五類土壤鹽化等級面積相互轉換的動態(tài)信息,轉移矩陣公式[29]為:

(1)

式中,S表示面積;i與j分別表示轉移前后的土壤鹽化等級;Sij表示轉移前第i類鹽化等級向轉移后第j類鹽化等級轉移的面積;n表示土壤鹽化等級總數(shù),本研究中n為5。

2.2.2GPP數(shù)據(jù)空間分辨率提高方法

光合有效輻射比率(FPAR)是衡量植被冠層對光合有效輻射能量吸收能力的一個重要參數(shù), 其估算精度直接影響陸地生態(tài)系統(tǒng)GPP的估算精度[30—32]。CASA模型是基于遙感技術估算FPAR的主要方法,其原理是利用NDVI與比值植被指數(shù)(SR)構建關系模型進行估算[33]。為精準反映研究區(qū)內(nèi)的植被總初級生產(chǎn)力,本文基于Landsat-8 OLI數(shù)據(jù)與Sentinel-2數(shù)據(jù)依次構建GPP 30m與10m 空間分辨率的FPAR調(diào)節(jié)因子,實現(xiàn)對MODIS 500m空間分辨率GPP數(shù)據(jù)的空間分辨率提高,以期更精確反映研究區(qū)內(nèi)植被總初級生產(chǎn)力的空間差異及其動態(tài)變化特征。技術路線圖見圖4。

圖4 技術路線圖Fig.4 Methodological workflow

基于Landsat-8 OLI數(shù)據(jù)與Sentinel-2數(shù)據(jù)分別計算NDVI與SR,計算公式如下:

(2)

(3)

式中,BNIR與BR依次代表遙感影像的近紅外波段與紅波段。

基于以上數(shù)據(jù)利用以下公式分別計算Landsat-8與Sentinel-2數(shù)據(jù)的FPAR[32]。

(4)

式中,FPAR是光合有效輻射,NDVImax與 NDVImin依次為研究區(qū)植被的NDVI最大值與最小值,SRmax與SRmin依次為研究區(qū)植被的SR最大值與最小值。

依據(jù)圖5,利用以下公式分別構建基于Landsat-8與Sentinel-2植被指數(shù)數(shù)據(jù)的GPP修正因子,修正因子在Arcgis 10.2軟件中采用Phthon語言二次開發(fā)計算所得。

圖5 三種數(shù)據(jù)空間分辨率關系示意圖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陸地標準產(chǎn)品分辨率(500m×500m)窗口內(nèi)對應基于Landsat-8生成修正因子值,FPARLmean是MOD17A2H陸地標準產(chǎn)品分辨率(500m×500m)窗口內(nèi)對應基于Landsat-8生成的修正因子平均值。FPARS是Landsat-8影像分辨率(30m×30m)窗口內(nèi)對應基于Sentinel-2生成的修正因子值,FPARSmean是Landsat-8影像分辨率(30m×30m)窗口內(nèi)對應基于Sentinel-2生成的修正因子平均值。

利用GPP修正公式對MOD17A2H陸地標準產(chǎn)品中提取的GPP數(shù)據(jù)進行修正。

GPPA=GPPAS=AS×GPPAL

(9)

GPPAL=AL×GPP

(10)

式中,GPPA為修正后的GPP數(shù)據(jù),空間分辨率為10m;GPPAS與GPPAL分明為利用Sentinel-2與Landsat-8數(shù)據(jù)修正后的GPP數(shù)據(jù)。

2.2.3GPP變化分析方法

為探討豎井排鹽工程等改良技術對農(nóng)作物長勢的影響效果,本研究利用公式11分析工程施工前后GPP的變化量(GPP Change,GC),公式為:

GC=GPP施工后-GPP施工前

(11)

(12)

(13)

3 結果與分析

3.1 土壤含鹽量變化結果分析

從圖3、表3中,發(fā)現(xiàn)土壤改良效果最明顯的區(qū)域為12—20號地塊,且不同土壤深度中不同鹽化等級轉入轉出面積存在差異,土壤含鹽量改良效果隨深度增加逐漸降低,這種現(xiàn)象與豎井排鹽工程導致地下水位下降,使得上層土壤中無機鹽溶于水后,隨水位逐漸下降至深層土壤有關。具體表現(xiàn)為:土壤深度0—20cm,強度鹽土、鹽土的面積為轉出狀態(tài),轉出面積共約64.69hm2;而土壤深度40—60cm,強度鹽土、鹽土的面積為轉入狀態(tài),轉入面積共約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中發(fā)現(xiàn),土壤深度20—40cm,中度鹽土的轉入量最大,面積約43.21hm2,強度鹽土與鹽土的轉入量次之,分別約19hm2、10.48hm2。其中,中度鹽土面積增加的主要來源是強度鹽土,面積約20.96m2,其次為輕度鹽土,面積約15.68hm2。強度鹽土與鹽土轉移面積的主要來源分別是中度鹽土、強度鹽土,轉移量分別約15.79hm2、8.81hm2。表6中發(fā)現(xiàn),土壤深度40—60cm,含鹽量較低的非鹽土、輕度鹽土、中度鹽土均為轉出狀態(tài),強度鹽土與鹽土均為轉入狀態(tài)。其中中度鹽土面積轉出量最大,約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 植被總初級生產(chǎn)力變化結果分析

3.2.1植被總初級生產(chǎn)力變化量分析

由圖6得,修正后的GPP數(shù)據(jù)在項目區(qū)內(nèi)的空間異質(zhì)性明顯,可以準確判斷不同地塊鹽堿地土壤改良工程的修復效果。為進一步探討研究區(qū)植被長勢的變化規(guī)律,該部分將未種植農(nóng)作物的地塊剔除,統(tǒng)計分析種植區(qū)域的GPP變化量?；诖?利用分位數(shù)法將2018年GPP數(shù)據(jù)與GPP提升/降低區(qū)域數(shù)據(jù)分別分為5級,GPP提升區(qū)域從Ⅰ級至Ⅴ級變化強度逐漸增加,而GPP降低區(qū)域從Ⅰ級至Ⅴ級變化強度逐漸降低,分級標準見表7。

表7 數(shù)據(jù)等級表/(g C/m2)Table 7 Table of data grading

圖6 GPP數(shù)據(jù)修正前后對比圖Fig.6 Comparison of GPP data before and after correction

表8、表9表明,研究區(qū)內(nèi)GPP提升區(qū)域面積明顯大于其他區(qū)域面積,面積約多60.36hm2,說明豎井排鹽工程治理效果顯著。2018年GPP數(shù)據(jù)的前四類,其提升區(qū)域面積明顯大于降低區(qū)域面積,而其第五類數(shù)據(jù)的提升區(qū)域面積小于降低區(qū)域面積,這說明鹽堿地改良技術對GPP低值區(qū)域的改良效果優(yōu)于GPP高值區(qū)域。

表8 GPP提升區(qū)域變化強度分級統(tǒng)計表/hm2Table 8 Hierarchical statistical table of Change intensity in the lifting area of GPP

表9 GPP降低區(qū)域變化強度分級統(tǒng)計表/hm2Table 9 Hierarchical statistical table of Change intensity in the lowering area of GPP

在2018年GPP數(shù)據(jù)的Ⅰ、Ⅱ等級,提升區(qū)域的Ⅳ、Ⅴ等級的面積稍大于Ⅰ、Ⅱ、Ⅲ等級,降低區(qū)域的Ⅳ、Ⅴ等級的面積稍大于Ⅰ、Ⅱ等級,這說明在GPP數(shù)據(jù)低值區(qū)域,提升區(qū)域變化強度大的面積要大于變化強度小的面積,而降低區(qū)域變化強度大的面積要小于變化強度小的面積。

在2018年GPP數(shù)據(jù)的Ⅳ、Ⅴ等級,提升區(qū)域變化強度從Ⅰ級至Ⅴ級,面積量逐漸減少,降低區(qū)域變化強度從Ⅰ級至Ⅴ級,面積量同樣有逐漸降低趨勢,這說明在GPP數(shù)據(jù)高值區(qū)域,提升區(qū)域變化強度大的面積要小于變化強度小的面積,而降低區(qū)域變化強度大的面積要大于變化強度小的面積。

在2018年GPP數(shù)據(jù)的Ⅲ等級,提升區(qū)域變化強度從Ⅰ級至Ⅴ級,面積程遞增趨勢,降低區(qū)域變化強度Ⅰ級、Ⅴ級的面積存在明顯小于其他等級面積的趨勢,這說明在GPP數(shù)據(jù)中值區(qū)域,提升區(qū)域變化強度大的面積要大于變化強度小的面積,而降低區(qū)域變化強度最小與最大的面積要小于變化強度中等的面積。

3.2.2植被總初級生產(chǎn)力變化量的空間異質(zhì)性分析

本研究基于熱點分析方法得到GPP提升與降低區(qū)域內(nèi)變化程度的空間分布圖,以期說明GPP變化強度的空間分布差異。圖7表明,提升區(qū)域的冷點區(qū)域分布較零散,在4—9號、18—24號地塊面積占比較大,說明此區(qū)域鹽堿地改良技術對植被GPP的提升效果不明顯;而提升區(qū)域的熱點區(qū)域集中分布在12—17、25—27、34號地塊,說明工程技術對此區(qū)域的植被GPP提高有明顯改善效果;除4、28、30、36號地塊不存在提升區(qū)域中的其他區(qū)域,其他地塊均含有,說明工程技術的改良效應在以上區(qū)域內(nèi)表現(xiàn)不強烈。降低區(qū)域的其他區(qū)域多分布在地塊之間的邊界區(qū)域,說明工程技術改良效果較差,反作用程度一般;冷點區(qū)域集中分布在6、12—14、22—23、27—28、30—31號地塊內(nèi),說明這些區(qū)域內(nèi)工程技術對植被GPP有明顯的降低作用,更說明改良效果最差,有嚴重的反作用;而其熱點區(qū)域多分布在3—4、18—19號地塊內(nèi),說明工程技術改良效果較差,反作用程度較輕。研究結果將為次年該區(qū)域農(nóng)作物種植提供參考,即重點關注降低區(qū)域內(nèi)植被的生長狀態(tài)。

圖7 GPP提升/降低區(qū)域變化強度分布圖Fig.7 Distribution map of change intensity in the lifting/lowering area of GPP

4 結論與討論

通過分析趙貴圪旦組鹽堿地改良前后土壤含鹽量與農(nóng)作物GPP變化情況,主要結論包括:

(1)12—20號地塊為土壤改良效果的顯著區(qū)域,且土壤深度0—20cm含鹽量的降低效果顯著,其中鹽土面積減少61.69hm2。表層土壤鹽分下滲,造成土壤深度20—40cm與40—60cm非鹽土、輕度鹽土、中度鹽土的面積轉出量增加,而強度鹽土與鹽土的面積轉入量增加顯著。

(2)農(nóng)作物種植區(qū)域內(nèi),GPP提高區(qū)域面積占比70%以上,且GPP提高區(qū)域約為降低區(qū)域面積的2.5倍,說明鹽堿地改良技術對70%以上的區(qū)域存在正向作用,且該技術對2018年GPP低值區(qū)域的改良效果優(yōu)于GPP高值區(qū)域。

(3)基于移動窗口法+逐層逐像元的修正原理,修正后的GPP數(shù)據(jù)能夠精準反映GPP變化的空間差異,具體表現(xiàn)為改良效果在農(nóng)作物種植的東部區(qū)域優(yōu)于西部區(qū)域,在最北部區(qū)域效果最差,且在6、12—14、22—23、27—28、30—31號地塊內(nèi),GPP降低效果顯著。

本研究通過構建逐層逐像元修正因子,將Landsat-8 OLI影像、Sentinel-2號影像、MODIS GPP數(shù)據(jù)融合,將GPP數(shù)據(jù)的空間分辨率從500m提高至10m?；诖?融合空間分辨率為亞米級的高分數(shù)據(jù),空間分辨率將為亞米級,將有助于進一步精確定位土壤改良效果較差區(qū)域,有利于生態(tài)效益的空間異質(zhì)性分析。另外,通過參考本文構建的GPP數(shù)據(jù)融合技術,不僅能夠?qū)崿F(xiàn)GPP數(shù)據(jù)的空間尺度轉換,而且為其他類型數(shù)據(jù)融合提供思路,從而實現(xiàn)多源數(shù)據(jù)的空間尺度轉換,故將該技術推廣至其他數(shù)據(jù)集的融合需進一步通過實驗驗證效果。

為獲取較高空間分辨率的GPP數(shù)據(jù),最基本方法為利用較高空間分辨率的GPP計算模型因子修正原始GPP數(shù)據(jù),僅采用一種數(shù)據(jù)修正GPP數(shù)據(jù)的研究較多[24],本研究引用Landsat、Sentinel兩種數(shù)據(jù),依次將GPP的空間分辨率提高至30m、10m,增加逐層降尺度過程,符合數(shù)據(jù)融合的優(yōu)勢,增加GPP數(shù)據(jù)的穩(wěn)定性和準確性。然而,本研究未揭示兩種數(shù)據(jù)修正結果與一種數(shù)據(jù)修正結果間的差異,將需要依據(jù)不同場景量化證明。