基于高光譜的射陽河口懸浮泥沙濃度定量反演研究

潘 潔,張 鷹,譚子輝

(1.南京林業(yè)大學 森林資源與環(huán)境學院,江蘇 南京 210037;2.南京師范大學 地理科學學院,江蘇 南京210042;3.山東省臨沂市蔬菜辦公室,山東 臨沂276000)

水體懸浮物含量是重要的水質(zhì)參數(shù),國內(nèi)外諸多研究建立了多種基于遙感技術(shù)反演水體懸浮物含量的理論和經(jīng)驗?zāi)Ｐ?常用的主要有對數(shù)模式、Gordon模式和指數(shù)模式[1]。此外,我國學者還提出了一些改進算法,李炎等[2]提出以NOAAAVHRRCH1和CH2波段反射率差的最大值條帶作為大氣校正參考點,提出了以海面-遙感器的光譜反射率斜率傳遞現(xiàn)象為基礎(chǔ)的斜率算法;汪小欽等[3]利用兩個時相的 TM 影像,用線性光譜混合分析法提取福建閩江口的懸浮物濃度;李四海等[4]在 2000年根據(jù)多時相NOAA/AVHRR衛(wèi)星遙感數(shù)據(jù)和準同步實測表層含沙量資料,利用斜率法、灰度法和泥沙指數(shù)法三種模式分別建立泥沙遙感定量模式,并對其效果和適用性進行了比較。近年來,高光譜影像也開始應(yīng)用到懸浮泥沙監(jiān)測領(lǐng)域中來,Hyperion影像兼?zhèn)溆胁ǘ钨Y源豐富和空間分辨率高的特點,在水體懸浮泥沙濃度遙感提取方面具有很大的應(yīng)用潛力[5]。

對河口而言,懸浮泥沙濃度不僅是河口地區(qū)港口與航道工程十分關(guān)心的問題,同時細顆粒泥沙又是各種營養(yǎng)鹽和污染物的重要載體,對河口水質(zhì)環(huán)境研究也相當重要。河口區(qū)的懸浮泥沙在徑流潮流相互作用、鹽淡水混合以及風浪等多種動力因素的作用下,其運動規(guī)律極其復(fù)雜,時空分布變化快。本研究利用射陽河口實測光譜及模擬Hyperion光譜數(shù)據(jù)與表層泥沙濃度進行相關(guān)性分析,構(gòu)建定量模型,從而實現(xiàn)了射陽河口水體懸浮泥沙濃度的定量反演。

1 數(shù)據(jù)獲取與處理

1.1 研究區(qū)概況

射陽河口位于江蘇沿海侵蝕性岸段與淤長性岸段之間,其河口北側(cè)屬于侵蝕性岸段,南側(cè)屬于淤長性岸段。由于水動力條件復(fù)雜,懸浮泥沙質(zhì)量濃度最大可達到11.08 kg/m3,最小僅為0.035 kg/m3,攔門沙和航道的位置也經(jīng)常發(fā)生變化[6](圖1)。

射陽河口作為射陽河的入海口,射陽河排污量大,由 1996年全省海洋調(diào)查結(jié)果可知,與全省其他入海河口相比,射陽河口水質(zhì)最差,主要以油類、無機氮磷及重金屬污染為主[7]。

1.2 水體光譜測量數(shù)據(jù)及懸沙質(zhì)量濃度獲取

水體光譜測量主要是為獲取離水輻射率、歸一化離水輻射率和遙感反射率等遙感參數(shù),進而通過這些參數(shù)反演得到水體信息。對于二類水體,表面以上測量法是目前唯一有效的方法[8]。

本研究對于射陽河口的水體光譜測量共進行了兩次。試驗1測量時間是2007年7月1日,如圖1所示,本次試驗設(shè)置5個測量點(S1,S2,S3,S4,S5),通過GPS定位,S1,S2,S3三個測點以及S3,S4,S5三個測點分別位于一條直線。其中,S1,S2,S3三個測點從9:00至15:00每隔1h連續(xù)測量水體光譜值。試驗2光譜測量時間是2008年5月29～31日,在射陽河口附近海域進行了三個航次的測量和采樣,共在60個測點測量了水體的光譜(圖1)。每天測量的時間為9:30～14:30。為了使所建立的模型在射陽河口及其附近海域更具有代表性,第三個航次的取樣點選擇在離河口較遠的海域,以確保樣本點不僅僅代表河口,還能代表河口附近海域。

圖1 射陽河口研究區(qū)2007年、2008年采樣點位置Fig.1 Sampling sites in 2007 and 2008

兩次試驗光譜測量所用光譜儀是ASD公司生產(chǎn)的FieldSpec地物光譜儀,該儀器測定的光譜范圍為282～1 090 nm,標準板是經(jīng)嚴格定標的反射率為0.3的灰板。在每個樣點上按上述觀測幾何條件分別測取水體、天空光和標準板的DN值各10條,并同時采取水樣。同時用 GPS定位并記錄當時的風速,以便于對影像進行大氣校正。

在測量水體光譜的同時,每個樣點(包括試驗 1中的 S4,S5樣點)同時采取水樣,采用烘干法獲取各樣品懸沙質(zhì)量濃度。兩次試驗所獲取的懸沙質(zhì)量濃度的量值范圍分別為 0.035～2.15 kg/m3(試驗 1)和0.04～3.25 kg/m3(試驗 2)。

1.3 Hyperion影像數(shù)據(jù)獲取

本研究選用的Hyperion影像是由美國地質(zhì)調(diào)查局(USGS)處理后生成的L1R產(chǎn)品,其VNIR波段和SWIR波段間的空間錯位已經(jīng)經(jīng)過糾正[9],成像時間為北京時間2007年7月11日。Hyperion數(shù)據(jù)采用HDF(Hierarchical Data Format)數(shù)據(jù)集的形式存儲,波段存儲格式為 BIL格式,由于后面的預(yù)處理程序需要在ENVI軟件環(huán)境下運行,因此在對數(shù)據(jù)進行處理前需要將其數(shù)據(jù)存儲格式轉(zhuǎn)化為 ENVI軟件能夠直接讀取的標準格式。

對 Hyperion影像的預(yù)處理包括:波段剔除、壞線去除、條紋去除、大氣校正和幾何校正等。考慮到輻射定標、波段噪聲與水體信息提取,Hyperion影像的 242個波段中,我們最后保留了 8～57波段(430～930 nm)共50個波段用于研究。利用像元灰度斜率閾值法可以實現(xiàn)影像各波段的壞線判別。在進行壞線修復(fù)以后,影像上還存在另外一種像元值異常情況,即垂直條紋。針對條紋的垂直分布特征,采樣全局歸一化法GNM(Global Normalization Method),通過像元的列平均值、標準差與波段平均值、標準差之間的差異對像元進行分波段線性化修正,可以消除垂直條紋的影響,利用ENVI IDL編制程序,實現(xiàn)對Hyperion影像的壞線修復(fù)和垂直條紋去除。采用ENVI軟件自帶的基于MODTRAN模型的大氣校正模塊FLAASH軟件來實現(xiàn)Hyperion影像的大氣校正。影像的幾何精校正采用經(jīng)過幾何精校正的的影像底圖進行影像到影像的校正[10]。

2 結(jié)果與分析

2.1 射陽河口含沙水體光譜特征分析

通過對射陽河口不同懸沙質(zhì)量濃度水體的光譜特征分析(圖2)可以看出,隨著懸沙質(zhì)量濃度增大,水體光譜反射率有升高的趨勢。在350～500 nm之間,含沙水體反射率相對較低;在560～720 nm之間有一個反射峰,當泥沙濃度較小時,其峰值主要在560～610 nm 之間,且峰值反射率較低,當泥沙濃度增大時,其峰移在690～720 nm之間,且峰值反射率較高。此外,在790～820 nm之間,還有一個反射峰,當泥沙濃度較低時,該峰值明顯低于第一個反射峰,隨著泥沙濃度的上升,該峰值迅速上升,并接近第一個反射峰的值;850 nm之后反射率趨于下降。

圖2 不同懸沙質(zhì)量濃度水體的反射光譜曲線Fig.2 Spectra of water samples with different suspended sediment contents

圖3 光譜反射率與懸浮泥沙含量相關(guān)性分析Fig.3 Correlation between the ASD-based reflectance and suspended sediment content

對試驗2中60個測點的懸浮泥沙濃度與對應(yīng)的實測光譜反射率值的進行相關(guān)分析,結(jié)果如圖3所示,由于在小于350 nm和大于910 nm的波段受噪聲干擾較大,所以此分析只取350～910 nm之間的數(shù)據(jù)進行。結(jié)果表明,在波長350～500 nm,懸沙質(zhì)量濃度與光譜反射率呈現(xiàn)負相關(guān),且相關(guān)系數(shù)小于 0.2;波長大于500 nm開始,懸沙質(zhì)量濃度與光譜反射率的正相關(guān)性逐漸增大,最大相關(guān)系性出現(xiàn)在波長 898～904 nm間,最大相關(guān)系數(shù)值為0.8858。

2.2 基于Hyperion高光譜數(shù)據(jù)的射陽河口懸浮泥沙濃度定量模型構(gòu)建

模型構(gòu)建前,有必要對實測的光譜數(shù)據(jù)按光譜響應(yīng)函數(shù)進行光譜重采樣,以模擬傳感器的光譜響應(yīng)特征。由于未能從NASA獲得Hyperion的光譜響應(yīng)函數(shù),而 Hyperion是高光譜影像,各波段半高寬(FWHM)很窄,可以用高斯函數(shù)模擬它的光譜響應(yīng)函數(shù)[11]。將模擬得到的Hyperion各波段反射率及波段組合因子分別與懸沙質(zhì)量濃度作相關(guān)性分析,通過相關(guān)性分析尋找敏感波段或波段組合來構(gòu)建射陽河口懸沙質(zhì)量濃度的定量反演模型。已有研究表明,一階或二階微分因子雖能體現(xiàn)高光譜數(shù)據(jù)的優(yōu)勢,但這一因子對傳感器的靈敏度和信噪比的要求很高,對大氣校正等影像預(yù)處理的要求也十分苛刻,如果通過影像不能獲得與水面實測光譜完全一致的光譜曲線,則用微分因子構(gòu)建的模型很難在影像反演中獲得令人滿意的結(jié)果[11]。因此,本研究不再考慮用微分因子作為敏感波段構(gòu)建反演模型。

由圖4所示,Hyperion光譜反射率與懸沙質(zhì)量濃度的相關(guān)性和實測光譜與懸沙質(zhì)量濃度的相關(guān)性變化趨勢基本相同,當波長大于 500 nm后,Hyperion光譜反射率與懸沙質(zhì)量濃度呈現(xiàn)正相關(guān),且最大相關(guān)系數(shù)位于中心波長 896 nm(最大相關(guān)系數(shù)值達到0.8916),這與實測光譜預(yù)測懸沙質(zhì)量濃度的敏感波段吻合。

前人的研究表明,單個波段的反射率難以全面地反映出不同泥沙濃度的光譜信息[4]。因此,對Hyperion的各波段進行了9種波段組合(表1)。通過9種波段組合分別與60個測點的懸浮泥沙濃度進行相關(guān)性分析,得出相關(guān)系數(shù)均不超過 0.8,并沒有單波段的最大相關(guān)系數(shù)高。由此可見,對泥沙質(zhì)量濃度為0.035～3.25 kg/m3范圍的河口水體,896 nm雖然僅代表了單波段信息,但其反射率對水體懸沙質(zhì)量濃度變化具有更靈敏的反應(yīng),更能反映含沙水體的光譜特征。

圖4 Hyperion光譜反射率與懸浮泥沙含量相關(guān)性分析Fig.4 Correlation between Hyperion reflectance and the suspended sediment content

圖5 Hyperion 中心波長 896nm波段光譜反射率與懸浮泥沙含量回歸分析Fig.5 Regression of the suspended sediment contents based on the Hyperion reflectance at 896 nm

表1 單波段因子和波段組合因子Tab.1 Factors of single band and bands combination

同時,通過對水體氮磷濃度、重金屬濃度與Hyperion光譜反射率的相關(guān)性分析,896 nm波段都表現(xiàn)出顯著的相關(guān)特征[12]。因此,以Hyperion 896 nm波段反射率與懸浮泥沙含量的相關(guān)性構(gòu)建射陽河口懸沙質(zhì)量濃度的定量模型(圖5),如式(1)所示:

其中,CSS指射陽河口懸沙質(zhì)量濃度(kg/m3)。模型R2值達到 0.799,顯著高于各種波段組合構(gòu)建的射陽河口懸沙質(zhì)量濃度反演模型。

2.3 模型檢驗

以2007年進行的試驗1中S1,S2和S3三個監(jiān)測位點的實測懸沙質(zhì)量濃度及光譜測量數(shù)據(jù)對上述模型進行檢驗。3個位點光譜測量數(shù)據(jù)自中午12:00至下午15:00,每隔1 h測量光譜,并在對應(yīng)時間取水樣,利用烘干法獲取懸浮濃度,共獲取檢驗數(shù)據(jù)12組。實測光譜同樣進行重采樣以模擬Hyperion影像各波段的反射率。模型檢驗采用常用的相對根均方差法(relative root mean square error,以E表示)對模擬值與觀測值之間的符合度進行統(tǒng)計分析,E值的計算見方程(2)。

其中,Pi和Oi分別為預(yù)測值和觀測值,為觀測值的平均值[13]。如表2所示,以模擬Hyperion光譜數(shù)據(jù)模擬射陽河口懸沙質(zhì)量濃度模型檢驗的E為38.7757%,表明模型具有較好的預(yù)測精度。

表2 射陽河口表層懸沙質(zhì)量濃度預(yù)測誤差Tab.2 Predicted errors of suspended sediment contents for Sheyang estuary

繪制觀測值與預(yù)測值之間的1 :1關(guān)系圖,以直觀的展示模擬值與觀察值的擬合度和可靠性[14],結(jié)果顯示模擬值與觀測值間具有較好的符合度(圖6)。

圖6 射陽河口懸浮泥沙質(zhì)量濃度預(yù)測值與實測值比較Fig.6 Modeled and observed values of the suspended sediment contents for Sheyang estuary

2.4 模型反演

將上述射陽河口懸沙質(zhì)量濃度回歸模型應(yīng)用于懸浮物濃度反演制圖,結(jié)果如圖7所示。

圖7 2007年7月11日Hyperion影像射陽河口懸沙質(zhì)量濃度反演結(jié)果Fig.7 Result of model simulation on suspended sediment content by Hyperion image

由圖7所示,射陽河口懸浮泥沙濃度由近岸向遠海呈現(xiàn)遞減的規(guī)律,在河口南側(cè)的懸沙質(zhì)量濃度高于北側(cè),主要是由于河口南側(cè)屬于淤長性海岸。

3 總結(jié)

本文利用射陽河口實測光譜及模擬Hyperion光譜數(shù)據(jù)與表層泥沙濃度進行相關(guān)性分析,優(yōu)選了896 nm波段作為研究水體懸浮物濃度的敏感波段并據(jù)此建立了定量反演模型,實現(xiàn)了射陽河口水體懸浮泥沙濃度的遙感制圖。盡管已有的研究表明單個波段的反射率難以全面地反映出不同泥沙濃度的光譜信息,但針對射陽河口水域,其懸沙類型為淤泥質(zhì)粉沙,顆粒較為細膩均一,基本沒有明顯的分級。懸浮泥沙水體對實測光譜反射率的敏感波段為898～904 nm,對光譜分辨率為10 nm的Hyperion影像敏感波段正好位于896 nm,以此構(gòu)建的指數(shù)模型,相關(guān)指數(shù)達到 0.89。模型的檢驗結(jié)果表明,相對RMSE值為38.78%,表明單波段896 nm對預(yù)測射陽河口的懸沙質(zhì)量濃度具有良好的精度。

由于模型構(gòu)建與檢驗的實測懸沙質(zhì)量濃度數(shù)據(jù)均來源于上午落潮后的4～5 h后取樣,該時間段與Hyperion影像的成像時間相匹配,此時水流場較穩(wěn)定,河口水體中的泥沙運移與擴散規(guī)律與影像解譯結(jié)果相似,即:近岸懸沙質(zhì)量濃度高,遠海濃度低;河口南側(cè)由于屬于淤長性海岸,懸沙質(zhì)量濃度比北側(cè)高。

[1]Richard L M,Brent A M.Using MODIS terra 250m imagery to map concentrations of total suspended mattering coastal water[J].Remote Sensing Environment,2004,93:259-266.

[2]李炎,李京.基于海面-遙感器光譜反射率斜率傳遞現(xiàn)象的懸浮泥沙遙感算法[J].科學通報,1994,44(17):1892-1897.

[3]汪小欽,王欽敏,鄔群勇,等.遙感在懸浮物質(zhì)濃度提取中的應(yīng)用——以福建閩江口為例[J].遙感學報,2003,7(1):54-57.

[4]李四海,惲才興.河口表層懸浮泥沙氣象衛(wèi)星遙感定量模式研究[J].遙感學報,2001,5(2):154-160.

[5]Bowers D,Gaffney S,White M.Turbidity in the southern Irish Sea[J].Continental Shelf Research,1998,18(5):487-500.

[6]張忍順,陸麗云,王艷紅.江蘇海岸侵蝕過程及其趨勢[J].地理研究,2002,21(4):469-478.

[7]江蘇省海洋污染基線調(diào)查隊.江蘇省海洋污染基線調(diào)查報告[M].南京:河海大學出版社,2001:10-56.

[8]唐軍武,田國良,汪小勇,等.水體光譜測量與分析Ⅰ:水面以上測量法[J].遙感學報,2004,8(1):37-44.

[9]譚炳香,李增元,陳爾學,等.EO-1 Hyperion高光譜數(shù)據(jù)的預(yù)處理[J],遙感信息,2005,6:36-41.

[10]張東,張鷹,李歡.海岸帶星載高光譜遙感影像預(yù)處理方法[J].海洋科學進展,2009,27(1):92-97.

[11]趙祥,梁順林,劉素紅,等.高光譜遙感數(shù)據(jù)的改正暗目標大氣校正方法研究[J].中國科學 D輯,2007,37(12):1653-1659.

[12]Pan Jie,Zhang Ying,Xu Yong.Assessment of water quality in Sheyang Estuary (China) using hyperspectral data[C]//Tan Zhenghua.2010 3rdInternational Congress on Image and Signal Processing Symposium.Yantai:Yantai University,2009:1-10.

[13]Michele R,Nicola L,Zina F.Evaluation and application of the OILCROP –SUN model for sunflower in southern Italy[J].Agricultural Systems,2003,78:17-30.

[14]Snyder R L,Spano D,Cesaraccio C,et al.Determining degree-day thresholds from field observations[J].International Journal of Biometeorology,1999,42(4):177-182.