PSO聚類算法的京津冀地區(qū)氣溶膠光學厚度反演

王書濤，王貴川，凡堃堃，吳 興，王玉田

燕山大學電氣工程學院，河北省測試計量技術(shù)及儀器重點實驗室，河北 秦皇島 066004

引 言

近年來，遙感衛(wèi)星影像反演已經(jīng)成為了獲取較大時空分布的氣溶膠光學厚度(aerosol optical depth, AOD)的主要手段，國內(nèi)外學者根據(jù)不同星載傳感器數(shù)據(jù)的特點和成像區(qū)域特征提出很多算法：暗目標法[1-2]、深藍算法[3]、多星協(xié)同算法[4]等。其中，暗目標法充分利用了遙感影像中濃密植被暗像元的信息特征，受到地表噪聲的干擾較小，反演流程也比較簡單，因此得到了廣泛應(yīng)用。我國學者對暗目標法進行改進以用于國產(chǎn)衛(wèi)星數(shù)據(jù)反演AOD，如HJ系列衛(wèi)星[5]、GF系列衛(wèi)星[6]、CBERS02B衛(wèi)星[7]等。

氣溶膠理化性質(zhì)的不確定性和時空分布的無規(guī)律性是從遙感影像中反演氣溶膠參數(shù)的難點。AOD反演的核心是如何剔除地表反射率貢獻，通過輻射模型如6S模型、MODTRAN模型等能夠?qū)崿F(xiàn)這一目標。模型的輸入?yún)?shù)中，氣溶膠類型是影響AOD反演精度的主要因素之一，由于不同地域下地理環(huán)境、氣候和社會發(fā)展情況的差異，模型中預(yù)設(shè)氣溶膠類型的針對性和適用性較差，所以確定研究區(qū)域氣溶膠類型是準確反演AOD的關(guān)鍵。

針對這些問題，提出了一種新方法：一方面調(diào)整暗目標法以適用于高分四號(GF-4)多光譜遙感數(shù)據(jù)，構(gòu)建類三元一次方程組；另一方面利用粒子群(particle swarm optimization, PSO)聚類算法確定符合京津冀地區(qū)的氣溶膠類型，將其作為6S模型的輸入?yún)?shù)建立查找表，然后逐像元查表求解方程組，實現(xiàn)對研究區(qū)域AOD的反演。為了校驗此方法AOD的反演精度，分別利用MODIS氣溶膠產(chǎn)品(MOD04_L2)和AERONET(AErosol RObotic NETwork)地基觀測數(shù)據(jù)與AOD反演結(jié)果進行對比，并進行分析和討論。

1 實驗部分

1.1 算法流程及暗目標法的調(diào)整

算法流程如圖1所示，首先針對暗目標法不適用于GF-4數(shù)據(jù)的情況，對GF-4/PMI傳感器的數(shù)據(jù)進行了分析，并根據(jù)已有數(shù)據(jù)和資料對暗目標法進行調(diào)整。利用長時序的AERONET地基數(shù)據(jù)和PSO聚類算法對京津冀地區(qū)的氣溶膠類型進行了聚類分析，并結(jié)合實驗驗證選取了符合實際情況的氣溶膠類型。將聚類得到的氣溶膠類型代入6S輻射傳輸模型建立查找表，結(jié)合調(diào)整之后的暗目標法建立類三元一次方程組。最后選取京津冀區(qū)域GF-4衛(wèi)星影像進行AOD反演，并與MOD04_L2產(chǎn)品和AERONET地基觀測數(shù)據(jù)進行對比來驗證算法的可行性。

圖1 算法流程圖Fig.1 Algorithm flow diagram

GF-4衛(wèi)星是我國第一顆靜止軌道遙感衛(wèi)星，衛(wèi)星采用地球同步軌道，其位與赤道106°E上空36 000 km處，作為我國高分辨率對地觀測重大專項的衛(wèi)星之一，GF-4衛(wèi)星上的PMI傳感器兼具中空間分辨率和高時間分辨率的成像特點，相對于其他同步軌道衛(wèi)星，不僅具有更加細膩的成像效果，還能對某一大面積區(qū)域進行長期而連續(xù)的觀測。GF-4衛(wèi)星上搭載的PMI傳感器參數(shù)如表1所示。

表1 GF-4衛(wèi)星PMI傳感器參數(shù)Table 1 Parameters of PMI sensor for GF-4 satellite

假設(shè)地球為均勻郎伯表面，則星載傳感器接收到的大氣頂層反射率(即表觀反射率)ρTOA可以表示為[8]

(1)

式(1)中，θs為太陽天頂角，θv為衛(wèi)星天頂角，φ為相對方位角；ρ0是路徑輻射項等效反射率；ρs為地表二項反射率；S為大氣下界的半球反照率；T(μs)和T(μv)分別是向下和向上整層的大氣透過率，兩者的乘積代表大氣總的透射率。

表觀反射率受到大氣環(huán)境下的等效反射率和地表反射率影響。因此，地表反射率ρs可以根據(jù)表觀發(fā)射率ρTOA計算得到，簡化如式(2)

(2)

暗目標法也叫濃密植被法，傳統(tǒng)暗目標法利用短波紅外波段與紅藍波段之間地表反射率之間的相關(guān)關(guān)系。由表1可知GF-4衛(wèi)星PMI傳感器缺少暗目標法所需的短波紅外波段，因此需要對暗目標法進行調(diào)整和改進。假設(shè)GF-4衛(wèi)星PMI多光譜數(shù)據(jù)中暗像元的紅藍波段地表反射率間具有良好的相關(guān)性，即

(3)

式(3)中，k為相關(guān)系數(shù)；b為相應(yīng)偏差。

利用2018年9月內(nèi)蒙古試驗場地測得的植被光譜數(shù)據(jù)，對其中212條不同的植被光譜數(shù)據(jù)分析，得到PMI傳感器紅藍波段等效地表反射率的散點分布和擬合結(jié)果如圖2所示，擬合之后兩者的關(guān)系為f(x)=1.643x-0.007，相關(guān)系數(shù)達0.926。

圖2 紅藍波段等效地表反射率的擬合曲線Fig.2 Fit curve of red-blue band equivalent surface reflectivity

(4)

1.2 粒子群聚類算法

粒子群優(yōu)化算法是由Kennedy和Eberhar對鳥群覓食行為的研究中分析總結(jié)出來的尋優(yōu)算法。算法的原理如下：在預(yù)設(shè)搜索空間中初始化一群具有潛在最優(yōu)解的粒子，這些粒子以一定的速度飛行，粒子搜索的方向和距離由這些粒子的速度決定，粒子的好壞由適應(yīng)度值篩選，粒子在搜索空間中飛行時，會結(jié)合自身和其他粒子的適應(yīng)度值不斷更新自身的適應(yīng)度值，完成目標尋優(yōu)[9]。

聚類可以理解為樣本的排列與組合，這種組合使得同一類中的樣本特征相似度最大，而不同類之間的樣本差異也最大。PSO聚類算法的過程中，每個粒子都看作是對樣本間聚類的劃分，都可以看成需求最佳聚類中心的可行解，因此可將粒子結(jié)構(gòu)描述為預(yù)定分類數(shù)K個聚類中心Ck=(ck1,ck2, …,ckN)(k=1, 2, …,K)的連接，并編碼構(gòu)造如式(5)

Xi=(c11,c12, …,c1N,c21,c22, …,c2N, …,ck1,ck2, …,ckN)

(5)

粒子i可看成K×N維向量，當所有粒子的位置確定時，其代表的聚類中心隨之確定。假設(shè)共有J個樣本，對粒子進行拆分，第j個樣本的所屬類的劃分采用最近鄰域法

(6)

(7)

則樣本Sj屬于k′類，Sk指的是第k個聚類中的所有樣本，numk是該類中的樣本數(shù)。對該粒子處的各樣本進行分類之后，計算適應(yīng)度函數(shù)值fitness(t)，fitness(t)采用類內(nèi)離散度之和Jc來表示

(8)

式(8)中，Skn指第k類中的第n個樣本，Ck指第k類的聚類中心，t是迭代次數(shù)。當粒子的適應(yīng)度值小于其個體和全局歷史最優(yōu)適應(yīng)度值時，則更新粒子的適應(yīng)度值。

1.3 氣溶膠類型對AOD的敏感性分析

6S模型模擬了輻射在太陽-地表目標-接收器這一傳輸過程。假設(shè)研究區(qū)域為單一地表，使用6S模型中的城市型和海洋型、以及OPAC數(shù)據(jù)集[10]中定義的大陸平均型和大陸污染型氣溶膠類型，研究不同氣溶膠類型對表觀反射率的影響，如圖3所示?？芍擜OD較小時，四種氣溶膠類型對表觀反射率影響的差異不大，當AOD較大時，差異逐漸增大。大陸平均型、大陸污染型和海洋型三種氣溶膠類型得到的表觀反射率均隨著AOD值的增大而逐漸增大，且增加的幅度不相同，但是城市型氣溶膠類型得到的表觀反射率隨AOD值的增大而減小。實驗結(jié)果說明氣溶膠類型對AOD的敏感程度有很大差異，這種差異來自于氣溶膠類型中各種類型粒子不同的吸收、散射特性。進一步證明，氣溶膠類型會對AOD反演造成很大誤差。

1.4 京津冀區(qū)域氣溶膠類型聚類分析

根據(jù)長時序持續(xù)觀測的AERONET地基數(shù)據(jù)發(fā)現(xiàn)，使用氣溶膠微物理特性和光化學特性參數(shù)共同表征氣溶膠類型，并利用聚類算法得出氣溶膠類型是可行的[11-12]。以京津冀地區(qū)為研究區(qū)域，選擇AERONET地基網(wǎng)絡(luò)中2015年1月1日到2017年1月1日的北京站點(39.977°N, 116.381°E)和香河站點(39.754°N, 116.962°E)氣溶膠的特性數(shù)據(jù)進行統(tǒng)計分析。

圖3 四種氣溶膠類型對表觀反射率的影響Fig.3 Effect of four aerosol types on apparent reflectivity

最佳聚類數(shù)量是聚類分析的關(guān)鍵參數(shù)之一。聚類數(shù)量m不斷增加時，使用PSO聚類算法充分迭代求得各樣本到自身所對應(yīng)聚類中心的歐氏距離之和的最小值，用能量E表示(即最優(yōu)適應(yīng)度fitness)，E與m的關(guān)系為E=f(m)，根據(jù)m對于E的敏感性和擬合關(guān)系進行分析，可得到最佳聚類數(shù)量m。圖4反映了兩者在實數(shù)平面坐標上的散點分布情況，將其放到對數(shù)坐標軸上(圖4右上方)，可知能量E和m直線擬合度較好，證明了在實數(shù)坐標空間內(nèi)E隨m變化而指數(shù)衰減。由圖4可以看出，聚類數(shù)量m足夠大時，新增聚類所占比例越來越小，能量E的衰減程度也越來越弱, 表明新增聚類對于氣溶膠類型差異性的影響越來越小。

圖4 能量E與聚類數(shù)目m的關(guān)系Fig.4 Relationship between energy E and cluster number m

引用半衰期序列來進一步觀測[13]，以E0=f(1)為基準，即所有的樣本只劃分為一種類型時對應(yīng)的能量，當能量依次減半時取各聚類數(shù)量對應(yīng)的E最為接近的情況，進而可以得到處于能量半衰期位置的序列p，p對應(yīng)聚類數(shù)目的下限。由于聚類分析之前要把所有的氣溶膠樣本歸一化，使得E大幅減少，因此選擇p=2，即第二個半衰期對應(yīng)的m=8作為最佳聚類數(shù)量。PSO聚類算法的參數(shù)設(shè)置如表2所示，根據(jù)PSO聚類算法得到8種氣溶膠類型聚類結(jié)果。試驗發(fā)現(xiàn)，這8種氣溶膠類型中的5種占比高達89.4%，滿足研究區(qū)域大多數(shù)的大氣情況，故只對聚類結(jié)果中的這5種進行討論。

表2 PSO聚類算法的參數(shù)Table 2 Parameters of PSO clustering algorithm

該地區(qū)氣溶膠粒子的體積尺度譜適合用雙模態(tài)對數(shù)分布描述[14-15]，所以描述氣溶膠物理特性的粒子譜分布參數(shù)主要有粒子半徑、體積濃度等，表征氣溶膠光學特性的參數(shù)有復(fù)折射指數(shù)、單次散射反照率以及不對稱因子。聚類所得5種氣溶膠類型體積尺度譜分布特征如圖5所示，復(fù)折射指數(shù)隨波長變化情況如圖6所示。C1，C2和C4這三種類型占據(jù)了該地區(qū)這兩年氣溶膠類型的70%。C1和C2的體積譜相似，二者中兩種模態(tài)粒子的體積濃度都比較小，而且二者細模態(tài)粒子峰值半徑均在0.15 μm左右，5種類型中最小。C2聚類結(jié)果占比最大，其細粒子和粗粒子的體積占比最低，復(fù)折射指數(shù)的虛部最大，表示很強的吸收能力，是清潔天氣的氣溶膠類型；C1的占比較C2小，其粗細粒子體積濃度比背景值均略有增加，而且粗粒子濃度增幅大于細粒子，為晴朗天氣的氣溶膠類型。C4和C5體積譜均近似于單峰分布，但表現(xiàn)相反，C4中細模態(tài)粒子占優(yōu)，所以細模態(tài)粒子的消光貢獻遠大于粗模態(tài)粒子，且復(fù)折射指數(shù)實部和虛部也明顯降低，即C4代表著霧霾、陰霾天氣；而C5體積尺度譜以粗模態(tài)為主導(dǎo)峰，這主要來源于冬季和春季的沙塵型氣溶膠，這種氣溶膠的實部較高，具有較強的散射能力。對于C3來說，其體積尺度譜呈現(xiàn)出駝型分布，即細粒子體積和粗粒子體積相當，濃度峰值均在0.15 μm3·μm-2以上，代表嚴重的灰霾天氣。

圖5 5種氣溶膠類型體積尺度譜分布Fig.5 Volumetric spectral distribution of five aerosol types

圖6 5種氣溶膠模式復(fù)折射指數(shù)(a)：實部；(b)：虛部Fig.6 (a) The birefringence index real part (a) and birefringence index imaginary part (b) of five types of aerosol

基于對5種氣溶膠類型特點的分析，C1，C2和C4型更符合京津冀地區(qū)夏、秋季節(jié)的氣溶膠類型特征。經(jīng)過實驗驗證，C1和C2型氣溶膠類型的反演結(jié)果相似，且二者的物理和光學特性非常相似，為此不分析C2，以C1，C4以及大陸型氣溶膠類型進行AOD的反演。

2 結(jié)果與討論

為了評價AOD反演結(jié)果的可信度，選取2017年9月12日的GF-4/PMI遙感數(shù)據(jù)進行AOD反演，GF-4衛(wèi)星的過境時間是10:47，MODIS數(shù)據(jù)的獲取時間是10:20，時間相差不超過30 min。將反演結(jié)果經(jīng)過投影轉(zhuǎn)換，并將分辨率重采樣到10 km，再與MOD04_L2進行比較分析，提取N=658個有效像元，去除異常值之后進行逐像元對比，反演結(jié)果和MOD04_L2產(chǎn)品對比及擬合結(jié)果如圖7所示。

由圖7可以看出，3種氣溶膠類型AOD反演結(jié)果在層次方面和變化方面與MOD04_L2產(chǎn)品基本保持一致。C4和大陸型氣溶膠反演的AOD情況較好，兩者均比MODIS氣溶膠產(chǎn)品獲取了更多細節(jié)性的信息，相關(guān)系數(shù)分別達到了0.794和0.811，均方根誤差也分別為0.260 7和0.122 4。而C4類型反演結(jié)果均方根誤差較大的原因在于其對AOD值變化的敏感特性，即氣溶膠濃度越高，AOD值越大，反演出來的結(jié)果也越與實際結(jié)果相符合。這充分體現(xiàn)了C4氣溶膠類型具有良好的張性和波動性。C4型氣溶膠粒子中的細粒子占據(jù)主導(dǎo)消光貢獻，隨著空氣中相對濕度的增大，氣溶膠粒子呈現(xiàn)吸濕性增長，增強了氣溶膠粒子的散射性能，此時C4型較其他兩種類型更具優(yōu)勢。

圖7 GF-4衛(wèi)星氣溶膠光學厚度反演結(jié)果及一致性分析(a)：原始數(shù)據(jù)；(b)：C1型反演結(jié)果；(c)：C4型反演結(jié)果；(d)：大陸型反演結(jié)果；(e)：MODIS溶膠產(chǎn)品；(f)：C1型反演結(jié)果分析；(g)：C4型反演結(jié)果分析；(h)：大陸型反演結(jié)果分析Fig.7 Inversion results and consistency of the aerosol optical depth of GF-4 satellite(a)：Original image；(b)：C1 aerosol；(c)：C4 aerosol；(d)：Continental aerosol；(e)：MODIS aerosol products；(f) Analysis of C1 type aerosol inversion results；(g)：Analysis of C4 type aerosol inversion result；(h)：Analysis of continental aerosol inversion results

AERONET地基觀測網(wǎng)使用的CET-318太陽光度計只提供波長340，380，440，500，675，870，1 020和1 640 nm處的AOD值，Angstrom[16]研究發(fā)現(xiàn)，AOD值與波長之間存在著指數(shù)衰減的關(guān)系，可計算得到550 nm波長處的AOD，并選取GF-4衛(wèi)星過境時間前后30min的插值結(jié)果求取平均值，當作AERONET的AOD觀測數(shù)據(jù)。另外，以影像中的地基站點為中心，周圍5 km×5 km像元窗口下的AOD反演結(jié)果求取平均值，作為GF-4衛(wèi)星多光譜影像的AOD反演結(jié)果。京津冀區(qū)域2018年不同時間不同氣溶膠類型的反演結(jié)果與AERONET地基觀測結(jié)果對比如表3。

表3 不同氣溶膠類型反演結(jié)果與AERONET地基數(shù)據(jù)對比Table 3 Comparison of inversion results of different aerosol types with AERONET ground-based data

相比之下，三種類型的氣溶膠在不同時相的AOD反演結(jié)果基本一致，除個別天數(shù)的AERONET地基數(shù)據(jù)缺失之外(如8月25日和9月9日的AERONET數(shù)據(jù)缺失)，各氣溶膠類型下的AOD的反演結(jié)果均略高于測量值。分析各氣溶膠類型反演結(jié)果與AERONET觀測值的絕對誤差時，發(fā)現(xiàn)9月4日影像的反演序列中的誤差比較大，在0.26～0.3之間，分析可知，當日北京與北京周邊有大量薄云覆蓋，本算法受此影響，將薄云判斷成了氣溶膠；另外，10月4日中三者的絕對誤差較為分散，原因在于C1型和大陸型對氣溶膠變化的不敏感，相比而言C4型氣溶膠能做到對氣溶膠粒子的有效識別。

為進一步討論AOD反演結(jié)果與地基觀測數(shù)據(jù)的關(guān)系，分別將C1、C4和大陸型氣溶膠的反演結(jié)果與AERONET地基觀測數(shù)據(jù)構(gòu)建了線性方程，相關(guān)系數(shù)分別為0.510，0.890，0.836。其中C4型的反演結(jié)果與地基觀測結(jié)果更接近，相關(guān)性為0.890，其次為大陸型氣溶膠，相關(guān)性為0.836。由此可知，聚類結(jié)果中以細粒子為主導(dǎo)的C4型氣溶膠較模型內(nèi)置的大陸型氣溶膠更符合研究區(qū)域的氣溶膠類型，其AOD反演結(jié)果也更接近實測值。

3 結(jié) 論

基于GF-4衛(wèi)星多光譜影像對京津冀區(qū)域的氣溶膠光學厚度的反演，通過調(diào)整暗目標法使之適應(yīng)于GF-4多光譜影像數(shù)據(jù)；分析模型輸入?yún)?shù)中氣溶膠類型對AOD反演結(jié)果的影響，使用PSO聚類算法對研究區(qū)域氣溶膠類型進行聚類分析。使用聚類產(chǎn)生的氣溶膠類型構(gòu)建查找表對研究區(qū)域AOD進行高精度反演，并將反演結(jié)果分別與MOD04_L2產(chǎn)品和AERONET地基觀測結(jié)果進行比較，結(jié)果顯示C4型氣溶膠類型對AOD具有很強的敏感性，與研究區(qū)域的氣溶膠特點更接近，反演結(jié)果與地基數(shù)據(jù)相關(guān)性為0.890；將反演結(jié)果與MOD04_L2產(chǎn)品進行一致性分析，相關(guān)性為0.794。結(jié)果表明以細粒子和吸濕性粒子為主的C4型氣溶膠類型滿足該地區(qū)霧霾時相下的氣溶膠狀況，能夠有效減小在AOD反演中氣溶膠類型引起的誤差，同時也證明了該聚類算法的有效性和可靠性。