基于HJ1A-HSI的龍口污水灌溉區土壤重金屬含量反演

蔡東全，吳泉源，曹學江，蒙永輝，周歷媛，劉 碩，張龍龍

（1.山東師范大學人口·資源與環境學院，山東　濟南250014；2.山東省煤田地質規劃勘察研究院，山東泰安271000；3.山東省地質環境監測總站，山東濟南250014）

基于HJ1A-HSI的龍口污水灌溉區土壤重金屬含量反演

蔡東全1，吳泉源1，曹學江2，蒙永輝3，周歷媛1，劉 碩1，張龍龍1

（1.山東師范大學人口·資源與環境學院，山東濟南250014；2.山東省煤田地質規劃勘察研究院，山東泰安271000；3.山東省地質環境監測總站，山東濟南250014）

污水灌溉在一定程度上減輕了農業用水問題，但是也造成了重金屬元素的積累。以HSI高光譜遙感影像為數據源，利用地理信息系統和偏最小二乘回歸（PLSR）分析方法，并結合實地采樣檢測的土壤重金屬含量數據，分析了重金屬離子在HSI數據中的光譜特征，建立了土壤重金屬含量與高光譜數據的偏最小二乘回歸模型，實現了對龍口污水灌溉區土壤重金屬含量的反演研究。結果表明:原始光譜數據經過一定的光譜變換后，就建模效果而言，Cu的建模效果最好，判定系數（R2）高達0.965 5，Mn、Ni和Pb的建模效果較好，R2分別為0.927 3、0.910 4 和0.912 7；驗證樣本中，Mn的R2為0.9128，說明該模型對Mn的反演能力最強，As、Cu和Ni的模型反演能力相當，R2分別為0.727 6、0.780 4和0.733 3。基于HJ1A-HSI數據并結合PSLR方法反演龍口污水灌溉區土壤重金屬含量具有方法可行性與結果的可靠性。

土壤重金屬；HSI；偏最小二乘回歸；污水灌溉區

合理的污水灌溉不僅可以緩解一個地區農業用水的緊張局面，而且可以充分利用污水中所含有的作物營養元素，降低農業生產成本，還可以提高水資源的重復利用率，并可作為污水處理回用的一種有效手段［1］。然而隨著灌溉時間的不斷延長和灌溉面積的不斷擴大，這種方式帶來了十分顯著的土壤重金屬積累問題。受到重金屬污染的土壤常常通過農作物和水進入人體，經過富集作用對人類和環境產生危害［2］。

傳統的土壤重金屬監測方法可以針對不同的重金屬元素種類取得較好的測量精度，但是只能對特定的采樣點或剖面進行監測，無法做到實時、動態地對重金屬污染進行監測研究，且效率低下、耗時費力［3-4］。高光譜遙感數據具有更精細的地物光譜曲線，更強的目標識別和信息提取能力［5］，使得以遙感手段定量反演和研究土壤的化學成分成為可能，也為定量反演土壤重金屬含量提供了一條新的思路。近年來，國內外的一些學者嘗試利用高光譜數據，對土壤重金屬含量進行定性與定量研究［6-8］，但以我國環境減災衛星HSI高光譜為數據源，定量反演土壤重金屬含量的研究則相對較少。因此，本文利用地理信息系統和偏最小二乘回歸方法，研究龍口污水灌溉地區土壤幾種主要重金屬鹽分參數與HSI高光譜影像的光譜反射率之間的關系，為快速、動態和大面積定量反演污水灌溉區土壤重金屬含量尋求一種新的方法。

1　研究區概況

龍口市地處膠東半島的西北部，地理坐標為東經120°13′20″～120°44′45″，北緯37°27′45″～37°47′30″。全市地勢南高北低，呈階梯形下降，東部和南部為低山丘陵地帶，西部和北部為濱海平原。受成土母質、地形、地貌、氣候和水文條件的影響，該地區形成了棕壤、褐土、潮土和姜黑土四大土壤類型。研究區內有效灌溉面積超過8 000 hm2，其中污水灌溉總面積約4 500 hm2，主要灌溉水源為龍口、黃城兩座污水處理廠排放的工業廢水和生活廢水。

2　數據獲取與預處理

2.1土壤數據的獲取

本研究以Landsat8數據7、4、1波段合成的遙感圖像為參考圖，并利用GPS定位技術于2013年11月5—15日在研究區內進行野外調查和采點，選擇32個100 m×100 m樣方區（見圖1），并要求樣方區在各種土地利用類型中盡可能規則分布，且樣方區周圍的土壤性質異質性較小、環境因子類似，便于分析統計。在每個樣方區內采用網格法采集表層土壤樣品10～15個，采樣深度為20 cm，并用數碼相機對采樣區進行拍照，實地記錄采樣點土壤表面狀況及景觀特征。將采集的土壤樣本在實驗室內自然晾干，研磨后過1 mm篩，剔出土壤中的石塊等雜質，按水土比5∶1配置土壤浸提液，在樣品中加入熱硝酸溶解后并采用電感耦合等離子體原子發射光譜法（ICP-AES）測定重金屬元素含量，最后得到包括鉻（Cr）、錳（Mn）、鐵（Ee）、鎳（Ni）、銅（Cu）、鋅（Zn）、砷（Sb）、鎘（Cd）、鉛（Pb）9種元素在內的42種土壤鹽分離子含量，去掉異常值后，每個樣方區采樣點重金屬元素含量取平均值，即得到所需的土壤重金屬含量數據。

圖1　采樣點示意圖Eig.1 Sketch diagram of the sampling location

2.2土壤樣本重金屬含量分析

本研究對32個土壤樣本的重金屬含量進行了簡單的統計分析，其結果見表1。

由表1可以看出:研究區土壤樣本中重金屬元素含量的變異系數普遍較小，其中Cu的變異系數最大，為9.79%；從偏度上看，除Zn以外，其他重金屬元素都屬于中高度偏態分布；另外，所有重金屬元素含量的最大值和最小值之間相差較大，說明研究區內重金屬元素的分布不均衡。

2.3高光譜遙感數據的獲取

遙感數據來源于我國環境小衛星（HJ-1 A）的HSI高光譜數據。我國環境小衛星在2009年3月30日開始正式交付使用，分為A星和B星，簡寫為HJ-1A/1B。兩顆星都配備了CCD相機，其中HJ-1A星搭載的超光譜成像儀（HSI）能夠完成對地展寬為50 km、地面像元分辨率為100 m、110～128個光譜譜段的推掃成像，具有±30°側視功能和星上定標功能，重訪周期為4 d，能夠迅速、大尺度地觀測大氣、土壤和植被等狀況，可為反演計算提供數據支持［9-10］。

表1　土壤樣本重金屬含量分析結果Table 1 Contents of the heavy metals in soil-samples

本研究高光譜遙感數據源選擇2013年10月覆蓋整個研究區的3景HJ-1A高光譜2級產品，要求所選取的影像清晰，無云或少云覆蓋，且與采集土壤樣本的時間盡量接近，用以提高土壤重金屬含量的反演精度。

2.4HSI高光譜數據的預處理

HSI高光譜數據在應用于重金屬含量反演之前，需要進行包括幾何糾正（誤差控制在半個像元內）、大氣校正、投影變化、拼接和裁剪等在內的幾何及光譜方面的預處理。本數據源屬于2級產品，在下載之前已經完成了波譜復原和系統粗校正等處理。產品在生成時選用了100的擴大因子，所以需要將像元DN值還原為標準輻射亮度，即

式中:Lλ為真實輻射亮度值［W/（m2·sr·μm）］；DN為原始圖像像元值［W/（m2·sr·μm）］。

通過對不同地物輻射亮度光譜曲線的觀察，發現光譜曲線的總體形態與趨勢相似，且受大氣中水汽的影響，不同地物在760 nm和820 nm附近均存在兩個明顯的波谷，如果直接利用未經校正的影像，在水汽的影響下將很難體現地物的真實物理特性，從而影響遙感反演的精確度。

ENVI5.1能夠充分滿足遙感影像在讀取、配準和大氣校正等方面的需求，因此本研究利用該軟件提供的ELAASH模型，并通過讀取影像頭文件獲取的模型輸入參數對HJ-1A數據進行大氣校正，最終得到115個波段的地表反射率圖像，此時的光譜曲線可以較好地體現出真實地物的光譜特征。HSI高光譜數據覆蓋450～950 nm波段，而在480 nm波段之前，即前10個波段的噪聲非常大，信息的準確性不高，因此信息量主要集中在480～950 nm波段之間［11］。

3　研究過程

3.1模型方法的選擇

高光譜遙感擁有較強的連續性波段、豐富的光譜信息和高分辨率等優點，它在特定光譜域內蘊含著近似連續的地物光譜信息，使地物的識別能力極大地提高［12-13］。然而由于高光譜數據的波段較多、信息冗余現象嚴重，使得各波段間存在著顯著的相關性，而利用傳統的最小二乘回歸方法分析高光譜數據建立的估計模型時，無法避免地會遇到多重共線性和自由度過小的問題，將造成模型方程系數矩陣的嚴重病態性，使模型參數的最小二乘估計失真。

偏最小二乘回歸（PLSR）方法是在傳統的最小二乘回歸分析方法的基礎上，引入主成分分析和方差分析的思想，為一種多變量回歸分析方法［14］?；谠摲椒?，利用高光譜數據對土壤重金屬含量進行反演計算，其主要思路是從高光譜原始數據及其變換形式中提取能夠最大解釋土壤鹽分含量變化的主成分，提取的主成分是各波段光譜數據的線性組合，各個主成分之間是正交關系，既最大程度地保存了信息量，又消減了這些主成分線性相關性的影響。

PLSR方法并不直接考慮因變量對自變量的回歸建模，而是從自變量中逐步提取多個對自變量系統X和因變量系統Y都具有最佳解釋能力的新綜合變量t1，亦稱之為主成分。首先建立yk對主成分的MLR回歸方程；然后還原為yk關于原自變量系統的PLS回歸方程，PLS關鍵是提取主成分；最后利用模型的交叉有效性Q2h判定回歸方程是否達到滿意的精度。為了判定模型的精度，引入判定系數（R2）、均方根誤差（RMSE）和相對均方根誤差［RMSE（%）］3個判定指標，R2越大、RMSE和RMSE越小，表明該模型穩定性和精度越好。均方根誤差（RMSE）亦稱標準誤差，其計算公式為

相對均方根誤差（RMSE）計算公式為

上式中:εi為土壤樣本重金屬元素含量實測值（mg/kg）；?εi為土壤樣本重金屬元素含量PLSR模型預測值（mg/kg）；ˉε為土壤樣本重金屬元素含量均值；n為樣本個數。

本文從32個采樣點中隨機抽取20個作為建模樣本，余下12個作為檢驗樣本，采用SIMCA-P軟件進行統計分析。

3.2光譜指標的提取

為了提高反演效果，將從遙感影像獲取的原始光譜反射率進行一系列的光譜變換，即以HSI影像的光譜反射率及多種數學變換形式為特征向量，提取土壤光譜信息。這些數學變換形式包括光譜去除包絡線（CR）、光譜反射率倒數的對數變換［Log（1/ R）］、光譜反射率的一階微分（RD1）和二階微分（RD2）變換。

光譜去除包絡線是一種常用的光譜處理方法，它是將反射率歸一化在0～1.0范圍之內，并將其歸一化到統一的背景值上，突出了光譜的反射和吸收特征［15-16］。將光譜反射率進行倒數的對數變換（log （1/R））是處理非線性問題的方法，它不僅趨向于增強可見光區的光譜差異，而且趨向于減少因光照條件變化引起的乘性因素的影響［17］。對反射率值進行微分變換有助于限制低頻噪聲對目標光譜的影響［18-19］。對于光譜反射率的微分計算，一般用光譜的差分作為微分的有限近似，光譜數據的一階微分計算公式近似為

式中:λi為波長（nm）；Δλ為相鄰的兩個波長λi+1和λi的間隔（nm）；ρ′（λi）為波長λi的一階微分。

3.3模型分析

為了研究各種重金屬離子在HSI影像中的光譜響應，將采樣點重金屬離子含量與影像的反射率及光譜變換形式建立偏最小二乘回歸模型，其建模結果見表2。

由表2可以看出:建模所用的光譜指標全部是光譜反射率的數學變換形式，包括微分變換、去除包絡線和反射率倒數的對數變換等；由于在前期對原始數據采用9點加權移動平均法進行平滑處理，所以選擇的主成分個數較少；對于建模效果而言，Cu的建模效果最好，判定系數（R2）高達0.965 5，Mn、Ni和Pb的建模效果較好，R2分別為0.927 3、0.910 4和0.912 7，As和Zn的R2分別為0.855 5 和0.832 7，Cr的建模效果較差，R2為0.640 5；驗證樣本中，Mn的R2最大，為0.912 8，說明該模型對Mn的預測能力最強，As、Cu和Ni的模型預測能力相當，R2分別為0.727 6、0.780 4和0.733 3，Cr和Zn的模型預測能力一般，R2分別為0.598 9 和0.622，Pb的模型預測能力較差，R2為0.515?？梢?，利用光譜反射率及其變換形式估算土壤重金屬離子含量大部分能夠獲得較大的判定系數，取得了較好的反演效果，說明HSI影像對土壤重金屬離子具有較強的探測能力。

表2　PLSR建模結果比較Table 2 Comparison of the PLSR modeling results

圖2為利用PLSR方法建立的土壤重金屬含量預測值與實測值比較的效果圖。實測值與預測值越均勻分布于1∶1直線兩側，說明基于偏最小二乘回歸方法建立的土壤重金屬含量估測模型反演效果越好。

通過對建模結果的比較，選擇相應的敏感波段，獲得了土壤重金屬元素含量的回歸方程（見表3），其中Xn表示重金屬元素含量第n個波段經過特定光譜變換形式后的數值。

圖2　土壤重金屬含量實測值與預測值的比較Eig.2 Comparation of measured and predicted values of heavy metal ions

表3　反射率不同變化形式與土壤重金屬含量回歸方程Table 3 Regression equations between different variations of the reflectivity and the soil heavy metal content

利用土壤各重金屬元素含量的回歸模型，并結合HSI遙感影像反演了研究區范圍內土壤重金屬含量分布，其結果見圖3。由圖3可以看出:龍口地區土壤的重金屬污染主要集中在西北部沿海地區，隨著污水灌溉年限的延長，污水中的重金屬離子逐漸在龍口地勢低洼的地區積累起來。

圖3　研究區土壤重金屬含量分布圖Eig.3 Distribution map of the contents of soil heavy metal ions

4　結 論

本文以HSI高光譜遙感影像為數據源，利用地理信息系統和偏最小二乘回歸（PLSR）分析方法，并結合實地采樣檢測的重金屬含量數據，分析了土壤重金屬離子在HSI數據中的光譜特征，并建立土壤重金屬含量與高光譜數據的偏最小二乘回歸模型，實現了對龍口地區土壤重金屬含量的反演研究，得到如下結論:

（1）利用HSI原始高光譜數據及其變換形式反演土壤重金屬含量，能夠獲得較好的建模效果，從而證明了環境小衛星HJ-1A影像對土壤重金屬離子具有較好的響應能力，其空間覆蓋能力可滿足將來大面積的土壤重金屬含量反演，具有一定的可行性。

（2）偏最小二乘回歸（PLSR）方法能夠有效地解決樣本中多重共線性和自由度過小等問題，在非正態分布、小樣本的情況下，仍然很穩定。

（3）從建模效果較好的重金屬出發，對龍口污水灌溉地區的土壤重金屬含量進行反演計算，初步嘗試了由點到面的研究，研究結果表明龍口污水灌溉地區的土壤重金屬污染主要集中在西北部沿海地區。

但是由于本次研究樣本數據較少，HSI高光譜數據空間分辨率較低，僅為100 m×100 m，在混合像元的影響下必然降低了反演精度，因此在進一步的研究中應該著重考慮光譜重建的問題。

［1］張淑麗.淺析污水灌溉環境污染的研究進展［J］.科技情報開發與經濟，2009，19（27）:128-129.

［2］陳建勇，張江山，鄭育毅.基于熵權的物元分析法在土壤重金屬污染評價中的應用［J］.安全與環境工程，2011（5）:57-60.

［3］龔紹琦，王鑫，沈潤平，等.濱海鹽土重金屬含量高光譜遙感研究［J］.遙感技術與應用，2010，25（2）:169-177.

［4］徐明星，吳紹華，周生路，等.重金屬含量的高光譜建模反演:考古土壤中的應用［J］.紅外與毫米波學報，2011，30（2）:109-114.

［5］劉煥軍，張柏，王宗明，等.基于反射光譜特征的土壤鹽堿化評價［J］.紅外與毫米波學報，2008，27（2）:138-142.

［6］劉華，張利權.崇明東灘鹽沼土壤重金屬含量的高光譜估算模型［J］.生態學報，2007，27（8），3427-3434.

［7］王維，沈潤平，吉曹翔.基于高光譜的土壤重金屬銅的反演研究［J］.遙感技術與應用，2011，26（3），348-354.

［8］Kooistra L，Wehrens R，Leuven R S E W，et al.Possibilities of visible-near-infrared spectroscopy for the assessment of soil contamination in river flood plains［J］.Analytica Chimica Acta，2001，446（1/2）:97-105.

［9］李傳榮，賈媛媛，胡堅，等.HJ-1光學衛星遙感應用前景分析［J］.國土資源遙感，2008，9（3）:1-9.

［10］中國資源衛星應用中心.HJ-1A、B衛星介紹［EB/OL］.http://www.cresda.com/n16/n1130/n1582/8384.html.

［11］馬馳.基于HJ-1A反演松嫩平原土壤鹽分含量［J］.干旱區研究，2014，31（2）:226-230.

［12］浦瑞亮，宮鵬.高光譜遙感及其應用［M］.北京:高等教育出版社，2000.

［13］宮鵬，黎夏，徐冰.高分辨率影像解譯理論與應用方法中的一些研究問題［J］.遙感學報，2006，10（1）:1-5.

［14］李世玲.基于PLS成分的變量篩選法［J］.信息與電子工程，2003，1（2）:31-35.

［15］何挺，王靜，林宗堅，等.土壤有機質光譜特征研究［J］.武漢大學學報，2006，31（11）:975-979.

［16］張雪紅，田慶久.基于連續統去除法的冬小麥葉片氮積累量的高光譜評價［J］.生態學雜志，2010，29（1）:181-186.

［17］Clark R N.Spectral properties of mixtures of montmorillonite and dark carbon grains:Implications for remote sensing minerals containing chemically and physically adsorbed water［J］.Journal of Geophysical Research，1983，88:10635-10644.

［18］胡軍杰，陳美蘭，吳連喜，等.鈾礦區植物的光譜特性分析［J］.東華理工大學學報（自然科學版），2011，34（4）:374-400.

［19］Huang Z，Turner B J，Stephen J D，et al.Estimating foliage nitrogen concentration from HYMAP data using continuum removal analysis［J］.Remote Sensing of Environment，2004，93: 18-29.

Quantitative Retrieval of Soil Heavy Metal Content in Longkou Wastewater Irrigation Area Based on HJ1A-HSI Images

CAI Dongquan1，WU Quanyuan1，CAO Xuejiang2，MENG Yonghui3，ZHOU Liyuan1，LIU Shuo1，ZHANG Longlong1
（1.College of Population，Resources and Environment，Shandong Normal University，Jinan 250014，China；2.Shandong Provincial Coal Geological Planning，Exploration and Research Institute，Taian 271000，China；3.Shandong Monitoring Center of Geological Environment，Jinan 250014，China）

Wastewater irrigation can alleviate water problems in agriculture to some extent，but it causes the accumulation of heavy metals.This paper adopts the HSI hyperpectral remote sensing images as the data source using geographic information system and PLSR analysis method.Combining with concentration data of heavy metals through field sampling and monitoring，this paper analyzes spectral characteristics of heavy metal ions in HSI data，and establishes the equation between concentration data of heavy metals and hyperspectral data to achieve the quantitative retrieval of soil heavy metals contents in Longkou wastewater irrigation area.After spectral transformations of the original spectral data，Cu element is optimal in terms of modeling performance，with the determination coefficient（R2）up to 0.965 5；Mn，Ni and Pb are suboptimal with the R2values of 0.927 3，0.910 4 and 0.912 7，respectively.In validation samples，the R2of Mn is 0.912 8，indicating the best retrieval performance.And then it is followed by As，Cu and Ni，with the R2values of 0.727 6，0.780 4 and 0.733 3，respectively.Therefore，it is feasibleand reliable to retrieve soil heavy metal contents in Longkou wastewater irrigation area by using the data of HJ1AHSI and method of PSLR.

soil heavy metals；HSI；partial least squares regression（PLSR）；wastewater irrigation

X53

A

10.13578/j.cnki.issn.1671-1556.2015.05.006

1671-1556（2015）05-0033-07

2015-03-30

2015-08-11

國家自然科學基金項目（41371395）；國家科技支撐計劃項目（2012BAB11B01）；龍口礦區及周邊海岸帶遙感監測研究項目（魯勘字［2012］110號）；海（咸）水入侵典型區數值模型與預警預報系統開發建設項目（SDZS-2015-DHJ01）

蔡東全（1989—），男，碩士研究生，主要研究方向為遙感地學應用。E-mail:cdq8906@126.com

吳泉源（1959—），男，博士，教授，主要從事遙感地學應用方面的研究。E-mail:wqy642052582@163.com