基于GOCI 數據黃海氣溶膠光學厚度反演方法研究

鄭小慎，張亞男，劉雨華

（天津科技大學 海洋與環境學院，天津 300457）

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

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

1 數據來源

1.1 地基數據

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

圖1 AERONET 全球站位分布圖

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

1.2 GOCI 數據

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

2 反演方法

2.1 GPRO 算法

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

圖2 GPRO 算法流程圖

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

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

2.2 GNFW 算法

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

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

圖3 基于GOCI 數據的氣溶膠光學特性反演流程圖

2.2.1 表觀反射率的反演

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

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

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

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

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

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

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

2.2.2 離水反射率反演

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

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

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

GNFW 算法中使用的是MODIS 離水反射率每天的數據產品，這相比于GPRO 算法中通過30 天影像處理，而獲取離水反射率更高效且精度更高。MODIS 數據的處理需要與GOCI 數據在時間上保持對應，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 之間經過黃海上空的遙感數據進行處理分析。

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

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

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

當確定輸入圖像的噪聲像元和非噪聲像元位置之后，對原始數據按照5*5 窗口進行搜索，若目標像元為噪聲像元，利用窗口內非噪聲像元的均值代替噪聲像元原有的值。若目標像元為非噪聲像元仍然維持原值。

2.2.3 輻亮度反演

將處理后反射率轉化為輻亮度，輸出結果將為后續利用查找表確定最終反演結果做準備。

2.2.4 氣溶膠光學厚度反演

通過以上操作可以得到GOCI 數據8 個波段對應的去掉離水反射率之后轉換得到的輻亮度結果。

衛星遙感的輻亮度大氣校正方程可表達為

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

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

其中，ωA（λ）為氣溶膠單次散射比，ωA（λ）=KAS（λ）/KAE（λ），KAS為氣溶膠散射系數，KAE為氣溶膠消光系數，取經驗值0.8。τA（λ）表示氣溶膠光學厚度，FS（λ）為衰減的太陽輻亮度，PA（λ，θ，θs）為氣溶膠相函數，θs為太陽天頂角，θ 為觀測角。

2.2.5 查找表搜索

2.2.5.1 數據輸入

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

2.2.5.2 查找表的建立和搜索

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

表2 Rstar5b 氣溶膠基本模態組成

由表2 可以看出，氣溶膠模態在內部混合比例上為固定值。

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

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

表3 Rstar5b 氣溶膠模態外部混合組成

模型中用FMF 表示混合模態（Type=AB）中A 的混合比例，也就說，FMF 越大，混合模態中A的比例越高。

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

表4 AOD'和FMF 的輸入數據

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

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

圖4 LUT 搜索流程圖

3 反演結果驗證

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

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

圖5 Anmyon 站位AERONET 數據與GPRO 算法和GNEW 算法反演結果對比

如圖5 所示，在夏季，GNFW 算法相較于GPRO 算法更接近AFRONFT 數據，并且GPRO 算法在AOD 大于0.25 時有明顯的低估。在冬季，在整個監測區間內GNFW 算法比GPRO 算法更接近AFRONFT 數據，并且主要的氣溶膠光學厚度主要集中在0.02～0.3 區間內。冬季當氣溶膠光學厚度在0.02～0.15 區間時，GNFW 算法與GPRO 算法都存在高估的情況，但是GNFW 算法仍更接近AFRONFT 數據。當AOD 大于0.15 時，GNFW 算法與GPRO 算法都存在低估的情況，但是GNFW算法仍更接近AFRONFT 數據。總體而言，GNFW算法相比于GPRO 算法更接近AFRONFT 數據。

表5 GNEW 算法和GPRO 算法的相關統計數據

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

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

4 結論

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