董怡琦，周明睿，劉力源，余一冬，陳 科，童亞拉，2

(1.湖北工業(yè)大學理學院，武漢 430068；2.湖北省能源光電器件與系統(tǒng)工程技術研究中心，武漢 430068)

自20世紀初以來，數值天氣預測作為一種預測方法，受到了廣泛的關注.隨著高性能計算機的誕生，高分辨率數值模式MM5的使用也成為一種熱門手段.Bjerknes[1]提出從數值預測整體來看，其中重要一個環(huán)節(jié)在于初始場，其數據是否精準對最終輸出影響極大.MM5中的資料同化是一種變分資料同化，屬于非線性優(yōu)化問題.近年來，為了改進變分同化準確性和收斂性效率，粒子群算法(PSO)、遺傳算法(GA)等經典群智能算法被引到資料同化中.王順風、白晨等[2-3]分別提出了兩種遺傳算法(GA)同化方案，鄭琴[4]從PSO著手，改變慣性權重的變換策略，選擇合適的學習因子，設計了動態(tài)慣性權重粒子群算法(PSODIWAF)，該算法能有效解決含有不連續(xù)開關過程的變分資料同化問題.為滿足對精度和同化時間的要求，各種改進的粒子群算法被利用到資料同化中，李成、張頂學[5-6]等對陷入局部最優(yōu)的問題進行了改進；劉明、袁小平[7-8]等采用不同的學習機制，通過效率更高的迭代公式，在收斂速度上作出了改進，減少了收斂時間；何莉等[9]重點從運行時間進行突破，采用多核處理能力提高了算法處理速度，其采用的是粒子隨機替換策略，為了保證整個群體能夠進行信息交流，同時也設置了共享區(qū)域；鄭琴[10]和張成興[11]分別設計了動態(tài)慣性權重粒子群算法(PSODIWAF)和時變雙重壓縮因子粒子群算法(PSOTVCF)將其用于資料同化.盡管在收斂精度上兩種算法都有了較大的改進，但同化時間仍然存在缺陷.因此，受文獻[12]的啟發(fā)，本文設計了一種改進的并行算法.該算法通過并行計算對粒子群完成分組，采用PSOTVCF進行深度搜索，同時將具有不連續(xù)“開關”過程的變分資料同化作為應用場景，建立了一種新的同化模型.在幾乎沒有降低同化精度的前提中，極大地提高了同化速度，減少了同化時間.

1 資料同化和粒子群算法簡介

1.1 資料同化的模式方程

文獻[4]中的偏微分方程被用作資料同化處理中的控制方程：

(1)

(2)

模擬參數化過程中的“on-off”開關.可將(1)式轉化為其數值模式：

(3)

式中，空間格點用i表示，空間步長用Δl表示，且有l(wèi)i=iΔl.Δt為時間步長，tk=kΔt，k為空間層.N=T/Δt為積分過程中總的時間層；M+1=(L/Δl)+1為空間離散點總數.

1.2 粒子群算法

在PSO中，D維空間的搜索個體用粒子表示，粒子擁有位置和速度兩個屬性.其每個時刻所在的位置被看作是待求解問題的一個候選解，粒子位置的變換認為是最優(yōu)解的尋找過程.以及其移動速度是受截止到上一時刻自身所經歷過的最優(yōu)位置和整個群體到達最優(yōu)位置的影響，每次迭代進行相應調整.

用下列公式對粒子進行操作：

(4)

(5)

迭代終止條件沒有固定的要求，通常規(guī)定為當前定位到的最優(yōu)位置已經達到了設定的最小適應閾值或者達到了設定的最大迭代次數.

2 改進的并行粒子群算法

2.1 設計思路

文獻[11]中設計了一種時變雙重壓縮因子粒子群算法(PSOTVCF)，在慣性權重的基礎上引入雙重壓縮因子進行迭代.雖然提高了同化精度，但也因此使收斂速度變慢，同化時間增加.為解決同化時間過長的問題，本文受文獻[12]的啟發(fā)，將粒子群分成四個子集，分配多核CPU中一個核給每個子集進行計算，迭代的同時粒子與其他子集進行信息交流，從而每個粒子都能較好地掌握整個種群的迭代狀況，并且利用PSOTVCF進行粒子的更新.

2.2 時變雙重壓縮因子粒子群算法的并行處理

一般地，并行問題的處理有兩種方法：(1)功能分解形式，對問題自身進行拆分，將一個大問題拆分成幾個子問題，然后再加以解決；(2)域分解形式，對問題區(qū)域進行分解，對分解后的區(qū)域并行處理.P2PSO中采用區(qū)域分解形式進行并行，同時引入協(xié)同進化的策略.首先根據CPU核數確定分為4組，然后在主線程中依據分組數設定4個線程，再對每個線程進行任務分配.并行中信息交流也有兩種方法：同步和異步.同步形式是4個線程內粒子都完成了位置速度更新后，再評估和更新全局最優(yōu)個體的位置和速度.異步形式[12]是只要有其中一個線程內完成了個體最優(yōu)粒子位置速度的更新，就馬上評估和更新全局最優(yōu)粒子信息.P2PSO算法屬于并行算法，本質上是模擬鳥群捕食行為.因鳥群捕食過程中，只要有一只鳥發(fā)現了食物，便會通過一些形式去告知其他的鳥，而不是等所有的鳥都找到食物再去通知，這就是“即時種群通信行為”.因此該算法在信息交流的并行化中采用異步形式，每個線程完成了組內粒子的更新后，申請獲取讀寫同步鎖來更新整個種群內全局最優(yōu)粒子信息.

2.3 并行時變雙重壓縮因子粒子群算法

算法引入PSOTVCF進行深度搜索，慣性權重ω用壓縮因子取代，緩解全局與局部探尋間的矛盾關系.其中通過加速因子計算得到壓縮因子，加速因子的計算公式為

(6)

其中，C1N是第一個加速因子的最大值，C1M是第一個加速因子的最小值，C2N是第二個加速因子的最大值，C2M是第二個加速因子的最小值.MAXITER是設定的最大迭代數，ITER是當前迭代次數.

PSOTVCF中粒子速度更新公式表示：

v(t+1)=χ(v(t)+φy(t))，

(7)

y(t+1)=χv(t)+(1-χφ)y(t)，

(8)

由此得到系統(tǒng)矩陣：

(9)

φ=C1+C2，

(10)

滿足條件φ>4，由收斂準則得到壓縮因子的定義：

(11)

其中，壓縮因子χ為正實數，且χ∈(0，1).

單一壓縮因子[13]在高維多峰函數應用中精度較低且平衡全局和局部搜索中能力極為限制，所以采用雙重壓縮因子.

速度計算公式轉換為：

v(t+1)=χ1(v(t)+C1Rand

(Pi-x(t)+C2Rand(Pg+x(t))))，

(12-1)

v(t+2)=χ2?(v(t+1)).

(12-2)

χ1和χ2均由公式(10)和(11)計算，其中χ1中加速因子取得是設定的初始值，而χ2中采用的是時變的加速因子，先由公式(6)計算加速因子再代入公式(10)(11)中.

設循環(huán)次數為N，得到最終簡化的粒子速度：

V(N)=(χ2?χ1)N-1(V(0)+φY).

(13)

時變雙重壓縮因子粒子群算法有效地改進了粒子群算法容易陷入局部最優(yōu)的問題.

在變分資料同化應用過程中，其中代價函數通過該公式進行計算：

(14)

其中

(15)

表示q0滿足物理約束和相容性條件，這里q為q0代入模式方程(1)的解.

式(14)離散化后得到離散化函數：

(16)

PSOTVCF算法中采用的適應度函數是：

(17)

其中

(18)

P2PSO的算法流程如下：

1)對算法參數進行初始化，設定整個群體粒子數P，慣性權重ω，加速因子，分組數C以及在各個分組中的最大迭代次數I.

2)讀寫同步鎖設定，按照設定的分組數依次建立線程，在各個線程中確定組內粒子數Pi和終止條件最大迭代次數I.

5)分組線程申請獲取讀寫同步鎖.

6)將更新后的組內最優(yōu)解pgi與此時整個種群的最優(yōu)解pg進行比較，若pgi優(yōu)于pg，則令pg=pgi.

7)分組線程釋放讀寫同步鎖.

9) 判斷是否滿足迭代終止條件，如果組內迭代數達到預先設置的最大迭代數I，則停止該線程的計算，若不滿足則轉為步驟(4).

10)直至所有線程迭代結束，產生結果.

3 數值實驗結果分析

3.1 仿真環(huán)境

采取文獻[11]的實驗數據和實驗分析方法，將三種算法——PSOTVCF、PSODIWAF和P2PSO分別應用于變分資料同化，在精度和時間上進行比較.其中，在P2PSO中，慣性權重ω從0.7線性遞減到0.1，C1N=2.98，C1M=2.78，C2N=1.55，C2M=1.35，種群中含有200個粒子，設定迭代最大次數為200，測試環(huán)境為硬件Intel Core i5，軟件MATLAB R2017a.

3.2 收斂精度

圖1表示的是在迭代200次時，三種算法在200次同化實驗后的平均收斂精度圖.圖中，自上而下的橫線分別表示PSODIWAF，P2PSO，PSOTVCF，縱坐標是收斂精度的對數log10J，其值越小表示同化質量越好，同化結果越準確.從實驗結果可直觀看出與PSODIWAF相比，P2PSO和PSOTVCF在同化實驗結果中更為準確.兩者平均收斂精度相近，但P2PSO略差與PSOTVCF.為了明確三種算法在運行過程中的差距緣由，需要繼續(xù)對收斂趨勢進行分析.

圖1 三種算法在200次同化實驗中平均收斂精度圖Fig.1 Average convergence accuracy of three algorithms in 200 assimilation experiments

圖2表示的是PSODIWAF、PSOTVCF、P2PSO在迭代200次時收斂精度的變化趨勢.圖中橫坐標代表迭代次數，縱坐標同圖1含義相同，即為收斂精度的對數log10J.星號標出的點分別為迭代次數50，100，150，200.折線從終點來看自上而下分別表示PSODIWAF，P2PSO，PSOTVCF.分段分析，在同化初期(迭代次數為50)，PSOTVCF與P2PSO結果大致相同，收斂精度均略低于PSODIWAF；同化中期(迭代次數為100)，三者同化質量相似；同化后期(迭代次數為150)，PSOTVCF和P2PSO后來居上，同化質量最好，PSODIWAF次之，并且粒子基本已經收斂；實驗終點(迭代次數為200)，PSOTVCF同化質量最好，其次是P2PSO，最后是PSODIWAF，其收斂數值依次為-13、-12.9、-11.2.從整體來看，P2PSO與原算法PSOTVCF在同化初期基本保持一致，自同化中期開始，兩者產生細微差距.因此P2PSO同化結果的質量高于動態(tài)權重粒子群算法，與PSOTVCF相似.并行算法的精度略低于原算法的原因，應該是由于在并行過程中信息交流的程度沒有單線程的高，從而導致收斂精度降低.

圖2 三種算法收斂精度變化趨勢對比圖Fig.2 Comparison of convergence accuracy trends of three algorithms

3.3 同化時間

將PSOTVCF、PSODIWAF和P2PSO應用于資料同化問題中，同時記錄三種算法在迭代次數為50，100，150，200時所花費的時間.因為粒子群算法具有一定的隨機性，為了使數據更加有說服力，每一個同化窗口進行6組實驗.將6組試驗數據平均值記錄于表1中.

表1 三種算法在迭代次數為50，100，150，200時平均同化時間(s)對比表Tab.1 Comparison of average assimilation time(seconds) of three algorithms with iteration times of 50,100,150,200 s

通過表中數據，直觀上可得到無論迭代次數是在50，100，150還是200，P2PSO的同化速度都要比PSOTVCF和PSODIWAF快，同化時間縮短一半左右.P2PSO迭代200次的時間比PSOTVCF和PSODIWAF迭代100次的時間還要短.數據上計算得出在200代迭代之后P2PSO的收斂速度和PSODIWAF相比提高了51.5%左右，和PSOTVCF相比提高了55.1%左右.因此，本文所設計的并行粒子群算法通過多個CPU的聯合計算對于提高同化速度，縮短同化時間具有顯著的效果.

4 結論與展望

本文著重在時間上改進了粒子群算法，結合計算機CPU的多核處理硬件知識，在不過度影響精度的情況下較大程度上加快了速度.在變分資料同化中采用所設計的算法，減少了同化時間，與時變雙重壓縮因子PSO和動態(tài)權重PSO算法相比在收斂速度上有了較大的提升，但在收斂精度上仍需改進.因此后續(xù)可針對因多線程導致粒子間信息交流程度略有降低的現象，對算法進一步改進，實現較大程度加快同化速度的同時也能夠進一步提高收斂精度.本文將含“開關”過程的變分資料同化作為該算法的應用場景，后續(xù)可嘗試著將P2PSO應用于其他場景下的資料同化，擴大所設計算法的適用空間.