基于地形校正后Landsat 8的土壤重金屬反演研究

王海瀟，王勇輝,2

(1.新疆師范大學(xué)地理科學(xué)與旅游學(xué)院，新疆 烏魯木齊 830054；2.新疆干旱區(qū)湖泊環(huán)境與資源實驗室，新疆 烏魯木齊 830054)

土壤中的重金屬滯留時間較長而且很難降解，因而易在植物體內(nèi)累積，進而通過整個食物鏈在人體中富集，造成潛在的健康風(fēng)險[1-2]。隨著社會經(jīng)濟的快速發(fā)展，人為活動也成為土壤污染最主要的原因，工業(yè)生產(chǎn)、交通運輸、礦業(yè)開采、農(nóng)藥濫用等均為增加土壤重金屬含量的主要途徑[3-4]。近幾年來，國內(nèi)外學(xué)者對土壤重金屬污染開展了廣泛研究[5-6]。傳統(tǒng)的重金屬原位監(jiān)測技術(shù)精度較高，但需耗費大量的人力物力，無法適用于大范圍重金屬污染監(jiān)測需求，而遙感技術(shù)具有快速、宏觀的特性，表現(xiàn)出其他方法無法比擬的優(yōu)點。如Malley和Williams利用湖泊沉積物的反射光譜快速預(yù)測重金屬含量[7]；Ren等[8]對礦區(qū)農(nóng)用地中的反射光譜與重金屬進行研究，發(fā)現(xiàn)反射光譜對As和Cu元素有一定的指示性；吳昀昭等[9]應(yīng)用農(nóng)田土壤光譜反射率估算土壤表層Hg元素的濃度；李巨寶等[10]采用偏最小二乘方法，構(gòu)建基于土壤反射光譜的Se，N，F(xiàn)e含量預(yù)測模型。地表反射率表征地球表面對太陽輻射的反射能力，是揭示地物目標(biāo)本質(zhì)的最有用信息[11]，而不同地區(qū)的像元由于受太陽、大氣和地形等多種因素的影響，接收到的太陽輻射有很大的差別，造成有些影像區(qū)處于陰影覆蓋下，而另一些卻處于過度感光狀態(tài)[12]，導(dǎo)致了“同物異譜”和“異物同譜”現(xiàn)象的出現(xiàn)，這對地物信息提取造成了不利影響。地形校正可以消除由地形起伏而引起的像元反射率的變化，使像元反射率更好地反映地物光譜特性[13]。目前利用地形校正模型消除地形影響，已成為目前高光譜遙感圖像步入應(yīng)用前的有效處理手段，郭云開等[14]和穆悅等[15]通過地形校正對地表反射率進行校正，并取得較好結(jié)果。因此本研究在Landsat 8 OLI數(shù)據(jù)的基礎(chǔ)上引入地形校正，旨在提高研究區(qū)重金屬反演精度。

瑪納斯湖濕地是干旱區(qū)典型的濕地，而土壤重金屬嚴(yán)重地危害了濕地的生態(tài)安全。為了更好地提高重金屬反演精度，本研究以Landsat8 OLI數(shù)據(jù)為基礎(chǔ)，利用DEM數(shù)據(jù)進行地形校正，并且對校正后的反射率進行倒數(shù)、導(dǎo)數(shù)、對數(shù)等數(shù)學(xué)變換，從而分別建立8種重金屬最優(yōu)的PLSR預(yù)測模型，并對土壤重金屬的空間分布及生態(tài)風(fēng)險進行評價，進而更好地利用遙感數(shù)據(jù)對土壤重金屬進行預(yù)測。

1 材料與方法

1.1 研究區(qū)概況

瑪納斯湖濕地位于歐亞腹地新疆維吾爾自治區(qū)昌吉回族自治州境內(nèi)，其所處位置屬于溫帶大陸性干旱氣候區(qū)，年平均氣溫6℃，年降水量100～150 mm，年蒸發(fā)量3 500 mm左右，是準(zhǔn)噶爾盆底西部的一個構(gòu)造沉降中心，是瑪納斯河地表徑流的尾閭，湖體受瑪納斯河季節(jié)性水量補給變化而游移變動?，敿{斯湖是一個咸水湖，湖區(qū)位于坐標(biāo)86.0695°～86.3088°E，44.3741°～44.5196°N，湖的南北跨度有28.62 km，東西跨度21.40 km，湖區(qū)面積約550 km2，平均深度約6 m。在氣候、地形和水文主要因素作用下，瑪納斯?jié)竦氐貐^(qū)形成以湖泊沼澤、沖洪積平原、湖積平原為主的地質(zhì)地貌類型。

1.2 數(shù)據(jù)來源及數(shù)據(jù)處理

于2015年11月10日至15日在瑪納斯湖濕地區(qū)域內(nèi)，選擇了3個有代表性的區(qū)域，分別是湖泊入湖口(A區(qū))、退化區(qū)(B區(qū))和恢復(fù)區(qū)(C區(qū))，每個區(qū)布設(shè)9個樣點，共布設(shè)27個采樣點，每個采樣點采集0～20 cm土層土壤樣本2個，共采集到54個土壤樣本。將采集的樣品在室內(nèi)風(fēng)干，剔除枯枝落葉、根莖、石子、動物殘體等雜質(zhì)，碾碎，用100目尼龍篩過篩后供測試用。土壤樣品中As 和 Hg采用原子熒光法(AFS)測定；重金屬Cu、Mn、Cr、Pb、Co、Ni、Zn采用離子體光譜(ICP)測定，將每個采樣點的測定數(shù)據(jù)取平均值。為了與野外實測數(shù)據(jù)相結(jié)合，選擇的影像數(shù)據(jù)為野外采樣時間段內(nèi)的Landsat 8 OLI遙感數(shù)據(jù)，獲取時間為2015年11月15日，DEM數(shù)據(jù)的來源為中國科學(xué)院數(shù)據(jù)中心的ASTER GDEM 數(shù)據(jù)，分辨率為30 m。

地形校正法是目前應(yīng)用最為廣泛的C校正模型，前期數(shù)據(jù)處理包括幾何精校正和COST模型大氣校正，幾何精校正誤差(RMS)小于0.5，然后用矢量邊界進行裁剪。利用DEM數(shù)據(jù)在ENVI軟件中計算得到坡度、坡向數(shù)據(jù)，同時為了消除大氣散射和地表相鄰點反射光折射造成的像元亮度值(DN值)和太陽入射角α之間的余弦關(guān)系，本研究在二階校正模型的基礎(chǔ)上利用IDL(Interactive data language)語言進行改進，并經(jīng)過二次校正才能得到校正結(jié)果。

1.3 PLSR模型的建立與驗證

PLSR集成了主成分分析、典型相關(guān)分析和普通多元線性回歸3種方法的優(yōu)點，克服了自變量之間的多重線性相關(guān)和樣本數(shù)量小于波長變量的問題，使構(gòu)建的模型更穩(wěn)定，有助于多元數(shù)據(jù)統(tǒng)計分析[16-17]。

1.4 生態(tài)風(fēng)險預(yù)警評估方法

生態(tài)風(fēng)險預(yù)警評估源于生態(tài)風(fēng)險評價，更強調(diào)對生態(tài)系統(tǒng)可能存在風(fēng)險的警示研究，狹義預(yù)警僅指對自然資源或生態(tài)風(fēng)險可能出現(xiàn)的衰竭或危機而建立的報警，而廣義預(yù)警則涵蓋了生態(tài)風(fēng)險的維護及防止危機發(fā)展的過程[19-20]。對研究區(qū)土壤重金屬進行生態(tài)風(fēng)險預(yù)警評估采用Rapant等[21]提出的生態(tài)風(fēng)險指數(shù)法，屬狹義預(yù)警，其表達式為:

式中，IERi表示超過臨界限量的第i種重金屬生態(tài)風(fēng)險指數(shù),CAi表示第i種重金屬的實測含量(mg·kg-1),CRi表示第i種重金屬的臨界限量(mg·kg-1),IER表示待測樣品生態(tài)風(fēng)險指數(shù)。其中CRi采用新疆土壤背景值[22]。

Rapant等[21]給出了相應(yīng)生態(tài)風(fēng)險劃分標(biāo)準(zhǔn)，參考相關(guān)研究提出的生態(tài)風(fēng)險程度描述[23-24]，將生態(tài)風(fēng)險預(yù)警指數(shù)值與警度進行了概念關(guān)聯(lián)[25]，如表1所示。

表1土壤生態(tài)風(fēng)險預(yù)警綜合判別標(biāo)準(zhǔn)

Table 1 Standard of soil ecological risk warning

風(fēng)險等級Risk level風(fēng)險指數(shù)Risk index預(yù)警類型Warning type風(fēng)險程度描述Description of risk level1IER≤0無警No warning生態(tài)系統(tǒng)服務(wù)功能基本完整,生態(tài)環(huán)境基本未受干擾,生態(tài)系統(tǒng)結(jié)構(gòu)完整,功能性強,系統(tǒng)恢復(fù)再生能力強,生態(tài)問題不顯著,生態(tài)災(zāi)害少。The ecosystem service function is basically complete; the ecological environment is basically un-disturbed; the ecosystem structure is complete; the function is strong; the system has strong re-covery and regeneration ability; the ecological problem is not significant; and the ecological dis-aster is less.20

2 結(jié)果與分析

2.1 反射率的地形校正

瑪納斯湖流域地勢南高北低、溝壑相間的特點決定了局部微地形是影響研究區(qū)表觀反射率、植被、土壤、氣候等各方面地域分異的主要控制因素。由圖1可知，以B1波段反射率和重金屬Cu為例，經(jīng)過地形校正后的B1波段的反射率值與27個實測土壤表層Cu含量的決定系數(shù)高于地形校正前，R2從0.46提高至0.52，擬合結(jié)果表明，利用DEM校正后的反射率反映土壤重金屬Cu的效果較DEM校正前的好。

2.2 土壤反射率(反射率變換形式)與土壤各重金屬含量相關(guān)性分析

將進行地形校正后的各波段反射率及其各種數(shù)學(xué)變換(倒數(shù)1/R、對數(shù)LogR、一階微分R′)與各土壤重金屬含量進行相關(guān)性分析，并且從各種數(shù)學(xué)變換中篩選出與各土壤重金屬含量決定系數(shù)最大的波段，其相關(guān)性均達到0.01顯著水平，結(jié)果如表2所示。從表2可知，相對于地形校正后原始反射率(R)而言，其倒數(shù)(1/R)、對數(shù)LogR與各重金屬(除Pb以外)的相關(guān)性都有所提高，而一階微分R′與Cu和Pb含量的相關(guān)性有所降低，與其余6種重金屬的相關(guān)性也有所提高。

2.3 土壤重金屬反演模型建立與檢驗

圖1 土壤重金屬Cu實測值和地形校正前(a)后(b)B1反射率擬合效果對比Fig.1 Relationships between Cu content and reflectance of B1 before (a) and after (b) topographic correction

變換形式Transformation form波段及決定系數(shù)Band and coefficientof determinationCrCuMnNiPbZnHgAsR波段 Band21434124決定系數(shù) R20.410.520.410.240.510.420.430.441/R波段 Band11312243決定系數(shù) R20.450.590.510.550.480.450.470.56LogR波段 Band21431356決定系數(shù) R20.420.570.430.460.350.520.490.51R'波段 Band24132512決定系數(shù) R20.460.470.530.480.460.480.460.47

表3 基于單波段的土壤各重金屬預(yù)測模型

表4 基于多波段的土壤各重金屬PLSR預(yù)測模型

2.4 土壤重金屬空間分布格局及生態(tài)風(fēng)險預(yù)警評估

2.4.1 土壤重金屬風(fēng)險空間分布 利用表4中的各重金屬的預(yù)測模型分別反演各重金屬空間分布，同時根據(jù)表1將IER值按其大小分為5個等級以表示不同程度的風(fēng)險程度，并利用ARCGIS制作了研究區(qū)2015年11月的土壤重金屬風(fēng)險等級分布圖，如圖2所示。

通過圖2可知，研究區(qū)土壤重金屬風(fēng)險等級總體上呈現(xiàn)從西南方向至東北方向逐漸減弱的趨勢，并且研究區(qū)土壤重金屬污染總體上處于輕警以上級別，生態(tài)服務(wù)功能已開始退化，生態(tài)環(huán)境受到一定破壞，生態(tài)系統(tǒng)結(jié)構(gòu)有所變化，但尚可維持基本功能，受干擾后易惡化，生態(tài)問題顯現(xiàn)，應(yīng)該加強對該地區(qū)的重金屬污染的治理。

圖2 土壤重金屬風(fēng)險等級分布Fig.2 Distribution map of risk grade of soil heavy metals

2.4.2 基于實測數(shù)據(jù)生態(tài)風(fēng)險預(yù)警評估 為了驗證土壤重金屬遙感反演的精度，通過研究區(qū)采集的實測樣點實測數(shù)據(jù)進行驗證，并且通過對實測數(shù)據(jù)進行不同分區(qū)(A、B、C三個區(qū))來分析研究區(qū)重金屬生態(tài)風(fēng)險，基于實測數(shù)據(jù)的土壤重金屬風(fēng)險評價如表5、表6所示。由表5可知，從單項指數(shù)看Cr、Cu、Mn、Ni、Pb、Zn、Hg、As處于無警級的概率分別為92.59%、77.78%、96.30%、3.70%、11.11%、48.15%、3.70%、25.93%，處于預(yù)警級的概率分別為7.41%、7.41%、3.70%、96.30%、33.33%、40.74%、3.70%、74.07%，Cu、Pb、Zn和Hg處于輕警級的概率分別為11.11%、55.56%、3.70%、3.70%，其他元素均未出現(xiàn)輕警級的樣點，處于中警級的重金屬元素主要有Cu和Hg，其概率分別為3.70%、37.04%，處于重警級的重金屬元素主要有Zn和Hg，其概率分別為4.41%、51.85%。其單項生態(tài)風(fēng)險指數(shù)處于無警以上級別的風(fēng)險概率大小順序為Hg=Ni>Pb>As>Zn>Cu>Cr>Mn。從綜合指數(shù)看研究區(qū)11.11%的樣點處在無警級別，屬于最低生態(tài)風(fēng)險，有33.33%的樣點處于預(yù)警級風(fēng)險，而樣點落在輕警、中警、重警級風(fēng)險的概率分別為25.93%、18.52%、11.11%，由表5、表6可知，屬于重警級污染的樣點主要分布在恢復(fù)區(qū)(C區(qū))，并且研究區(qū)土壤重金屬生態(tài)風(fēng)險排序為：恢復(fù)區(qū)(C區(qū))>退化區(qū)(B區(qū))>湖泊入湖口(A區(qū))。

通過比較基于遙感的土壤重金屬生態(tài)風(fēng)險預(yù)警與基于實測值的土壤重金屬生態(tài)風(fēng)險預(yù)警可知，兩者的分析結(jié)果較為一致，結(jié)果表明可以用遙感的手段來反演該研究區(qū)的重金屬分布情況。

表5 研究區(qū)重金屬生態(tài)風(fēng)險分級評價

表6 研究區(qū)不同分區(qū)重金屬生態(tài)風(fēng)險分級評價

3 結(jié)論與討論

目前已有不少學(xué)者利用遙感技術(shù)對土壤重金屬污染監(jiān)測進行研究探索，并取得一定的成果[26]，但是很少有研究會考慮去除地形對反射率的影響。本研究在Landsat 8的基礎(chǔ)上，引入DEM數(shù)據(jù)對反射率進行地形校正，經(jīng)過地形校正后的B1波段的反射率值與27個實測土壤表層Cu含量的決定系數(shù)高于地形校正前，R2從0.46提高至0.52，擬合結(jié)果表明，利用DEM校正后的反射率反映土壤重金屬Cu的效果較DEM校正前的好。諸多學(xué)者的研究表明，將反射率進行適當(dāng)?shù)臄?shù)學(xué)變換可以有效抑制遙感影像中噪聲對地物光譜的影響，改善反射率與土壤某些元素的相關(guān)性，提高模型反演精度[27]，因此本研究在地形校正后原始反射率的基礎(chǔ)上進行倒數(shù)、對數(shù)和導(dǎo)數(shù)等變換，研究發(fā)現(xiàn)相對于地形校正后原始反射率(R)而言，其倒數(shù)(1/R)、對數(shù)LogR與各重金屬(除Pb以外)的相關(guān)性都有所提高，這與郭云開等[28]的研究結(jié)果較為一致。

在重金屬的反演模型中，基于多波段的PLSR預(yù)測模型整體上相對于單波段要好，同時利用土壤各重金屬的最佳預(yù)測模型分別反演相應(yīng)的土壤重金屬含量，并引入土壤重金屬生態(tài)風(fēng)險指數(shù)用于評價研究區(qū)的土壤重金屬生態(tài)風(fēng)險并進行預(yù)警，研究發(fā)現(xiàn)研究區(qū)土壤重金屬風(fēng)險等級總體上呈現(xiàn)從西南方向至東北方向逐漸減弱的趨勢，土壤重金屬生態(tài)風(fēng)險排序為恢復(fù)區(qū)(C區(qū))>退化區(qū)(B區(qū))>湖泊入湖口(A區(qū))，這與實地驗證情況較為一致，原因可能為C區(qū)為恢復(fù)區(qū)，靠近克拉瑪依市，其受到克拉瑪依市生活廢水、工業(yè)生產(chǎn)廢水、農(nóng)業(yè)污水灌溉的影響，使得該區(qū)的土壤重金屬風(fēng)險高于A區(qū)和B區(qū)，而A區(qū)由于為湖泊如湖口，受到湖水對土壤的沖刷作用，使得土壤重金屬風(fēng)險最低。研究區(qū)土壤重金屬主要來源可能是來自工業(yè)生產(chǎn)、生活用水中的重金屬成分，這與劉娜等[29]、馬輝英等[30]的研究結(jié)論基本一致。

為了驗證基于遙感的土壤重金屬生態(tài)風(fēng)險預(yù)警的預(yù)測精度，將研究區(qū)土壤重金屬含量實測數(shù)據(jù)也通過重金屬生態(tài)風(fēng)險指數(shù)進行計算，通過比較發(fā)現(xiàn)兩者的分析結(jié)果較為一致，結(jié)果表明可以用遙感的手段來反演該研究區(qū)的重金屬分布情況，同時研究區(qū)土壤重金屬污染總體上處于輕警以上級別，生態(tài)服務(wù)功能已開始退化，應(yīng)該加強該地區(qū)的重金屬污染治理工作。