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

王海瀟，王勇輝,2

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

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

瑪納斯湖濕地是干旱區典型的濕地，而土壤重金屬嚴重地危害了濕地的生態安全。為了更好地提高重金屬反演精度，本研究以Landsat8 OLI數據為基礎，利用DEM數據進行地形校正，并且對校正后的反射率進行倒數、導數、對數等數學變換，從而分別建立8種重金屬最優的PLSR預測模型，并對土壤重金屬的空間分布及生態風險進行評價，進而更好地利用遙感數據對土壤重金屬進行預測。

1 材料與方法

1.1 研究區概況

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

1.2 數據來源及數據處理

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

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

1.3 PLSR模型的建立與驗證

PLSR集成了主成分分析、典型相關分析和普通多元線性回歸3種方法的優點，克服了自變量之間的多重線性相關和樣本數量小于波長變量的問題，使構建的模型更穩定，有助于多元數據統計分析[16-17]。

1.4 生態風險預警評估方法

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

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

Rapant等[21]給出了相應生態風險劃分標準，參考相關研究提出的生態風險程度描述[23-24]，將生態風險預警指數值與警度進行了概念關聯[25]，如表1所示。

表1土壤生態風險預警綜合判別標準

Table 1 Standard of soil ecological risk warning

風險等級Risk level風險指數Risk index預警類型Warning type風險程度描述Description of risk level1IER≤0無警No warning生態系統服務功能基本完整,生態環境基本未受干擾,生態系統結構完整,功能性強,系統恢復再生能力強,生態問題不顯著,生態災害少。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 結果與分析

2.1 反射率的地形校正

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

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

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

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波段及決定系數Band and coefficientof determinationCrCuMnNiPbZnHgAsR波段 Band21434124決定系數 R20.410.520.410.240.510.420.430.441/R波段 Band11312243決定系數 R20.450.590.510.550.480.450.470.56LogR波段 Band21431356決定系數 R20.420.570.430.460.350.520.490.51R'波段 Band24132512決定系數 R20.460.470.530.480.460.480.460.47

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

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

2.4 土壤重金屬空間分布格局及生態風險預警評估

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

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

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

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

通過比較基于遙感的土壤重金屬生態風險預警與基于實測值的土壤重金屬生態風險預警可知，兩者的分析結果較為一致，結果表明可以用遙感的手段來反演該研究區的重金屬分布情況。

表5 研究區重金屬生態風險分級評價

表6 研究區不同分區重金屬生態風險分級評價

3 結論與討論

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

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

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