基于GOCI 數(shù)據(jù)黃海氣溶膠光學(xué)厚度反演方法研究

鄭小慎，張亞男，劉雨華

（天津科技大學(xué) 海洋與環(huán)境學(xué)院，天津 300457）

大氣氣溶膠是大氣中的微量成分，在大氣中氣溶膠的含量雖然不高，但對調(diào)控地球輻射收支平衡及全球氣候模式有著十分重要的作用，氣溶膠成為研究大氣污染和氣候變化的重要參數(shù)。海洋上空氣溶膠光學(xué)性質(zhì)的研究對于估計電磁波的輻射、分析全球氣候變化和提高定量遙感精度有重要意義（劉玉光，2009；沈永平 等，2013）。Griggs（1975）和Fraser（1976）利用單通道原理建立的海洋氣溶膠光學(xué)厚度（Aerosol Optical Depth，AOD）反演算法，對于海表反射率較大的區(qū)域反演效果并不理想。多通道遙感算法反演氣溶膠主要是利用通道反射特征避免下墊面反射的影響來提高反演精度。Durkee 等（1986）用AVHRR 的630 nm 和830 nm波段來反演AOD，Remer 等（2005）驗證了660 nm和870 nm MODIS 海上AOD 在全球范圍內(nèi)的反演有效性，結(jié)果證明MODIS 數(shù)據(jù)的反演效果很好。劉毅等（2008）利用MODIS 氣溶膠業(yè)務(wù)產(chǎn)品驗證渤海、黃海、東海和日本以南海域的反演結(jié)果，分析結(jié)果表明基于MODIS 數(shù)據(jù)的雙通道算法和業(yè)務(wù)算法對黃海的反演結(jié)果并不理想。多角度氣溶膠遙感的原理是基于大氣和地表對大氣頂層衛(wèi)星信號貢獻(xiàn)的比率隨觀測角度的不同而不同的特點，將地表與大氣信號區(qū)分開來。Veefkind 等（1998；2000）利用雙角度算法（ATSR-DV）算法反演了不同區(qū)域的氣溶膠并得到了較好的反演結(jié)果。

針對黃海海域，氣溶膠光學(xué)特性的反演主要受到以下三點因素的影響，（1）黃海海域上空不但受到海洋水汽蒸發(fā)的影響，也受到包括沿岸工業(yè)排放、化石燃料及西北內(nèi)陸黃沙等陸源物質(zhì)的影響。因此，其上空的氣溶膠包括非吸收性氣溶膠和吸收性氣溶膠。（2）區(qū)域下墊面受到內(nèi)陸河流影響較為明顯，近岸海域為典型的渾濁水體（牟冰，2014；王正，2013），下墊面的反射影響反演精度。（3）由于氣溶膠的空間分布不均勻，時間變化較大，低時空分辨率的遙感數(shù)據(jù)不能在實際觀測應(yīng)用中起到作用。因此，減少以上因素的影響有助于提高黃海上空氣溶膠光學(xué)特性的反演精度。在此背景之下，本文提出了一種針對黃海氣溶膠光學(xué)厚度的新反演算法。該算法基于高時空分辨率衛(wèi)星GOCI數(shù)據(jù)，通過增加吸收性氣溶膠模態(tài)，減小由于水體反射引起的誤差，達(dá)到提高反演精度的目的。

1 數(shù)據(jù)來源

1.1 地基數(shù)據(jù)

AFRONFT（Aerosol Robotic Network） 是一個由NASA 和PHOTONS 共同建立的氣溶膠地基遙感觀測網(wǎng)絡(luò)（郭陽潔，2014），該網(wǎng)絡(luò)現(xiàn)已經(jīng)覆蓋了全球主要的區(qū)域，目前全球共有超過500 個站點。通常，將AFRONFT 所測得的氣溶膠光學(xué)厚度作為真值，來檢驗遙感反演的氣溶膠光學(xué)厚度的準(zhǔn)確性，圖1 為AFRONFT 全球站位圖。

圖1 AERONET 全球站位分布圖

AFRONFT 氣溶膠數(shù)據(jù)包括1 級、1.5 級和2級數(shù)據(jù)產(chǎn)品， 1 級產(chǎn)品未經(jīng)篩選和最終的校正；1.5 級產(chǎn)品經(jīng)過了云污染篩選，未經(jīng)過最終校正；2級產(chǎn)品經(jīng)過了校正、云污染篩選和質(zhì)量控制。黃海AFRONFT 站位少，且部分站位數(shù)據(jù)的可用性差，實驗選擇有數(shù)據(jù)質(zhì)量保證且位于GOCI 刈幅內(nèi)的Anmyon 站位（36.54°N，126.33°F），能夠驗證反演算法對渾濁水體的適用性。下載網(wǎng)址為https://aeronet.gsfc.nasa.gov。

1.2 GOCI 數(shù)據(jù)

2010 年韓國成功發(fā)射地球靜止軌道衛(wèi)星COMS，這顆衛(wèi)星上搭載著世界上第一個靜止軌道GOCI 海洋水色傳感器，每天能夠提供從00:16 UTM到7:16 UTM 的分辨率為500 m 的8 幅衛(wèi)星遙感圖像。由于時空分辨率較高，GOCI 數(shù)據(jù)對于監(jiān)測黃海的海洋環(huán)境變化十分有意義，本文選用的GOCI遙感數(shù)據(jù)包括GOCI L1B 數(shù)據(jù)以及GDPS 處理的L2A 級數(shù)據(jù)，GOCI 數(shù)據(jù)直接從韓國COMS 衛(wèi)星官網(wǎng)上免費下載，網(wǎng)站為http://kosc.kiost.ac.kr/。

2 反演方法

2.1 GPRO 算法

GOCI 傳感器的優(yōu)勢在于利用靜止軌道衛(wèi)星獲取高時空分辨率的遙感影像，GPRO 算法為AOD反演的業(yè)務(wù)化方法，通過去除離水反射率的影響反演AOD。反演流程如圖2 所示。

圖2 GPRO 算法流程圖

對于渾濁水體，GPRO 僅利用去除離水反射率的490 nm 反射率反演550 nm 處的AOD。與AFRONFT 數(shù)據(jù)比較，GPRO 反演算法相對于MODIS 業(yè)務(wù)算法的反演精度更高（Lee et al，2010）。

GPRO 反演算法主要有三點不足：第一，對渾濁水體采用30 d 清潔窗確定離水反射率存在一定的誤差；第二，GPRO 算法中并沒有具體給出氣溶膠Type；第三，GOCI 業(yè)務(wù)化算法對云去除過程過于嚴(yán)格導(dǎo)致在渾濁水體上空反演會出現(xiàn)偏差。

2.2 GNFW 算法

針對GPRO 算法中對離水反射率的處理以及查找表中氣溶膠模態(tài)設(shè)置中存在的缺陷，本文提出一種基于GOCI 數(shù)據(jù)的氣溶膠光學(xué)特性反演算法（GNFW 算法）。該算法的主要流程如圖3 所示。

該算法主要由三部分組成，第一部分為表觀反射率與輻亮度之間的相互轉(zhuǎn)換反演，為了去除由于太陽天頂角及觀測點與太陽之間距離的變化而引起的誤差，該算法首先將GOCI 的輻亮度數(shù)據(jù)轉(zhuǎn)化為表觀反射率。該算法利用去除離水反射率的反射數(shù)據(jù)來減少由于水體輻射信號引起的干擾。第二部分為離水反射率除噪聲反演，針對離水反射率產(chǎn)品存在明顯條帶噪聲的問題，該算法提出了插值法去條帶噪聲以提高離水反射率的準(zhǔn)確性和可用性。第三部分為查找表的建立和搜索。該算法的輸出結(jié)果為氣溶膠光學(xué)厚度（AOD），細(xì)顆粒比例（Fraction of AOD contributed by fine dominated model，F(xiàn)MF） 及氣溶膠模態(tài)（Type）。

圖3 基于GOCI 數(shù)據(jù)的氣溶膠光學(xué)特性反演流程圖

2.2.1 表觀反射率的反演

基于下載的GOCI L1B 數(shù)據(jù)，根據(jù)大氣校正方程（張民偉，2009）將衛(wèi)星接收的輻亮度轉(zhuǎn)化為表觀反射率，式（1） 為表觀反射率的計算公式（Teillet et al，2001），

其中，ρ（λ）表示λ 波長處的表觀反射率；d 表示日地距離（天文單位）；θ 表示太陽天頂角；ESUN（λ）表示λ 波長處平均太陽輻照度；L（λ）表示λ 波長處傳感器接收到的輻亮度，單位為W·m-2·sr-1·μm-1。

上式中日地距離是遙感影像成像時間的函數(shù)，表達(dá)式如下：

其中，JD 表示遙感成像儒略日。

波段平均太陽輻照度是大氣層外的太陽光譜輻照度與傳感器各波段光譜響應(yīng)函數(shù)的積分，表達(dá)式如下：其中，S（λ）表示λ 波長處的光譜響應(yīng)函數(shù)，BMSI（λ）表示λ 波長處大氣層外太陽光譜輻照度。

采用Wehrli85 太陽光譜曲線計算GOCI 接收的波段平均太陽輻照度（張楠，2013；張璐 等，2014）。Wehrli85 太陽光譜范圍為0.2～10 μm，GOCI 的波譜響應(yīng)區(qū)間主要集中在300～1 000 nm 范圍內(nèi)，因此，對300～1 000 nm 區(qū)間內(nèi)的大氣層外太陽光譜輻照度和光譜響應(yīng)函數(shù)進行積分，最終獲得GOCI 傳感器接收的大氣層外平均太陽光譜輻照度（Mean solar exoatmospheric irradiances，F(xiàn)SUN），如表1 所示。

表1 GOCI 接收的平均太陽輻照度（W·m-2·sr-1·μm-1）

2.2.2 離水反射率反演

大洋水體下墊面比較均勻，離水輻射在可見光和近紅外波段受到水體的影響較小，針對大洋水體的氣溶膠反演已經(jīng)很成熟并且反演精度較高。但是渾濁水體受到的影響比較復(fù)雜，影響其上空氣溶膠散射的主要因素是離水反射率，計算離水反射率成為氣溶膠光學(xué)特性反演的重要過程。

2.2.2.1 NIR-SWIR 算法反演離水反射率

對于渾濁水體區(qū)域，NIR 離水反射率算法處理后會出現(xiàn)明顯的誤差（Lavender et al，2005；Ruddick et al，2000；Siegel et al，2000），利 用SWIR 在近岸渾濁水體反射率基本為零的特性來反演離水反射率更為準(zhǔn)確。GOCI 沒有SWIR 波段，用MODIS 數(shù)據(jù)來彌補GOCI 波段設(shè)置的缺陷，使用NIR-SWIR 算法提高反演離水反射率精度，間接提高AOD 的反演精度。利用渾濁水體系數(shù)對水體進行區(qū)分，若為渾濁水體采用SWIR 波段進行處理，若為非渾濁水體采用NIR 波段進行處理。

GNFW 算法中使用的是MODIS 離水反射率每天的數(shù)據(jù)產(chǎn)品，這相比于GPRO 算法中通過30 天影像處理，而獲取離水反射率更高效且精度更高。MODIS 數(shù)據(jù)的處理需要與GOCI 數(shù)據(jù)在時間上保持對應(yīng)，GOCI 中心掃描時間為0:16、1:16、2:16、3:16、4:16、5:16、6:16、7:16 UTC，MODIS Terra 通過黃海上空的時間為1:30～3:30 UTC 之間，因此，實驗采用MODIS Terra 1:30～3:30 UTC 之間經(jīng)過黃海上空的遙感數(shù)據(jù)進行處理分析。

2.2.2.2 離水反射率反演去噪聲處理

MODIS NIR-SWIR 算法處理后的離水反射率表現(xiàn)出明顯的條帶狀噪聲。采用插值法對MODIS離水反射率影像上的條帶噪聲進行處理，利用式（4）確定最終的噪聲像元位置。

其中，l 表示數(shù)據(jù)矩陣的行，*表示數(shù)據(jù)矩陣l 行的所有列元素，fl 表示最近非噪聲行，T 為檢測噪聲的閾值。

當(dāng)確定輸入圖像的噪聲像元和非噪聲像元位置之后，對原始數(shù)據(jù)按照5*5 窗口進行搜索，若目標(biāo)像元為噪聲像元，利用窗口內(nèi)非噪聲像元的均值代替噪聲像元原有的值。若目標(biāo)像元為非噪聲像元仍然維持原值。

2.2.3 輻亮度反演

將處理后反射率轉(zhuǎn)化為輻亮度，輸出結(jié)果將為后續(xù)利用查找表確定最終反演結(jié)果做準(zhǔn)備。

2.2.4 氣溶膠光學(xué)厚度反演

通過以上操作可以得到GOCI 數(shù)據(jù)8 個波段對應(yīng)的去掉離水反射率之后轉(zhuǎn)換得到的輻亮度結(jié)果。

衛(wèi)星遙感的輻亮度大氣校正方程可表達(dá)為

其中，L（λ）表示衛(wèi)星探測的輻亮度；LR（λ）表示由于大氣中分子瑞利散射引起的輻亮度；LA（λ）表示由于氣溶膠散射引起的輻亮度；LRA表示由于大氣中空氣分子與氣溶膠之間多次散射引起的輻亮度；Lr（λ）表示由于太陽耀斑引起的輻亮度；T（λ，θ）表示大氣的直接透射率；Lwc（λ）表示由于白帽引起的輻亮度；t（λ，θ）表示大氣的漫透射率；Lw（λ）表示離水輻亮度；λ 表示傳感器波長；θ 表示衛(wèi)星天頂角。

對于GOCI 數(shù)據(jù)8 個可見光波段的單次散射來說，如果不考慮多次散射引起的輻亮度，消除了每個波段中瑞利散射，太陽耀斑，白帽和離水輻亮度，氣溶膠散射引起的輻亮度簡化為：

其中，ωA（λ）為氣溶膠單次散射比，ωA（λ）=KAS（λ）/KAE（λ），KAS為氣溶膠散射系數(shù)，KAE為氣溶膠消光系數(shù)，取經(jīng)驗值0.8。τA（λ）表示氣溶膠光學(xué)厚度，F(xiàn)S（λ）為衰減的太陽輻亮度，PA（λ，θ，θs）為氣溶膠相函數(shù)，θs為太陽天頂角，θ 為觀測角。

2.2.5 查找表搜索

2.2.5.1 數(shù)據(jù)輸入

Rstar5b 是由 Nakajima 和Tanaka（Nakajima et al，1986；1988；Stamnes et al，1988）提出的用于模擬大氣—陸地—海洋之間輻射傳輸過程的輻射傳輸模型。該模型曾用于ADFOS/OCTS 及ADFOS-Ⅱ/GLI 產(chǎn)品反演算法中。Rstar5b 中假定陸地或海洋之上的平行平面為若干均勻?qū)?，模擬波段范圍為0.2～200 μm。

2.2.5.2 查找表的建立和搜索

Rstar5b 中包括水蒸氣、冰粒子、類塵埃、水溶性粒子、海鹽粒子、煙塵粒子、火山灰粒子、75 %H2SO4、黃沙在內(nèi)的9 種基本的氣溶膠顆粒。氣溶膠基礎(chǔ)模態(tài)由這9 種基本氣溶膠顆粒按照不同比例進行內(nèi)部混合，Rstar5b 中氣溶膠基本模態(tài)組成比例如表2 所示。

表2 Rstar5b 氣溶膠基本模態(tài)組成

由表2 可以看出，氣溶膠模態(tài)在內(nèi)部混合比例上為固定值。

實驗中針對黃海上空氣溶膠模態(tài)的分析，選用類塵埃、煙塵、海鹽、城市、黃沙這五種氣溶膠基本模態(tài)進行外部混合。氣溶膠的外部混合方式按照表3 所示的方式，其中√表示對該基本模態(tài)進行混合，空白表示不進行處理。

如表所示，分別用代碼3、4、8、9、11 表示類塵埃、煙塵、海鹽、城市、黃沙氣溶膠模態(tài)，外部混合方式為將代碼依次連接起來，依次表示為Type=38，Type=39，Type=48，Type=49，Type=98，Type=811，Type=911。

表3 Rstar5b 氣溶膠模態(tài)外部混合組成

模型中用FMF 表示混合模態(tài)（Type=AB）中A 的混合比例，也就說，F(xiàn)MF 越大，混合模態(tài)中A的比例越高。

因此，Rstar5b 模型的輸入變量可根據(jù)以上說明建立，特別說明的是其中方位信息可根據(jù)圖像角度文件獲得，對于輸入變量中的自變量AOD＇及FMF 的組成分別如表4 所示。

表4 AOD'和FMF 的輸入數(shù)據(jù)

如表5 所示，Rstar5b 的輸入文件個數(shù)取決于太陽天頂角、出射光線角、相對方位角、AOD＇以及FMF 的個數(shù)。由于以上參數(shù)進行組合后數(shù)據(jù)量極大，因此實驗中對原Rstar5b 模型進行修改，使得模型可以循環(huán)運行并逐個輸出對應(yīng)文件。

查找表（Look up table，LUT） 建立時對于給定角度參數(shù)，每種LUT 包括11 種不同的FMF 參數(shù)，每種FMF 中包括10 種不同的AOD＇。查找表搜索流程如圖4 所示。在確定像元在LUT 中的起始位置后，對每種LUT 提取相應(yīng)的110 行進行逐個計算，確定最佳LUT。實驗中設(shè)計的氣溶膠模態(tài)有7 種，將上述過程循環(huán)7 次獲取7 個相應(yīng)的FMF、490 nm AOD 及490 nm 輻亮度值。將每種LUT 對應(yīng)的490 nm 處輻亮度插值結(jié)果與觀測輻亮度進行對比，將最小差值作為最佳的LUT，同時輸出相應(yīng)LUT 的FMF 值和AOD 值。因此，輸出結(jié)果包括每個像元的FMF、AOD 及氣溶膠模態(tài)。

圖4 LUT 搜索流程圖

3 反演結(jié)果驗證

基于Anmyon 站位的AFRONFT 數(shù)據(jù)對GNFW算法、GPRO 算法的反演結(jié)果進行了驗證，其中GNFW 為本文提出的基于GOCI 數(shù)據(jù)氣溶膠光學(xué)厚度反演算法，GPRO 算法為GOCI 數(shù)據(jù)業(yè)務(wù)化AOD反演算法?；贕OCI 數(shù)據(jù)時間與質(zhì)量，選擇Anmyon 站位無云天氣狀況下，2014 年和2015 年夏季（6-8 月）、2015 年和2014 年冬季（12-2 月）1.5 級AFRONFT 數(shù)據(jù)進行驗證。

在Anmyon 站位基于AFRONFT 數(shù)據(jù)分別對比GPRO 反演數(shù)據(jù)和GNFW 反演數(shù)據(jù)，結(jié)果如圖5所示，（a）圖為夏季，（b）圖為冬季。

圖5 Anmyon 站位AERONET 數(shù)據(jù)與GPRO 算法和GNEW 算法反演結(jié)果對比

如圖5 所示，在夏季，GNFW 算法相較于GPRO 算法更接近AFRONFT 數(shù)據(jù)，并且GPRO 算法在AOD 大于0.25 時有明顯的低估。在冬季，在整個監(jiān)測區(qū)間內(nèi)GNFW 算法比GPRO 算法更接近AFRONFT 數(shù)據(jù)，并且主要的氣溶膠光學(xué)厚度主要集中在0.02～0.3 區(qū)間內(nèi)。冬季當(dāng)氣溶膠光學(xué)厚度在0.02～0.15 區(qū)間時，GNFW 算法與GPRO 算法都存在高估的情況，但是GNFW 算法仍更接近AFRONFT 數(shù)據(jù)。當(dāng)AOD 大于0.15 時，GNFW 算法與GPRO 算法都存在低估的情況，但是GNFW算法仍更接近AFRONFT 數(shù)據(jù)?？傮w而言，GNFW算法相比于GPRO 算法更接近AFRONFT 數(shù)據(jù)。

表5 GNEW 算法和GPRO 算法的相關(guān)統(tǒng)計數(shù)據(jù)

表5 列舉了GNFW 算法和GPRO 算法在夏季和冬季反演結(jié)果的相關(guān)統(tǒng)計數(shù)據(jù)，由于驗證GOCI和AFRONFT 的數(shù)據(jù)是隨機選擇的，因此結(jié)果可以說明GNFW 算法在冬季和夏季都表現(xiàn)出較GPRO算法更好的反演效果。綜上所述，GNFW 算法在渾濁水體上空能較為準(zhǔn)確地反演AOD，同時GNFW算法的反演精度相比于GPRO 算法有提高。

GNFW 算法誤差產(chǎn)生的可能原因，主要包括以下三點，第一，輻射傳輸模型中大氣及氣溶膠顆粒參數(shù)設(shè)定不準(zhǔn)確，這其中包含大氣組分濃度、大氣壓力廓線、大氣質(zhì)量、溫度廓線、粒子尺度分布等（王欣睿等，2016）；第二，由于水體環(huán)境變化使得離水反射率隨時間出現(xiàn)變化。第三，去除離水反射率時對大氣漫透射率的估計出現(xiàn)誤差。

4 結(jié)論

黃海海域氣溶膠光學(xué)特性反演精度受到陸源物質(zhì)的影響顯著，近岸水體上空的氣溶膠模態(tài)為吸收性氣溶膠和非吸收性氣溶膠共存的狀態(tài)。本文在GOCI 數(shù)據(jù)業(yè)務(wù)化氣溶膠反演算法GPRO 算法基礎(chǔ)上，提出了改進的GNFW 算法。針對黃海海域氣溶膠光學(xué)特性的反演，該算法基于Rstar5b 氣溶膠模態(tài)內(nèi)部比例，選用類塵埃、煙塵、海鹽、城市、黃沙這五種氣溶膠基本模態(tài)進行外部混合形成黃海上空氣溶膠模態(tài)類型，兼顧了黃海上空非吸收性氣溶膠和吸收性氣溶膠共存的狀態(tài)。對于渾濁水體，利用MODIS 波段設(shè)置的優(yōu)勢，進行離水反射率反演，進行GOCI 數(shù)據(jù)中離水反射率的去除，對于氣溶膠反演精度的提高有重要作用。GOCI 數(shù)據(jù)的高時空分辨率在氣溶膠觀測中起到了重要作用。該算法能較好地反演渾濁水體上空的氣溶膠光學(xué)厚度，相比于GPRO 算法在反演精度上有所提高。