膠州灣海域Landsat8/OLI數(shù)據(jù)處理中多種大氣校正方法的評(píng)價(jià)

劉曉燕， 申 辰， 崔文璽， 楊 倩1, ， 禹定峰1, *，

高 皜1, 3， 楊 雷1, 3， 周 燕1, 3， 趙新興3

1. 齊魯工業(yè)大學(xué)(山東省科學(xué)院)海洋儀器儀表研究所， 山東 青島 266100

2. 中國(guó)海洋大學(xué)信息科學(xué)與工程學(xué)部海洋技術(shù)系， 山東 青島 266010

3. 齊魯工業(yè)大學(xué)(山東省科學(xué)院)海洋技術(shù)科學(xué)學(xué)院， 山東 濟(jì)南 250300

引 言

離水輻亮度與恰好水面之上的向下輻照度之比被定義為水體光譜的遙感反射比(Rrs)， 高精度的海洋遙感反演研究需要精確的Rrs數(shù)據(jù)[1]。 海洋水色衛(wèi)星傳感器接收到的總信號(hào)中水體貢獻(xiàn)的離水輻亮度不足10%， 我們需要對(duì)接收的信號(hào)進(jìn)行大氣校正， 這也是海洋光學(xué)遙感的關(guān)鍵技術(shù)之一。 精確的大氣校正是得到高精度的Rrs數(shù)據(jù)、 實(shí)現(xiàn)高精度水色信息反演的前提。 目前近岸Ⅱ類渾濁水體的大氣校正算法有亮像元法[2]、 暗像元法[3]、 神經(jīng)網(wǎng)絡(luò)法[4]、 光譜匹配法[5]、 MUMM法[6]、 光譜優(yōu)化法[7]等。 以往通常通過與現(xiàn)場(chǎng)數(shù)據(jù)集的誤差分析來評(píng)價(jià)衛(wèi)星海洋光學(xué)數(shù)據(jù)(包含Rrs(λ)數(shù)據(jù))質(zhì)量， 誤差分析結(jié)果可以表達(dá)出海洋光學(xué)觀測(cè)數(shù)據(jù)質(zhì)量的信息[8-9]。 然而， 研究表明單獨(dú)的波段對(duì)比散點(diǎn)圖的實(shí)際評(píng)價(jià)意義并不能代表完整的Rrs(λ)光譜的質(zhì)量情況， Wei等[10]在基于大量高質(zhì)量的現(xiàn)場(chǎng)高光譜遙感反射比數(shù)據(jù)的基礎(chǔ)上建立了一個(gè)可以客觀評(píng)估Rrs(λ)光譜質(zhì)量的質(zhì)量保證系統(tǒng)(QA)， 并使用NOMAD遙感反射比數(shù)據(jù)集、 沿海和大洋的衛(wèi)星遙感反射比數(shù)據(jù)對(duì)該系統(tǒng)進(jìn)行進(jìn)一步測(cè)試， 結(jié)果顯示有問題的或可能錯(cuò)誤的Rrs(λ)光譜可以通過該評(píng)價(jià)系統(tǒng)很好地識(shí)別出來。

在過去的二十多年里， 世界各國(guó)發(fā)射了各種各樣的太空海洋水色傳感器， 以進(jìn)行每天的全球海洋觀測(cè)。 主要觀測(cè)平臺(tái)是極地軌道衛(wèi)星， 例如海洋觀測(cè)寬視場(chǎng)傳感器(SeaWiFS)、 中分辨率成像光譜儀(MODIS)、 中分辨率成像光譜儀(MERIS)、 可見紅外成像輻射計(jì)套件(VIIRS)和海洋陸地顏色儀(OLCI)等。 此外還有具有重大歷史意義的地球靜止海洋水色成像儀(GOCI)。 雖然這些傳感器數(shù)據(jù)的分辨率(SeaWiFS: 1 km; MODIS: 1 km; MERIS： 300 m; VIIRS: 400 m; OLCI: 300 m; GOCI: 500 m)可以滿足寬闊海域的遙感應(yīng)用需求， 但是對(duì)于光學(xué)性質(zhì)復(fù)雜的陸海環(huán)灣狹窄海域則需要更高分辨率的衛(wèi)星傳感器數(shù)據(jù)[11]。 美國(guó)NASA的陸地衛(wèi)星計(jì)劃Landsat最初是為了觀測(cè)陸地而研發(fā)， 但是其高空間分辨率(30 m)的優(yōu)勢(shì)在海洋遙感監(jiān)測(cè)中表現(xiàn)突出。 眾多科技人員應(yīng)用Landsat衛(wèi)星數(shù)據(jù)對(duì)地球上的河流、 湖泊、 內(nèi)陸環(huán)灣的生物光學(xué)特性等進(jìn)行了研究[12-13]。

膠州灣海域是位于山東半島的一個(gè)半封閉淺水海灣， 面積近500 km2， 隸屬黃海西部， 是北溫帶海灣生態(tài)系統(tǒng)的重要代表。 灣內(nèi)營(yíng)養(yǎng)鹽高， 是魚、 蝦、 貝等的生長(zhǎng)繁殖的優(yōu)良海域， 隸屬山東省海洋牧場(chǎng)養(yǎng)殖區(qū)的規(guī)劃范圍， 其水體光學(xué)性質(zhì)復(fù)雜。 對(duì)該海域?qū)嵭兴|(zhì)參數(shù)遙感監(jiān)測(cè)， 為海洋牧場(chǎng)養(yǎng)殖單位提供水質(zhì)環(huán)境監(jiān)測(cè)參數(shù)及其預(yù)警有重要意義。 Landsat8/OLI的高空間分辨率使得它成為該研究的不可或缺的衛(wèi)星數(shù)據(jù)。 由Landsat8/OLI原始數(shù)據(jù)得到水體的遙感反射比數(shù)據(jù)， 需要進(jìn)行大氣校正， 目前可選擇的大氣校正算法頗多， 本文基于QA Score光譜質(zhì)量評(píng)價(jià)體系對(duì)Landsat8/OLI數(shù)據(jù)處理中的五種大氣校正算法(Seadas-Default， Acolite DSF， Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR)在膠州灣海域的適應(yīng)性進(jìn)行了評(píng)價(jià)分析。 膠州灣海域的地理位置見圖1。

1 數(shù)據(jù)介紹

Landsat是美國(guó)NASA的陸地衛(wèi)星計(jì)劃， 1972年發(fā)射Landsat1， 至今已發(fā)射8顆。 目前， 在役服務(wù)的有Landsat8, 于2013.02.12發(fā)射成功。 Landsat8上攜帶了2個(gè)傳感器： 運(yùn)行陸地成像儀(operational land imager， OLI)和熱紅外傳感器(thermal infrared sensor， TIRS)。 本工作應(yīng)用的是Landsat8/OLI數(shù)據(jù)， 波段設(shè)置見表1。 OLI包括9個(gè)波段， 空間分辨率為30 m， 其中包括一個(gè)15 m的全色波段， 成像寬幅為185 km×185 km。 數(shù)據(jù)下載地址： http://www.gscloud.cn/sources。 其在膠州灣海域的條帶號(hào)為120， 行列號(hào)為35。

圖1 膠州灣海域地理位置圖

表1 Landsat8/OLI波段設(shè)置

選用了2013年—2017年膠州灣海域Landsat8/OLI的部分晴空數(shù)據(jù)(表2)， 基于QA Score光譜質(zhì)量評(píng)價(jià)體系對(duì)五種大氣校正算法在該海域Landsat8/OLI的大氣校正處理結(jié)果進(jìn)行了評(píng)價(jià)分析。

表2 本文選用的Landsat8/OLI數(shù)據(jù)日期

2 算 法

2.1 大氣校正算法

2.1.1 Acolite大氣校正算法

Acolite(Atmospheric Correction for OLI lite)是由比利時(shí)皇家自然科學(xué)研究所RBINS(Royal Belgian Institute of Natural Sciences)開發(fā)的， 可用于Landsat(5/7/8)和Sentinel-2(A/B)衛(wèi)星數(shù)據(jù)處理的專業(yè)軟件。 Acolite既可以對(duì)衛(wèi)星數(shù)據(jù)進(jìn)行大氣校正， 又可以輸出從水體反射率計(jì)算得出的一系列參數(shù)。 它最初應(yīng)用IDL開發(fā)實(shí)現(xiàn)(2014—2017)， 從2018年起已轉(zhuǎn)換為Python(2018年—至今)。 本工作使用其更新于2021年1月14日的最新版本對(duì)Landsat8/OLI數(shù)據(jù)進(jìn)行輻射校準(zhǔn)和大氣校正， 從原始數(shù)據(jù)獲得水體反射率值。 Acolite目前包含兩種大氣校正算法， 默認(rèn)的“暗光譜擬合”(dark spectrum fitting， DSF)算法[14-15]， 和“指數(shù)外推”(exponential extrapolation， EXP)算法[16]。

DSF大氣校正算法不假設(shè)哪個(gè)波段的水體反射率可以忽略不計(jì)， 甚至可以使用非水目標(biāo)來估算大氣路徑反射率(ρpath)。 對(duì)于每個(gè)波段， 根據(jù)擬合的直方圖中前一千個(gè)像素的偏移量估計(jì)最暗的目標(biāo)， 創(chuàng)建“暗光譜”(ρdark)。 然后根據(jù)觀測(cè)到的暗光譜ρdark和估算得到的ρpath之間的均方根誤差的最小值去選擇大陸或海洋氣溶膠模型。 并通過從SWIR波段估算太陽(yáng)閃爍反射率來進(jìn)行閃爍校正， 以獲得無(wú)閃爍的表面反射率[14]。

EXP大氣校正算法利用瑞利校正的兩個(gè)波段的反射率來估計(jì)氣溶膠的反射率(ρa(bǔ)m)。 默認(rèn)情況下使用兩個(gè)SWIR波段(此方法記為Acolite SWIR)，ρa(bǔ)m使用2.2 μm波段估計(jì)， 并使用指數(shù)函數(shù)和1.6與2.2 μm的瑞利校正反射系數(shù)的比值外推到可見和近紅外波段。 除此之外EXP大氣校正算法還包括Red/NIR(此方法記為Acolite Red/NIR), 和NIR/SWIR(此方法記為Acolite NIR/SWIR)兩種波段選項(xiàng)以供選擇。

2.1.2 Seadas大氣校正算法

NASA Ocean Color 目前為全球用戶免費(fèi)提供Seadas軟件用于海色遙感衛(wèi)星數(shù)據(jù)的處理與分析。 Seadas軟件在遙感數(shù)據(jù)的大氣校正方面提供了多種大氣校正算法供用戶進(jìn)行選擇。 在此選用了Seadas7.5的默認(rèn)大氣校正算法(即NASA標(biāo)準(zhǔn)大氣校正算法， 記為Seadas-Default)對(duì)Landsat8/OLI數(shù)據(jù)進(jìn)行大氣校正處理， 并得到校正后的遙感反射比Rrs數(shù)據(jù)。

Seadas的默認(rèn)大氣校正算法最初是由Gordon和Wang開發(fā)， Stumpf等人將其擴(kuò)展應(yīng)用到二類水體[17]。 隨后Ahmad等[18]對(duì)其進(jìn)行了修正。 該算法是一個(gè)迭代過程， 迭代過程中考慮近紅外波段水體反射率(ρw)非零的情況， 并假設(shè)近紅外波段水體光譜由懸浮物決定[式(1)]。 同時(shí)使用了80個(gè)由AERONET觀測(cè)資料和矢量輻射傳輸程序建立的海洋-大氣系統(tǒng)的氣溶膠模型。 NASA采用該算法為Seadas軟件的默認(rèn)大氣校正算法對(duì)海色遙感數(shù)據(jù)進(jìn)行處理得到L2產(chǎn)品提供給用戶。

ρw(λ)=f(λ, SPM)

(1)

2.2 光譜質(zhì)量評(píng)價(jià)算法

QA光譜質(zhì)量評(píng)價(jià)體系主要取決于光學(xué)水體類型的聚類分析法， 其建立及應(yīng)用的步驟如下：

(1)將現(xiàn)場(chǎng)測(cè)量的大量的Rrs譜進(jìn)行歸一化處理， 得到nRrs。

(2)

式(2)中，N為從1到9的總波段數(shù)，λi對(duì)應(yīng)波長(zhǎng)412， 443， 488， 510， 531， 547， 555， 667和678 nm。nRrs值在0～1之間， 且保留了原始遙感反射比Rrs譜的譜形， 波段比與原始Rrs譜的比值相同。 Wei[10]等使用了全球較多處海域(包含渤海、 黃海海域)的958組現(xiàn)場(chǎng)高光譜測(cè)量數(shù)據(jù)對(duì)其進(jìn)行水體聚類分析。 確定了最優(yōu)聚類數(shù)量為23。

(2)得到23種光學(xué)水體九波段的nRrs的譜。 通過nRrs譜得到每種水體每個(gè)波段的nRrs譜的上邊界、 下邊界值， 以及其平均值。 這三組數(shù)值是QA評(píng)價(jià)體系的關(guān)鍵[10]。

(3)

3 結(jié)果與討論

Landsat8/OLI晴空數(shù)據(jù)經(jīng)五種大氣校正算法處理后得到同時(shí)空情況下的Rrs光譜的QA得分的概率統(tǒng)計(jì)曲線圖(圖2)。 橫坐標(biāo)表示經(jīng)過QA評(píng)價(jià)體系計(jì)算的Landsat8/OLIRrs光譜的得分， 縱坐標(biāo)表示每個(gè)得分出現(xiàn)的概率。

圖2 五種大氣校正算法的QA得分統(tǒng)計(jì)概率分布圖

統(tǒng)計(jì)不同大氣校正算法得到Rrs數(shù)據(jù)的QA得分概率詳細(xì)情況(表3)。 可以看出Seadas Default處理得到的Rrs得分值為1的概率(83.95%)要遠(yuǎn)大于Acolite DSF(49.66%)， Acolite SWIR(4.13%)， Acolite Red/NIR(7.25%)， Acolite NIR/SWIR(1.38%)四種大氣校正算法。 Acolite DSF算法優(yōu)于Acolite EXP的三種大氣校正算法(Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR)。 Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR在膠州灣海域?qū)andsat8/OLI數(shù)據(jù)的大氣校正效果較差， 不建議選用。

圖3—圖7顯示了膠州灣海域2015.03.10(春季)， 2016.09.20(近似代表夏季)， 2014.11.18(秋季)， 2016.02.09(冬季)五種大氣校正算法所得的Rrs(λ)光譜QA得分的空間分布結(jié)果。 可以看出春、 夏、 秋、 冬四個(gè)季節(jié)， 基于Seadas Default算法得到的Rrs(λ)光譜QA得分的評(píng)價(jià)結(jié)果都要優(yōu)于其他四種大氣校正算法。 對(duì)于Acolite大氣校正算法而言， Acolite DSF算法相比Acolite EXP(Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR)算法具有明顯的優(yōu)勢(shì)。 Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR三種大氣校正算法所得的Rrs(λ)光譜QA得分在膠州灣海域的空間分布上雖然存在一定的差異性， 但整體評(píng)價(jià)而言， 其大氣校正效果都較差。

表3 不同大氣校正算法Rrs數(shù)據(jù)的QA得分概率統(tǒng)計(jì)

鑒于Seadas Default算法和Acolite DSF算法處理Landsat8/OLI數(shù)據(jù)所得Rrs光譜的QA得分為1的概率大于其他三種算法， 我們對(duì)Seadas Default和Acolite DSF兩種大氣校正算法所得的同時(shí)空條件下(19 616組數(shù)據(jù))的443， 483， 561和655 nm波段Rrs數(shù)據(jù)進(jìn)行了比較(圖8)。 結(jié)果顯示Seadas Default算法和Acolite DSF算法所得的443， 483， 561和655 nm波段Rrs數(shù)據(jù)之間存在較大的差異性。

圖3 膠州灣海域Acolite DSF 大氣校正算法處理Landasat8/OLI數(shù)據(jù)的Rrs(λ) QA得分空間分布圖

圖4 膠州灣海域Acolite SWIR 大氣校正算法處理Landasat8/OLI數(shù)據(jù)的Rrs(λ) QA得分空間分布圖

圖5 膠州灣海域Acolite Red/NIR大氣校正算法處理Landasat8/OLI數(shù)據(jù)的Rrs(λ) QA得分空間分布圖

圖6 膠州灣海域Acolite NIR/SWIR大氣校正算法處理Landasat8/OLI數(shù)據(jù)的Rrs(λ)QA得分空間分布圖

圖7 膠州灣海域Seadas Default大氣校正算法處理Landasat8/OLI數(shù)據(jù)的Rrs(λ)QA得分空間分布圖

為了對(duì)兩種大氣校正算法所得各個(gè)波段遙感反射比數(shù)據(jù)進(jìn)行定量分析。 鑒于其他研究結(jié)果顯示MODIS衛(wèi)星數(shù)據(jù)在該海域的普適性[20]， 我們搜集了膠州灣海域同時(shí)空條件下的MODIS衛(wèi)星Rrs數(shù)據(jù)， 將其與Seadas Default和Acolite DSF兩種大氣校正算法處理Landsat8/OLI所得的443， 483， 561和655 nm波段數(shù)據(jù)進(jìn)行匹配(匹配原則是時(shí)間不超過3 h， 空間距離不超過300 m)， 共得到了1 853組匹配數(shù)據(jù)對(duì)， 得到了如圖9所示的對(duì)比結(jié)果。 其中(a)， (c)， (e)和(g)是Landsat8/OLI數(shù)據(jù)經(jīng)過Seadas Default大氣校正算法所得的443, 483, 561和655 nm波段遙感反射比數(shù)據(jù)與Modis數(shù)據(jù)的對(duì)比結(jié)果圖， (b)， (d)， (f)和(h)是同時(shí)空條件下的Landsat8/OLI數(shù)據(jù)經(jīng)過Acolite DSF大氣校正算法所得的443, 483,561和655 nm波段遙感反射比數(shù)據(jù)與Modis數(shù)據(jù)的對(duì)比結(jié)果圖。 計(jì)算其相關(guān)系數(shù)r、 均方根誤差RMSE， 并定義平均百分比偏差表達(dá)式ε[式(4)]對(duì)其進(jìn)行誤差分析(表4)。 結(jié)果表明， 在各個(gè)波段的Seadas Default 所得的大氣校正結(jié)果都要優(yōu)于Acolite DSF。

圖8 Landsat8/OLI數(shù)據(jù)分別經(jīng)Seadas Default和Acolite DSF大氣校正算法所得的

圖9 由Seadas Default算法得到的Landsat8/OLIRrs數(shù)據(jù)與Modis Rrs數(shù)據(jù)在443, 483, 561, 655 nm波段的對(duì)比圖[(a)， (c)， (e)， (g)]； 以及由Acolite DSF算法得到的Landsat8/OLI Rrs數(shù)據(jù)與Modis Rrs數(shù)據(jù)的對(duì)比圖[(b)， (d)， (f)， (h)]

(4)

表4 圖9中各對(duì)比結(jié)果的誤差參數(shù)

4 結(jié) 論

水體的遙感反射比光譜數(shù)據(jù)是海洋水色遙感反演海洋生物地球光學(xué)參數(shù)的關(guān)鍵。 對(duì)傳感器數(shù)據(jù)進(jìn)行大氣校正是我們得到精確的水體遙感反射比光譜數(shù)據(jù)的關(guān)鍵因素。 Landsat的高空間分辨率使得它在海洋遙感監(jiān)測(cè)應(yīng)用領(lǐng)域有著不可替代的優(yōu)勢(shì)。 目前支持Landsat數(shù)據(jù)進(jìn)行大氣校正的算法有很多， 哪種算法可以得到較高質(zhì)量的遙感發(fā)射比光譜， 需要我們對(duì)其進(jìn)行評(píng)價(jià)分析。 基于QA Score光譜質(zhì)量評(píng)價(jià)體系對(duì)Landsat數(shù)據(jù)處理中的五種大氣校正算法在膠州灣海域的適應(yīng)性進(jìn)行了評(píng)價(jià)分析。 結(jié)果表明在膠州灣海域Seadas的NASA標(biāo)準(zhǔn)大氣校正算法(Seadas Default)得到的Rrs數(shù)據(jù)QA得分為1的概率都要高于其他四種。 Acolite的暗光譜擬合算法(Acolite DSF)優(yōu)于Acolite的指數(shù)外推的三種大氣校正算法(Acolite SWIR， Acolite Red/NIR， Acolite NIR/SWIR)。 Acolite的指數(shù)外推的三種大氣校正算法在該海域內(nèi)對(duì)Landsat8/OLI數(shù)據(jù)的大氣校正效果較差。 Landsat8/OLI分別經(jīng)Seadas Default和Acolite DSF大氣校正算法所得的443, 483, 561, 655nm波段遙感反射比數(shù)據(jù)與Modis數(shù)據(jù)的比較結(jié)果(圖9)， 也從定量分析方面印證了SeadasNASA標(biāo)準(zhǔn)大氣校正算法在膠州灣海域?qū)andsat8/OLI數(shù)據(jù)的大氣校正優(yōu)于Acolite DSF。 鑒于上述研究結(jié)果， 建議在膠州灣及其鄰近海域應(yīng)用Landsat8/OLI數(shù)據(jù)進(jìn)行生物光學(xué)算法的遙感應(yīng)用研究實(shí)踐時(shí)， 首選Seadas的NASA標(biāo)準(zhǔn)大氣校正算法。 本研究結(jié)果也為研究人員在其他海域?qū)πl(wèi)星數(shù)據(jù)進(jìn)行大氣校正處理的算法選擇提供了新的判定思路。