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

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

引 言

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

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

針對這些問題，提出了一種新方法：一方面調整暗目標法以適用于高分四號(GF-4)多光譜遙感數(shù)據(jù)，構建類三元一次方程組；另一方面利用粒子群(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反演結果進行對比，并進行分析和討論。

1 實驗部分

1.1 算法流程及暗目標法的調整

算法流程如圖1所示，首先針對暗目標法不適用于GF-4數(shù)據(jù)的情況，對GF-4/PMI傳感器的數(shù)據(jù)進行了分析，并根據(jù)已有數(shù)據(jù)和資料對暗目標法進行調整。利用長時序的AERONET地基數(shù)據(jù)和PSO聚類算法對京津冀地區(qū)的氣溶膠類型進行了聚類分析，并結合實驗驗證選取了符合實際情況的氣溶膠類型。將聚類得到的氣溶膠類型代入6S輻射傳輸模型建立查找表，結合調整之后的暗目標法建立類三元一次方程組。最后選取京津冀區(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

假設地球為均勻郎伯表面，則星載傳感器接收到的大氣頂層反射率(即表觀反射率)ρ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)暗目標法利用短波紅外波段與紅藍波段之間地表反射率之間的相關關系。由表1可知GF-4衛(wèi)星PMI傳感器缺少暗目標法所需的短波紅外波段，因此需要對暗目標法進行調整和改進。假設GF-4衛(wèi)星PMI多光譜數(shù)據(jù)中暗像元的紅藍波段地表反射率間具有良好的相關性，即

(3)

式(3)中，k為相關系數(shù)；b為相應偏差。

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

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

(4)

1.2 粒子群聚類算法

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

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

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

(5)

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

(6)

(7)

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

(8)

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

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

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

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

根據(jù)長時序持續(xù)觀測的AERONET地基數(shù)據(jù)發(fā)現(xiàn)，使用氣溶膠微物理特性和光化學特性參數(shù)共同表征氣溶膠類型，并利用聚類算法得出氣溶膠類型是可行的[11-12]。以京津冀地區(qū)為研究區(qū)域，選擇AERONET地基網(wǎng)絡中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ù)量是聚類分析的關鍵參數(shù)之一。聚類數(shù)量m不斷增加時，使用PSO聚類算法充分迭代求得各樣本到自身所對應聚類中心的歐氏距離之和的最小值，用能量E表示(即最優(yōu)適應度fitness)，E與m的關系為E=f(m)，根據(jù)m對于E的敏感性和擬合關系進行分析，可得到最佳聚類數(shù)量m。圖4反映了兩者在實數(shù)平面坐標上的散點分布情況，將其放到對數(shù)坐標軸上(圖4右上方)，可知能量E和m直線擬合度較好，證明了在實數(shù)坐標空間內E隨m變化而指數(shù)衰減。由圖4可以看出，聚類數(shù)量m足夠大時，新增聚類所占比例越來越小，能量E的衰減程度也越來越弱, 表明新增聚類對于氣溶膠類型差異性的影響越來越小。

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

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

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

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

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

圖6 5種氣溶膠模式復折射指數(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型氣溶膠類型的反演結果相似，且二者的物理和光學特性非常相似，為此不分析C2，以C1，C4以及大陸型氣溶膠類型進行AOD的反演。

2 結果與討論

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

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

圖7 GF-4衛(wèi)星氣溶膠光學厚度反演結果及一致性分析(a)：原始數(shù)據(jù)；(b)：C1型反演結果；(c)：C4型反演結果；(d)：大陸型反演結果；(e)：MODIS溶膠產(chǎn)品；(f)：C1型反演結果分析；(g)：C4型反演結果分析；(h)：大陸型反演結果分析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ù)衰減的關系，可計算得到550 nm波長處的AOD，并選取GF-4衛(wèi)星過境時間前后30min的插值結果求取平均值，當作AERONET的AOD觀測數(shù)據(jù)。另外，以影像中的地基站點為中心，周圍5 km×5 km像元窗口下的AOD反演結果求取平均值，作為GF-4衛(wèi)星多光譜影像的AOD反演結果。京津冀區(qū)域2018年不同時間不同氣溶膠類型的反演結果與AERONET地基觀測結果對比如表3。

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

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

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

3 結 論

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