Shapley-模糊神經網絡方法在華南臺風衛(wèi)星云圖的長時效滾動預測中的應用*

黃小燕 何 立 趙華生 黃 穎 吳玉霜

1.廣西壯族自治區(qū)氣象科學研究所，南寧，530022

2.廣西壯族自治區(qū)氣象臺，南寧，530022

1 引 言

衛(wèi)星云圖以其時空分辨率高、覆蓋面廣的特點在氣象領域得到了廣泛應用，成為除地面氣象觀測資料、空中風資料、空中溫、壓、濕資料等常規(guī)資料以外的重要的非常規(guī)氣象資料。當前衛(wèi)星云圖的應用研究主要是在實時的監(jiān)測分析和目標識別等方面（Welch，et al，1989；張韌等，2004；王敏等，2010；張志清等，2017；李博等，2018；岳治國等，2018；鄭建宇等，2018；丁璐等，2018；郝曉靜等，2019；王新等，2020；崔林麗等，2020；Zhang，et al，2020），人工判讀仍是衛(wèi)星云圖分析的主要方法之一，這既包含一定程度的主觀因素，不利于衛(wèi)星云圖豐富信息的充分提取，同時也有礙天氣預報制作的自動化與定量化。如果能實現(xiàn)對云圖移動和變化狀況進行一定時效的預測，則可從一定程度上彌補實時衛(wèi)星云圖應用在時間上的缺陷，能夠顯著增強衛(wèi)星云圖資料在天氣預報中的實用性和及時性，使衛(wèi)星云圖在災害天氣的預警、預報工作中發(fā)揮更大的作用。

中外對于云圖的預測研究工作雖然已經進行了很久，也取得了較大的研究進展，但是主要以系統(tǒng)保持穩(wěn)定為前提條件，采用一些線性方法來基于云圖局部特征匹配和前后時次云圖運動矢量關系進行線性外推，在對云的發(fā)展趨勢預測和預報時效上存在一定的不足，限制了這些短時效云圖預測的可用度（Endlich，et al，1971；Smith，1975；Arking，et al，1978；Li，1998；Genkova，et al，1999；王登炎，2000；蘭紅平等，2009；蔡叔梅等，2011；梁立為等，2015）。而且，到目前為止，對云圖未來發(fā)展演變的預測研究大多是基于一些短時效（1—3 h）的預測，超過3 h的更長時效的未來云圖狀況變化的云圖預測理論和方法目前并不完善，模型預測效果也不理想 （Hamill，et al，1993；白 潔等 ，1997；王 雷 等 ，1999；龔克等，2000；陳剛毅等，2005；黃勇等，2005；劉科峰等，2006，2008；王偉等，2014）。

由于云的形成是一個非常復雜的過程，云圖隨時間的變化受到云內部和外部大氣環(huán)境的綜合影響，因而云圖與大氣各要素場存在非線性關系，而且云總是在不斷發(fā)生移動、形變、擴展收縮甚至分裂、融合以及更復雜的生、消過程中，使對云圖追蹤過程比計算機視覺中對一般的線性或近線性物體的追蹤問題要難得多。因此，傳統(tǒng)的線性預測方法存在一定的局限性，特別是超過3 h的未來云圖狀況變化，僅依靠線性外推或不涉及云圖移動變化的物理過程、動力和熱力因子，可能難以獲得好的預測結果。一些學者嘗試運用非線性智能計算預報方法來解決這一問題，王繼光等（2007）引入遺傳算法構造了云圖非線性動力預測模型，預報時效達到了42 h，云圖預測模型能夠較為客觀合理地描述特定季節(jié)和區(qū)域云演變的基本結構，主要特征也描述較好。劉科峰等（2008）則提出了奇異值分解和徑向基網絡相結合的云圖預測思想，對云圖特征值進行動力模型重構和模型參數(shù)反演的云圖演變非線性動力預測方法，該方法能合理地描述云運動的基本特征和演變趨勢。但上述方法都沒有考慮大氣中水汽相態(tài)變化和云物理過程，制約了對云發(fā)展變化本質的刻畫和描述。何如等（2010）、金龍等（2011）利用數(shù)值預報模式的物理量預報產品作為各預報分量的預報因子，并采用系統(tǒng)降維計算處理方法分別建立相應的時間系數(shù)神經網絡預報模型，得出未來時刻的衛(wèi)星云圖預報圖。但是這種非線性衛(wèi)星云圖預報模型在因子處理方法和預報模型上都還是處于探索階段，云圖的預測準確率也有待進一步的提高。因此，文中嘗試根據云團發(fā)生、發(fā)展直至消亡的整個生命過程會明顯受到大氣環(huán)境場各個物理量因子影響的特點，并且云態(tài)變化更多的是非平穩(wěn)、非線性的不規(guī)則變化，從云的演變過程與環(huán)境場影響因子的相關出發(fā)，建立二者的非線性統(tǒng)計關系模型，進一步設計采用k-近鄰互信息估計的分步式變量選擇算法（Kraskov，et al，2004；韓敏等，2012），通過兩步過程分別實現(xiàn)相關因子的選擇與弱變量的剔除，并結合運用合作對策的Shapley值（Shapley，1953；雷勛平等，2012）與模糊神經網絡個體生成相結合的方法，設計一種新的非線性統(tǒng)計集合預報方法，來提高單一統(tǒng)計預報方程的預報技巧和預報穩(wěn)定性，對衛(wèi)星云圖的未來發(fā)展演變進行預測，探索衛(wèi)星云圖較長時效的預報方法。

2 資料與預報建模設計

2.1 資 料

華南區(qū)域受臺風影響較為頻繁，因此文中在建立衛(wèi)星云圖非線性集合預報模型時云圖選用的是2013—2016 年華南區(qū)域（15°—30°N，100°—120°E）有臺風影響過程的資料。

衛(wèi)星云圖資料選用的是FY-2E的云頂亮溫（TBB）資料，云頂亮溫是氣象衛(wèi)星獲取的紅外探測通道的數(shù)據，它是生成紅外云圖以及其他顯示云圖中最原始的定量資料。云頂亮溫是以相對于黑體溫度而言，它的值一般低于0℃，黑體溫度越低表明云頂越高、對流越活躍，降水的潛在性越大。因此，可以根據云頂亮溫以及一系列的相關數(shù)據，通過網絡模型學習方法，為強降水的預測提供可能。文中云頂亮溫投影方式為等經緯度投影；云圖大小為401×301 像素；星下點分辨率為 0.05°×0.05°（5 km）；采樣間隔為 6 h。

與華南臺風衛(wèi)星云圖對應的預報因子來自于數(shù)值模式預報產品，選用歐洲氣象資料中心（ERAInterim）的全球再分析資料（https://www.ecmwf.int/），空間分辨率為 0.75°×0.75°經緯度，時間分辨率6 h（00、06、12 和18時，世界時，下同），200、500、700、850 hPa 4個高度的溫度、風場、位勢高度場、速度場、渦度、散度、水汽通量、水汽平流、相對濕度、濕度以及各種物理量指數(shù)等，共4年（2013—2016年）。所選取的數(shù)值預報模式產品資料范圍與云圖資料范圍一致。在衛(wèi)星云圖的物理量預報因子初選時，針對提取出的每張云圖的n個反映云圖主要特征的預報分量，分別與每個要素場做場相關計算普查，通過設定適宜的閾值得到云圖各個預報分量的初選預報因子。

2.2 預報建模設計

在華南區(qū)域臺風衛(wèi)星云圖6、12、18、24、48和72 h共6個時效的滾動預報建模時，主要的設計思想包括構造一系列云圖的時間序列，根據分辨率的不同每一幅云圖包含了若干的像素點，若將每個點作為研究對象建立預測模型，計算量將非常龐大，同時也會增加模型的噪聲干擾和計算誤差。為此，根據Takens（1981）的相空間重構理論把原始的云圖觀測資料看成是某動力系統(tǒng)的離散值，求解與之相反的數(shù)值問題，即可重構出表現(xiàn)該云區(qū)移動的動力模型。即從云圖資料序列中反演重構描述云區(qū)移動的動力預報模型。云圖資料序列視為三維時間-空間場，研究引入經驗正交函數(shù)分解（EOF）方法，對云頂亮溫的數(shù)值矩陣采用EOF展開方法進行時、空分離，將提取出來的空間場典型模態(tài)對應的時間系數(shù)序列作為模型輸出的預報分量。進一步根據云圖建模樣本作EOF展開后的主分量的方差和累積方差貢獻選出能夠較好地反映云圖主要特征的前若干個主分量對應的時間系數(shù)作為n個預報分量，分別建立非線性集合預報模型。取這些空間特征向量對應的時間系數(shù)作為建立云圖預測模型的預報分量，可以減少噪聲影響和誤差積累，提高預報的精度。

云系的移動和生消變化等與云系周圍大氣環(huán)境條件密不可分，目前以歐洲預報中心為代表的數(shù)值預報模式產品的精度已經得到公認。在此合理地提取反映云系變化的主要特征，并找出這些主要特征與其未來變化密切相關的物理量預報因子，構建合適的數(shù)學預報建模方法，就可能建立較長時效的衛(wèi)星云圖預測模型。文中以歐洲預報中心數(shù)值預報產品為基礎，構造對應的華南臺風衛(wèi)星云圖預報因子。進一步采用k-近鄰互信息估計的分步式變量選擇算法，通過兩步過程分別實現(xiàn)相關變量的選擇與弱相關變量的剔除，從高維預報因子數(shù)據集中提取少數(shù)幾個包含主要信息的變量作為衛(wèi)星云圖變化智能計算方法集成預報模型的輸入，從而構建網絡結構規(guī)模小、有效預報信息量大的非線性集成預報模型，提高衛(wèi)星云圖客觀預報的精度。

目前的神經網絡大多用于回歸估計，集成的輸出通常由各網絡的輸出通過簡單平均或加權平均產生。當神經網絡個體較多，集成規(guī)模較大時，優(yōu)化權值會適得其反，適宜采用簡單平均方法；而在集成規(guī)模較小，或者數(shù)據集中噪音較多時，優(yōu)化權值將會提高學習系統(tǒng)的泛化能力（Sollich，et al，1996）。Jin等（2008）也通過密度函數(shù)的積分推導，證明了當神經網絡集成個體之間的差異越大時，神經網絡集成預報的精度將越高。另外，由于云系移動或演變過程主要是受到云團內部和云系周圍大氣環(huán)境條件的綜合影響，因此其變化特征主要為非平穩(wěn)的非線性變化。為此，華南區(qū)域臺風衛(wèi)星云圖的預測模型主要采用合作對策的Shapley值與模糊神經網絡相結合的方法，設計一種新的非線性統(tǒng)計集合預報方法，確定模糊神經網絡集成預報個體的權重系數(shù)，以此提高單一統(tǒng)計預報方程的預報技巧和預報穩(wěn)定性。

最后，計算得到華南區(qū)域臺風衛(wèi)星云圖獨立樣本序列的時間系數(shù)的預測值。與該云圖樣本通過EOF得到的空間特征向量經過時、空反演重構，得到云圖獨立樣本二維云頂亮溫矩陣的預測值，并對結果進行分析檢驗。圖1為基于Shapley-模糊神經網絡的華南區(qū)域臺風衛(wèi)星云圖建模設計流程。

圖1 基于 Shapley-模糊神經網絡的華南區(qū)域臺風衛(wèi)星云圖建模設計流程Fig.1 Flow chart of the forecast model based on the Shapley-fuzzy neural network for satellite cloud images of typhoons in South China

3 方 法

3.1 k-近鄰互信息估計的分步式變量選擇算法

針對每張云圖都存在的高維因子集問題，采用k-近鄰互信息估計的分步式變量選擇算法，通過兩步過程分別實現(xiàn)相關變量的選擇與弱相關變量的剔除，替代原有的多元線性回歸技術。該算法采用Kraskov等（2004）提出的基于k-近鄰互信息計算方法，避免了直接對變量進行概率密度估計，適用于高維變量的互信息計算。其基本思想是：在X和Y構成的空間Z=(X,Y)中，以 εi/2為點zi=(xi,yi)到其k-近鄰的距離（通過最大值范數(shù)來度量）， εx(i)/2為點zi=(xi,yi)到x軸上相應點的距離，同理可得到εy(i)/2。統(tǒng)計到點xi的距 離 嚴 格 小 于 εi/2的 數(shù) 目（nx(i)），同樣對變量Y作相同的處理得到ny(i)。則變量X、Y的互信息通過下式計算

式中，N為隨機變量X的樣本數(shù)， εi(i=1,2,···,N)為任意小的數(shù)（文中取0.01）， ψ (?)為雙Γ函數(shù)，且滿足ψ(1)=?0.5772516， ψ (?+1)= ψ(?)+1/?；符號 〈? 〉表示對其中的所有變量i=1,2,···,N取平均；同時，近鄰數(shù)目k的取值決定估計密度函數(shù)的光滑性，k值越大，函數(shù)曲線越光滑，同時互信息量I?的統(tǒng)計誤差越小，而系統(tǒng)誤差卻相反。為了平衡這兩種誤差，經過多次計算試驗，文中k=8。

分布式算法的變量選擇過程分為相關變量的選擇與弱相關變量的剔除兩個步驟，理論上預報模型的輸入矩陣只包含強相關變量，而不應包含無關變量和弱相關變量。因此，分步算法首先考慮輸入變量與輸出變量的相關程度，分別得到強相關變量和弱相關變量的子集，由于弱相關變量中包含冗余變量，需要進一步對子集進行篩選，剔除冗余變量。如果輸入變量Xi滿足

式中， δ1(δ1∈ [0,1])為相關系數(shù)閾值，則表明Xi包含了關于Y的一定量的信息，即Xi為Y的相關變量。如果輸入變量Xi不 滿足式（2），可將Xi視為無關變量，將其從輸入子集中剔除。一般而言，當特征維數(shù)較高時， δ1可以設置得較大，從而有效降低特征維數(shù)。因為文中試驗的預報因子特征維數(shù)較高，經綜合考慮和試驗，δ1=0.35較合適。

通過k-近鄰互信息估計的分步式變量選擇算法從大量云圖備選因子中選擇出有限的預報因子后，可以充分挖掘預報因子集的有效預報信息，增強預報方法的物理基礎，提高預報模型的預報精度。

3.2 Shapley-模糊神經網絡集合預報個體生成方法和原理

3.2.1 模糊神經網絡集合預報個體生成

考慮到云變化的復雜性和非線性，在建立云圖未來變化的預測模型時采用了模糊神經網絡計算模型作為基礎模型。模糊神經網絡本質上也是一種多層前饋網絡，文中所采用的模糊神經網絡結構為數(shù)值型4層前饋網絡，包括輸入層、隸屬度生成層、推理層和反模糊化輸出層（王士同，1998），其中模糊規(guī)則形式為：

文中在模糊神經網絡的隸屬度生成層中，采用比較常見的高斯函數(shù)作為隸屬度函數(shù)，該函數(shù)具有較高的擬合和隸屬性，比較適合數(shù)據量比較大、預報精度要求較高的情況下運用。高斯隸屬度函數(shù)僅由2個參數(shù)確定，公式表達為

式中，m是模糊分割數(shù)，aij和 σij分別確定隸屬函數(shù)的中心和寬度。

在模糊邏輯中，按照“和”運算，一般有3種計算方法：①按照取最小方式，②按照相乘方式，③按照相加方式。研究結果（Zeng，et al，1996）表明，相乘方式的性能優(yōu)于取最小方式。文中在推理層中選擇第②種方法，即各節(jié)點的輸出值分別為該節(jié)點所有輸入的代數(shù)乘積。

在輸出層中，采用如下的反模糊化網絡輸出。

式中， ωj（j= 1,···,m）為連接權。

模糊神經網絡是利用前饋網絡的BP算法來訓練調整參數(shù)，在網絡訓練過程中，設隸屬度函數(shù)為式（3），首先利用構建的預報因子學習矩陣對網絡進行訓練，其主要計算步驟歸結為

（1）設初始值i＝1，學習率 α=0.9，網絡訓練誤差為e1。

（2）初始時刻，對網絡的連接權、隸屬函數(shù)的中心值及寬度值均用隨機數(shù)進行初始化。

（3）采用BP算法對參數(shù)和權值進行訓練調整。

（4）計算網絡的實際輸出與期望輸出的誤差e2。

（5）當e2＞e1時，返回到步驟（3），否則網絡訓練結束，利用訓練得到的網絡參數(shù)和連接權進行預測，得到預測個體。

3.2.2 合作對策 Shapley 值在個體中的分配

多個模糊神經網絡預測個體生成后會存在相互間的過擬合和泛化性問題，為了改進模糊神經網絡個體預報能力，在進行華南區(qū)域臺風衛(wèi)星云圖變化較長時效預測模型研究時，主要設計采用Shapley值方法，以組合預測有效度的平方和作為合作準則，按合作對策的Shapley值在各模糊神經網絡預報個體中進行分配，從而確定各預報個體在模糊神經網絡預報中的權重系數(shù)，提高模糊神經網絡預報個體的差異度，建立一種新的衛(wèi)星云圖變化的機器學習客觀集成預報模型。Shapley值為如何決策一個n人討價還價博弈中每個參與者所得的分配比例提供了一種很好的方法。近年來，學者們根據實際需要，對該方法做了相應改進（陳菁等，2011；陳啟明等，2012）。陳啟明等（2012）將各單項預測模型看作組合預測合作對策的局中人，按合作對策的Shapley值確定各單項預測模型在組合預測中的權重系數(shù)。

具體假設存在m個模糊神經網絡對某一個云圖主分量Yt(t=1,2,···,N)進行了預測，其預測結果記為M={yi,i=1,2,···,m}，則M為組合預測方法的局中人集合。

設M的所有子集為2M，則M中的任一子集s∈2M形成組合預測方法的一個聯(lián)盟，若干個局中人結成聯(lián)盟后，這個聯(lián)盟作為一個整體進行組合預測就是希望盡可能提高預測有效性。

設M={1,2,···,m}，s?M，v(s)為定義在 2M集合上的 實 值 函 數(shù) ，令v(s)=F(s)，滿 足v(Φ)=0，v(M)≥則稱組合預測方法為一合作m人對策。記為 Γ =[M,v]。v(s)稱為m人對策的特征函數(shù)，其中，F(xiàn)(s)表示聯(lián)盟s進行組合預測所得的預測有效度的平方和。

稱v(s∪{i})?v(s)為第i種各單項預測方法對聯(lián)盟s合作的“貢獻”，其中s?M。記 φi(v)為第i種各單項預測方法的平均“貢獻”，有

式中，s為M中包含 {i}的所有子集合； |s|為子集s中局中人的單項預測模型的個數(shù)。稱φ(v)=(φ1(v),φ2(v))為合作m個模糊神經網絡個體對策 Γ =[M,v]的Shapley值，可以證明

式中，各單項模糊神經網絡預報個體在組合預測方法中的平均貢獻 φi(v)之和等于合作的總成果，各單項模糊神經網絡預報個體在組合預報方法中的加權系數(shù)的大小應根據合作中的平均貢獻來確定。考慮到v(M)為組合預測有效度，將 φi(v)作如下歸一化處理，可求得各單項模糊神經網絡預報個體的權重系數(shù)

4 華南臺風衛(wèi)星云圖的非線性智能預報試驗

4.1 預報建模樣本和試驗預報樣本

臺風所引發(fā)的強降水是目前大氣科學研究的重點和難點，若能較為準確地預報未來云系的發(fā)生、發(fā)展情況，就能較好地判斷臺風強降水過程是否發(fā)生。因而，在建立衛(wèi)星云圖非線性集合預報模型時，主要選擇的是華南區(qū)域有臺風影響過程的云圖資料作為預報試驗對象。根據上述云圖資料的選擇標準以及云圖資料的完整性，主要選擇2013—2015年的臺風過程6 h間隔共196張云圖資料作為初次的預報建模樣本建立預報模型，后續(xù)隨著預報樣本的增加，建模樣本也相應增加。為了提高預報建模研究的客觀性，以相同的標準，選擇2016年的臺風影響過程，共86張云圖作為預報模型的獨立樣本進行實際的滾動預報檢驗。

4.2 預報對象的處理

由于每一張衛(wèi)星云圖的像素點較多（120701個），如果對全部每個點都進行預報建模，計算量巨大，耗費時間長，對實際業(yè)務預報不適用。因此，采用自然正交函數(shù)分解方法對云圖的實際云頂亮溫數(shù)據場時間序列作時-空正交分解，將時、空要素場轉化為若干空間的基本模態(tài)和相應的時間系數(shù)序列的線性組合，對時間系數(shù)進行建模預測，對預測后的時間系數(shù)與空間模態(tài)進行重構，得到新預測的衛(wèi)星云圖。

首先，采用EOF方法，對建模云圖樣本作分解

即任一時刻的云圖可以表示為

式中，v(x,y)為EOF的主分量，ξi(t)為相應于各主分量的展開時間系數(shù)，i為自然正交展開的階數(shù)。表1給出了云圖建模樣本作EOF展開后，前30個主分量展開的方差和累積方差貢獻。可以看到，前30個主分量累積方差貢獻接近75%，與原始云圖的平均相關率為81.3%。進一步增加10個主分量（第31—40個）進行計算分析可以發(fā)現(xiàn)，其累積方差貢獻僅提高4.07個百分點，對應的相關率基本保持不變；再繼續(xù)增加到50個主分量，累積方差貢獻也僅增加了3.09個百分點，相關率基本沒有提高。分析可知，后續(xù)的方差貢獻率增加越來越緩慢，但是增加的主分量越多，意味著反演結果的誤差也會增大。因此，綜合比較分析，前30個主分量已經能夠較好地反映云圖的主要特征，故取前30個主分量進行預報試驗。

這30個主分量能否較好地還原出原始衛(wèi)星云圖的主要畫面特征呢？進一步采用分解出的這30個主成分的時間系數(shù)與空間相應的空間向量進行重構合成繪圖，并與原始衛(wèi)星云圖進行比較。運用ENVI軟件進行衛(wèi)星云圖的繪制，為了進行客觀的對比分析，在云圖的顏色設定過程中，以每一張原始圖片的最大、最小值作為基準，從2.1節(jié)介紹的云頂亮溫的設定可知，溫度越低，表明云頂越高，對流越活躍，從而降水的潛在性也越大，因而溫度低即值越小的數(shù)據用白色表示，最大值設定為黑色，從白色到黑色一共分為256級，如此設定得到每一張云圖的標準色標。圖2和3分別為2013年1306號臺風“溫比亞”和2014年1409號超強臺風“威馬遜”的4個時次的原始云圖與提取的30個主成分向量合成云圖的對比?？梢园l(fā)現(xiàn)，前30個主分量描繪出的云圖輪廓與原始云圖的輪廓基本相似，能較好地反映原始云圖的主要特征。為此，以這前30個主分量對應的時間系數(shù)T1,T2,T3,···,T30作為30個預報分量，分別建立華南區(qū)域臺風衛(wèi)星云圖的非線性集合預報模型。再根據式（10）將預報出的時間系數(shù)與相應的空間向量合成，即可得出未來時刻的預報云圖。

4.3 預報因子的處理

由云圖觀測資料可以清楚地看到，大量的云系移動、生消變化都非常快，特別是一些強對流云系，可能在較短的時間內發(fā)生顯著變化。這種云系隨時間的快速變化，其根本原因是云系內部的動力、熱力、水汽狀況發(fā)生了變化，同時也與云系周圍大氣環(huán)境條件密不可分。因此，要想能夠有效地預測云系在相對長時間間隔（如24 h后）的狀況，必須在預報模型中考慮上述的云系內部、外部大氣物理因子狀況，建立云系的移動、變化狀況與導致云系變化的大氣物理量因子的非線性關系。

按上述方法普查可以得到的初選云圖預報因子數(shù)量眾多，如果全部將這些預報因子放入模型中進行建模，則預報模型的輸入結構將過于龐大，容易出現(xiàn)過學習的現(xiàn)象。為了減少預報模型的輸入節(jié)點，同時又盡可能保留所有預報因子攜帶的預報信息，在此采用k-近鄰互信息估計的分步式變量選擇算法，通過兩步過程分別實現(xiàn)相關變量的選擇與弱相關變量的剔除，對預報因子群進行降維去噪處理。具體做法為

表1 云圖前 30 個主分量展開的方差和累積方差貢獻 （%）Table 1 Variances and cumulative variance contributions of the first 30 principal components of the satellite cloud images （%）

（1）采用式（1）對初選的云圖預報因子群進行互信息計算。

（2）通過式（2）判定該預報因子與預報對象（云圖時間系數(shù)）是否為相關變量，滿足則選入，無關變量則剔除。

圖2 臺風“溫比亞”4 個時次 （2013 年 7 月 1 日 12 （a、e）、18 （b、f） 時，7 月 2 日 00 （c、g）、06 （d、h） 時 ）實況云圖 （a—d）與 30 個主成分重構云圖 （e—h）Fig.2 Four times （12:00 （a, e），18:00 （b, f） UTC 1 July 2013，00:00 （c, g），06:00 （d, h） UTC 2 July 2013） of observational cloud images （a－d） of Typhoon "Rumbia" and predicted cloud images （e－h(huán)） reconstructed by 30 principal components

續(xù)圖2Fig.2 Continued

圖3 超強臺風“威馬遜”4 個時次 （2014 年 7 月 18 日 12 （a、e）、18 （b、f）時、19 日 00 （c、g）、06 （d、h）時 ）實況云圖 （a—d）與 30 個主成分重構云圖 （e—h）Fig.3 Four times （12:00 （a, e），18:00 （b, f） UTC 18 July 2014，00:00 （c, g），06:00 （d, h） UTC 19 July 2014） of observational cloud images （a—d） of super typhoon "Rammasum" and predicted cloud images （e—h） reconstructed by 30 principal components

續(xù)圖3Fig.3 Continued

（3）對所有云圖預報因子做k-近鄰互信息估計的分步式變量選擇算法篩選后得到的所有因子組合在一起構成新的云圖預報因子集，新組成的高相關云圖預報因子群進一步進行顯著性檢驗，保留通過顯著性檢驗的預報因子。

具體以6 h的預報時效為例，在對6 h的30個云圖時間系數(shù)預報分量的智能計算預報模型研究時，首先根據歐洲預報中心數(shù)值預報模式6 h間隔的各種物理量相關要素場普查初選的云圖預報因子。針對第1個云圖時間系數(shù)，首先進行歐洲預報中心模式中的某一要素某一層次的每一個格點（共567個）分別與云圖時間系數(shù)的相關計算，將成片（≥10個格點）穩(wěn)定高相關（相關系數(shù)≥0.3）格點作為預報因子的選擇區(qū)，在區(qū)內選3個相鄰格點的最大平均值作為該相關區(qū)的代表值，作為待選因子。共在21個云圖物理量場中提取預報因子，包括各層次的散度場、垂直速度場、渦度場、比濕場等，以及各種云圖與水汽相變相關的物理量指數(shù)—K指數(shù)場、ky指數(shù)場，SI指數(shù)場等。經過初選因子的相關計算后，第1個云圖時間系數(shù)預報分量的預報因子的相關系數(shù)最大(0.7042)，入選的因子總數(shù)量超過2000個?？梢园l(fā)現(xiàn)，隨著云圖時間系數(shù)預報分量的增加，相關在不斷下降，到第30個時，初選因子的最大相關系數(shù)僅為0.4033，相關程度下降近一半。分析可知，前3個云圖時間系數(shù)預報分量對于云圖的預報精度有更大貢獻。

進一步用k-近鄰互信息估計的分步式變量選擇算法、設定 δ1=0.35閾值計算判識預報對象（云圖時間系數(shù)）的高相關變量，此時第1個主分量的預報因子已經只剩下247個，而第30個主分量的數(shù)量僅為64個。進一步進行顯著性檢驗后，第1主分量選入預報模型僅36個，而第30個主分量的因子數(shù)量為12個。

值得注意的是，進行顯著性檢驗后的預報因子集可以很好地控制因子間的復共線性關系，同時能分解出有物理意義的云圖預報因子，剔除影響預報效果的噪音因子，從而控制選入最終預報模型的預報因子數(shù)量。一般來說，適宜的預報因子矩陣結構進入預報模型能得到較高的預報精度，但是模型輸入過于龐大，則容易產生過擬合等問題（金龍等，2004）。

4.4 Shapley-模糊神經網絡華南臺風衛(wèi)星云圖的智能計算滾動預報試驗

利用3.2節(jié)Shapley-模糊神經網絡集合預報建模的集合預報個體生成方法，分別建立預報時效為6、12、18、24、48和 72 h的各 30個云圖時間系數(shù)預報分量的智能計算集合預報模型。在此，模糊神經網絡預測個體生成的參數(shù)設置統(tǒng)一為：各云圖預報分量的預報因子數(shù)量為網絡的輸入節(jié)點，輸出節(jié)點為1，3個推理層節(jié)，網絡訓練次數(shù)設定為1000次，學習因子取0.9，總體誤差定為0.001。訓練過程中，為了提高計算效率，設定了兩種終止訓練的條件，一是計算誤差小于目標誤差（0.001），二是訓練次數(shù)達到1000次，計算過程滿足其中的條件之一，就終止網絡訓練。

具體以6 h預報時效為例，針對衛(wèi)星云圖的第一個時間系數(shù)，采用Shapley值方法，以組合預測有效度的平方和作為合作準則，按合作對策的Shapley值對各模糊神經網絡預報個體進行合理分配，確定各預報個體在集成預報中的權重系數(shù)，提高集成預報個體的差異度，并充分考慮各單項預報模型對提高預測云圖有效度的貢獻，從而建立云圖第一個時間系數(shù)的預報模型進行云圖預測建模。在計算過程中，按照時間順序，逐次滾動進行，即用建模樣本（196個）確立預報模型后，采用該預報模型首先對第一個獨立預報樣本進行預報，得到預測結果后，加入到建模樣本中（此時建模樣本增加1，為197個），接著對第二個獨立樣本進行預報檢驗，依此類推，完成最后一個（第282個）獨立預報樣本的預報，如此就得到云圖第一個時間系數(shù)的全部86個獨立預報樣本的預測值。按該方法和預報步驟，對其余的29個云圖時間系數(shù)分量進行預報，得到最終的30個云圖時間系數(shù)的全部獨立預報樣本的預測值，使用該預測結果，與相應的30個云圖空間向量進行合成，可以得到預報樣本的每一張云圖的預報結果。按照此方法，完成所有預報時效的云圖滾動預測。

4.5 云圖預測效果檢驗分析

4.5.1 相關性評估分析

通過上述的建模預報試驗，可以分別得到預報時效為 6、12、18、24、48和 72 h的各 30個時間系數(shù)預報分量的Shapley-模糊神經網絡集合預報結果，分別與對應的空間系數(shù)進行合成，得到預測云圖。檢驗預報效果可用每一張預測云圖與實況云圖的像素點相關進行評估。這里進一步分別統(tǒng)計這6個預報時效的每一張預測云圖的像素點與實況云圖像素點的相關系數(shù)（圖4a）?？梢钥吹剑路桨傅?個預報時效的預報中，這些華南區(qū)域臺風影響過程86張預測云圖與實況云圖的相關系數(shù)主要分布在0.4—0.8，相關系數(shù)為0.5—0.6的個數(shù)以及0.6—0.7個數(shù)最多，其次為相關系數(shù)0.7—0.8的個數(shù)，0.4以下相關較差的云圖個數(shù)較少?？傮w而言，新方案預測云圖與實況的相關總體較好，相關系數(shù)0.5以上的個數(shù)超過總樣本的半數(shù)。

圖4 新方案預測云圖與實況云圖云頂亮溫值的相關 （a.相關系數(shù)分布，b.平均相關系數(shù)）Fig.4 Correlation between cloud images predicted by the new scheme and cloud top brightness temperatures of the observed cloud images （a.distribution of correlation coefficient，b.average correlation coefficient）

進一步分別統(tǒng)計6個預報時效的平均相關系數(shù)（圖4b）發(fā)現(xiàn)，隨著預報時效的延長，相關逐漸下降，這與實際情況吻合。主要原因與預報因子隨著預報時效的延長相關逐漸減弱有關。但是值得注意的是，云圖預測時效達到了72 h，這樣長的預報時效是有意義的。

為了客觀分析基于k-近鄰互信息估計的分步式變量選擇的Shapley-模糊神經網絡集合預報方法對華南區(qū)域臺風衛(wèi)星云圖預報的優(yōu)越性，將其預報結果與傳統(tǒng)采用的線性方法進行比較。針對同樣的建模樣本和獨立預報檢驗樣本，采用同樣的方案選出的云圖預報因子進行多元線性回歸的預報建模運算?？梢园l(fā)現(xiàn)（表2），在每一張預測云圖的像素點與實況云圖像素點的相關系數(shù)分布場上，多元線性回歸方案的6個預報時效的預報中，這些臺風過程86張預測云圖與實況云圖的相關系數(shù)主要分布于0.3—0.7，相關為0.4—0.5的個數(shù)最多，其次為0.5—0.6個數(shù)，0.7以上相關較強的云圖個數(shù)較少。對比分析可知，在預測云圖與實況云圖的相關上，新方案的相關性明顯高于多元線性回歸方案。6個預報時效的平均相關系數(shù)的對比情況同樣如此（圖5），兩個方案的預報結果隨著預報時效的增加，相關都是逐步下降的，但是新方案的預報精度均明顯高于多元線性回歸方案。

另外，兩種預報方案都存在預測云圖與原始云圖相關系數(shù)低于0.3的預報較差個例，以預報效果最佳的6 h預報時效進行分析發(fā)現(xiàn)，新方案中這樣的個例有8個，其時序分布基本不連續(xù)。分析這些云圖個例的特征可以發(fā)現(xiàn)，共同表現(xiàn)為云型比較松散，對流云體積小并且分布零散，這些個例較多為臺風前期發(fā)展和后期消散階段，說明預報方案對該時段分散凌亂云的預報能力相對較弱。而線性回歸方案在這方面的表現(xiàn)卻不太一致，其預報較差的個例基本與新方案的個例不同，9個個例的時序連續(xù)分布較多。云型的分布特征差別也較大。后者中大多數(shù)云圖的有云區(qū)和無云區(qū)的分布較為分明，對流云體積較大，云系處于較為旺盛的階段，但是線性回歸方案對這些云系的預報能力卻不足。在預報因子一致的情況下，兩種方案預報表現(xiàn)卻幾乎不一樣，這除了與線性和非線性方案的預報性能有關外，是否還與各方案對入選的預報因素的敏感性有關，值得進一步深入研究。

表2 新方案和多元線性回歸方案預測云圖與實況云圖云頂亮溫值的相關系數(shù)樣本數(shù)對比 （單位：個）Table 2 Comparison of the number of samples for the correlation coefficients of the cloud images predicted by the new scheme and by the multiple linear regression scheme with the cloud top brightness temperatures of observed cloud images （unit：piece）

圖5 不同預報時效新方案及多元回歸預測云圖與實況云圖云頂亮溫值的相關系數(shù)Fig.5 Distribution of correlation coefficient between the cloud images predicted by the new scheme and the cloud top brightness temperatures of the observed cloud images

4.5.2 實際預報個例分析

為了更直觀地對比分析非線性新方案和線性方案預報性能的差異，進一步針對預報檢驗的獨立樣本進行預測云圖與實況云圖的繪制，繪制方案與4.2節(jié)介紹的繪制方案相同。圖6和7為2016年影響華南區(qū)域的臺風“妮妲”實況云圖與預測重構云圖的對比。其中圖6為“妮妲”發(fā)展旺盛期的衛(wèi)星云圖，圖7為其減弱接近消失的云圖。

具體分析可知，圖6a為2016年8月2日18時的實況云圖，圖6b為30個主成分未做預報的擬合云圖，兩者相似度很高，說明提取的30個主成分分量是可以很好地展示原始云圖的主要特征的。圖6c—h分別為預報時效 6、12、18、24、48和72 h的各30個時間系數(shù)預報分量的Shapley-模糊神經網絡集合預測云圖?？梢钥吹剑?個時效預測云圖與實況云圖在云的基本輪廓、紋理特征分布以及云系強、弱方面相似度都很高。特別是在臺風的主體云團即廣西西南部到越南和北部灣一帶的云圖的預報上，預測云圖基本都能有所體現(xiàn)，形狀較好，強度也較接近實況。而多元線性回歸方案所得到的預測云圖（圖6i—n），則明顯在云系強弱上有所欠缺，清晰度也不夠，即黑白不分明，云系輪廓和紋理都較難分辨。

圖6 臺風“妮妲”實況云圖與預測重構云圖 （2016 年 8 月 2 日 18 時實況云圖 （a），30 個主成分重構云圖 （b），非線性方案預測云圖 （c—h 對應 6—72 h）；線性方案預測云圖 （i—n 對應 6—72 h））Fig.6 Observed and predicted cloud images of Typhoon "Nida" at 18：00 UTC 2 August 2016（a.observed cloud images，b.cloud images predicted by 30 principal components，c—h.cloud images predicted by nonlinear scheme （6—72 h），i—n.cloud images predicted by linear scheme （6—72 h））

續(xù)圖6Fig.6 Continued

續(xù)圖6Fig.6 Continued

分析預測的像素點上的數(shù)值可以發(fā)現(xiàn)，實況云圖的最小值為1760（對應白色，云圖發(fā)展旺盛，對流較強），最大值是2940（對應黑色，為無云狀態(tài)）。新方案所預測的云圖的像素值最小值的6個預報時效為1793—1973，最大值則分布在2896—3043，與實況云圖在最小值上的差距為33—213，差距較小。因此，能較好地顯示臺風云圖對流旺盛的區(qū)域。多元線性回歸方案的最大、最小值的預報與實況云圖差距較大，云圖的像素值最小值的6個預報時效在2183—2392，最大值則分布在2814—2955，與實況云圖在最小值上的差距較大，為423—632，因而在強度上預報都偏弱，不能較好地顯示云圖發(fā)展強盛的部分。

圖7 臺風“妮妲”實況云圖與預測重構云圖 （2016 年 8 月 3 日 00 時 ）（a.云圖，b.30 個主成分重構云圖，非線性方案預測云圖 （c—h 對應 6—72 h），線性方案預測云圖 （i—n 對應 6—72 h））Fig.7 Observed and predicted cloud images of Typhoon "Nida" at 00：00 UTC 3 August 2016，observed cloud images （a），cloud images predicted by 30 principal components （b），cloud images predicted by nonlinear scheme （c—h correspond to 6—72 h）， cloud images predicted by linear scheme （i—n correspond to 6—72 h）

續(xù)圖7Fig.7 Continued

續(xù)圖7Fig.7 Continued

分析圖7的情況亦如此，新方案6個時效預測云圖（圖7c—h）與實況云圖在云的基本輪廓、紋理特征分布以及云系強弱方面相似度比多元線性回歸方案（圖7i—n）都更好。這是“妮妲”的一個消亡減弱階段圖，實況圖上顯示其較強的主體云團主要位于越南和北部灣上空，廣西境內還存在一些殘留云系。新方案6個時效預測云圖（圖7c—h）在位置上預報把握較理想，強度則比實況云圖略強。而多元線性回歸方案（圖7i—n）在強度上的預報則同樣表現(xiàn)偏弱，圖片清晰度較差?？傮w上，新方案預測的重構云圖效果優(yōu)于多元線性回歸方法。

5 結 語

為更好地利用大量的衛(wèi)星云圖觀測資料提高對臺風暴雨的預報能力，解決并提高對臺風強降水云系變化的預報精度，加大對未來云系變化的預報時效，根據華南區(qū)域臺風衛(wèi)星云圖生消、演變的非線性、時變性和非平穩(wěn)性特點，從尋找對云團未來變化有顯著影響的大氣物理量預報因子出發(fā)，進行未來6—72 h共6個預報時效的云圖變化趨勢集合預報建模。

該方案的主要特點是：在云圖預報建模過程中，采用經驗正交函數(shù)分解方法進行華南臺風云圖的特征抽取，并以時間系數(shù)作為預報建模的預報分量，建立從云圖中抽取的主要特征向量與數(shù)值預報模式產品的非線性映射關系，使各預報分量預報模型具有很好的物理基礎。并進一步對預報因子采用k-近鄰互信息估計的分步式變量選擇算法，通過兩步過程分別實現(xiàn)相關變量的選擇與弱相關變量的剔除，充分挖掘了預報因子集的有效預報信息。而在數(shù)學預報建模方法上，設計了類似于數(shù)值預報模式集合預報方法的Shapley-模糊神經網絡集合預報方法。該方法不同于大氣科學中傳統(tǒng)的統(tǒng)計預報方法，是一種基于智能計算方法的非線性統(tǒng)計集合預報建模方法。從研究的華南區(qū)域實際臺風衛(wèi)星云圖預測試驗結果看，該方法對各預報分量進行實際預報時，不需要調整預報建模的各項參數(shù)，各預報分量預報模型具有很好的穩(wěn)定性和普適性，這為開展云圖預測提供了客觀實用的新方法。

另外，還進一步做了基于相同的云圖預報因子，針對同樣的建模和預報樣本采用多元線性回歸方案進行和新方案一致的云圖預測。對比結果表明，新方案預測的云圖與實況云圖相關高的樣本更多，重構云圖的基本輪廓、紋理特征分布、清晰度以及云系強弱方面都優(yōu)于多元線性回歸方案。

值得注意的是，文中的云圖預報時效達到了72 h，具有業(yè)務實用預報意義。