王維一，裴韜，秦承志

（1.中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國家重點實驗室，北京 100101；2.中國科學(xué)院大學(xué)，北京 100049）

0 引言

隨著空間信息技術(shù)和模糊理論的發(fā)展，模糊C均值聚類（Fuzzy C－Means，F(xiàn)CM）作為地學(xué)應(yīng)用中的一種常用算法，被廣泛應(yīng)用于多圖層?xùn)鸥竦乩頂?shù)據(jù)分析中［1－4］。然而，F(xiàn)CM算法屬于計算密集型算法，計算耗時長，現(xiàn)今隨著地學(xué)中高分辨率數(shù)據(jù)的應(yīng)用、輸入圖層數(shù)的增加以及研究尺度的擴展，地學(xué)計算的數(shù)據(jù)量急劇增加，面對如此海量數(shù)據(jù)，傳統(tǒng)的串行算法在計算效率和性能方面越來越不能滿足需要，因此急需將串行FCM并行化。

自從Kwok等提出并實現(xiàn)了基于MPI的并行FCM算法用于商業(yè)數(shù)據(jù)以來［5］，地學(xué)專家、學(xué)者對FCM的并行化算法進(jìn)行了不少研究。Gong等基于MPI的FCM并行算法使人對遙感影像進(jìn)行模糊聚類分析獲得了令人滿意的線性加速比［6］；Petcu等針對多光譜遙感影像，提出了一種基于MPI的FCM并行算法，也得到較好的線性加速比［7］。但是這些方法只適合于規(guī)則的柵格數(shù)據(jù)。Liu等提出了一種基于圖形處理器（Graphics Processing Unit，GPU）的FCM 并行算法［8］。

由于地學(xué)應(yīng)用研究區(qū)域的不規(guī)則性，并行算法采用按區(qū)域大小均勻劃分（如按行、列劃分和棋盤式劃分［9］）方法導(dǎo)致各節(jié)點負(fù)載不均衡，并行效率低。針對這個問題，王少文等提出了按計算強度（即單位像元的計算量）劃分的方法，并將該方法應(yīng)用于反距離加權(quán)插值和空間熱點探測［10］。然而該方法實現(xiàn)案例較少，且主要針對單圖層。

本文針對多圖層?xùn)鸥竦乩頂?shù)據(jù)的FCM算法進(jìn)行并行化，采用按計算強度均勻劃分?jǐn)?shù)據(jù)的方式，并通過實驗與采用按區(qū)域大小均勻劃分?jǐn)?shù)據(jù)方式實現(xiàn)的FCM并行算法進(jìn)行了性能比較。

1 模糊C均值聚類算法

FCM是用隸屬度確定每個數(shù)據(jù)點屬于某個聚類程度的一種聚類算法。Bezdek于1973年提出了該算法，是早期硬C均值聚類（HCM）方法的一種改進(jìn)［11］。FCM與HCM的主要區(qū)別在于FCM用模糊劃分，使得每個給定數(shù)據(jù)點用值在［0，1］間的隸屬度確定其屬于各個組的程度，而HCM用的是｛0，1｝硬劃分。

FCM算法可表示成數(shù)學(xué)規(guī)劃問題：

其中，U為模糊隸屬矩陣，P為聚類中心，uij為模糊隸屬矩陣U中表示第j個元素屬于第i個中心的隸屬度值，m為加權(quán)指數(shù)，dij表示第j個元素與第i個中心的距離，n表示數(shù)據(jù)集個數(shù)，c為聚類中心數(shù)。

求解上述數(shù)學(xué)規(guī)劃的過程如下：

步驟一：初始化，設(shè)定聚類中心數(shù)c，2≤c≤n，n是數(shù)據(jù)集中元素的個數(shù)，最大迭代次數(shù)Lmax，設(shè)定迭代停止閾值ε＞0，初始化聚類中心P（0），設(shè)定迭代計算器k＝0。

步驟二：按式（3）更新模糊隸屬度矩陣U（k）：

步驟三：按式（4）更新聚類中心P（k＋1）：

步驟四：用一個矩陣范數(shù) ‖·‖ 比較 P（k）與P（k＋1），若‖P（k）－P（k＋1）‖≤ε，則停止迭代，并輸出模糊隸屬度矩陣U和聚類中心P，否則令k＝k＋1，轉(zhuǎn)到步驟二進(jìn)行下一次迭代。

兩首譯詩只在幾個關(guān)鍵詞的選擇和使用上不同，就構(gòu)建出截然不同的兩種語境，兩種情景模型，這是由于譯者對原詩的理解不同。從心理學(xué)的角度看，由于第一首譯詩的譯者與詩人對“遼西”的認(rèn)知環(huán)境不同，直接導(dǎo)致兩人對同一首詩的理解迥然不同。在詩人看來，“遼西”既是地名又會激活擴散到“戰(zhàn)爭”的意義，而第一首譯詩的譯者沒有這種背景知識，無法進(jìn)行激活擴散，所以無法翻譯出原詩的情景模型。

2 并行FCM算法設(shè)計

FCM算法具有計算密集型特點。在地學(xué)應(yīng)用中多用于多圖層?xùn)鸥駭?shù)據(jù)，采用數(shù)據(jù)并行模式［9］，將數(shù)據(jù)根據(jù)地理范圍劃分為數(shù)據(jù)塊，分配給不同的進(jìn)程，各個進(jìn)程之間協(xié)同并行完成聚類任務(wù)，達(dá)到高效求解的目的。

并行FCM聚類算法的具體步驟如下：

（1）數(shù)據(jù)并行讀?。簩?shù)據(jù)均衡劃分分配成numProcss（進(jìn)／線程數(shù)）份，各進(jìn)程并行讀取數(shù)據(jù)。

（2）主進(jìn)程初始化聚類中心，并將聚類中心發(fā)送到各子進(jìn)程中。

（3）各進(jìn)程按式（3）計算本地模糊隸屬度矩陣U，uij為模糊隸屬矩陣U中表示第j個元素屬于第i個中心的隸屬度值。

（4）各進(jìn)程根據(jù)式（5）、式（6），分別求取計算聚類中心（式（4））中分子Num與分母Den：

（5）利用全局規(guī)約函數(shù)AllReduce（），對各進(jìn)程中分子與分母歸并求和，求和結(jié)果分別記為Num－Sum和DenSum，然后根據(jù)式（7）求解新的聚類中心P（k＋1）：

（6）迭代步驟3－5，直到目標(biāo)函數(shù)收斂到設(shè)定的閾值。

3 顧及數(shù)據(jù)劃分負(fù)載均衡的FCM并行算法

3.1 數(shù)據(jù)劃分中的負(fù)載均衡問題

柵格地理數(shù)據(jù)的并行模式多采用數(shù)據(jù)并行模式，該模式最基礎(chǔ)的是數(shù)據(jù)劃分，傳統(tǒng)柵格數(shù)據(jù)劃分主要是按區(qū)域大小劃分（如按行、列和棋盤劃分［9］），該方法較為簡單直觀，然而由于地理學(xué)的研究區(qū)不規(guī)則，該劃分方法將使得分配給各個進(jìn)程進(jìn)行計算的數(shù)據(jù)量不均勻，從而直接導(dǎo)致負(fù)載不均衡。圖1所示的柵格數(shù)據(jù)按區(qū)域大小均勻劃分成3份，由于研究區(qū)外的柵格不參與計算，使得各進(jìn)程的計算量明顯不均勻，進(jìn)程0和進(jìn)程2的計算量明顯小于進(jìn)程1的計算量（3個進(jìn)程中計算柵格數(shù)百分比分別為20.83%、49.53%、29.64%）。

圖1 按區(qū)域大小劃分方式示例Fig.1 Data partition based on size of the area

3.2 按計算強度的數(shù)據(jù)劃分方法

針對上述問題，王少文等提出了按計算強度劃分的方法，并將該方法應(yīng)用于反距離加權(quán)插值和空間熱點探測［10］。按計算強度劃分主要思路是通過構(gòu)建計算強度估計函數(shù)，預(yù)測各計算單元的計算量，從而指導(dǎo)數(shù)據(jù)劃分，實現(xiàn)任務(wù)負(fù)載均衡［10］。該方法首先按式（8）統(tǒng)計柵格數(shù)據(jù)各行的計算量，然后進(jìn)行數(shù)據(jù)劃分，使得各進(jìn)程計算量均勻。

式中：Fy表示柵格數(shù)據(jù)第y行的計算量，m、n分別表示柵格數(shù)據(jù)的行數(shù)和列數(shù)，g（x，y）表示像元（x，y）的計算強度。

記F為計算量Fy組成的數(shù)組（數(shù)組中元素個數(shù)為m），這樣按計算強度的數(shù)據(jù)劃分就抽象成如下問題：對數(shù)組F進(jìn)行連續(xù)劃分，使得各塊元素之和中最大值最小。對于該問題，通常采用模擬退火［12］、遺傳算法［13］和蟻群算法［14］的方法尋找全局最優(yōu)解，然而對于現(xiàn)在CPU的強大計算能力，最優(yōu)解相對于次優(yōu)解對并行效率的提高并不明顯，且模擬退火、遺傳算法和蟻群算法相對復(fù)雜，故本文以下介紹一種相對簡單的求解次優(yōu)解的方法對數(shù)據(jù)進(jìn)行劃分。

3.3 按計算強度進(jìn)行數(shù)據(jù)劃分的并行FCM算法設(shè)計

本文采用按計算強度劃分?jǐn)?shù)據(jù)的方式來解決FCM并行算法中的任務(wù)負(fù)載不均衡問題。由于FCM算法中每一個研究區(qū)柵格像元的計算量一樣，研究區(qū)外柵格像元不參與計算，因此，本文給出較為簡單的計算強度估計函數(shù)：

在此基礎(chǔ)上，選擇一種相對簡單的求次優(yōu)解的方法劃分柵格數(shù)據(jù)。假設(shè)將數(shù)據(jù)劃分成K分，則只需要進(jìn)行K－1次劃分，每一次劃分都盡量保證是對所要劃分的數(shù)據(jù)進(jìn)行均勻劃分，即找出每一次劃分的位置tu（u＝1，2，3，…，K－1），使得式（10）的值最小。數(shù)組劃分示意圖如圖2。

式中：u表示第u次劃分（u＝1，2，3，…，K－1），tu表示第u次劃分的位置，當(dāng)u＝1時，記tu－1為0。

圖2 數(shù)組劃分Fig.2 Array divided

以此方法將圖1中柵格數(shù)據(jù)按計算強度均勻劃分成3份（圖3），3個數(shù)據(jù)塊中計算柵格數(shù)百分比分別為33.33%、33.33%、33.33%，可見按計算強度均勻劃分能實現(xiàn)各進(jìn)程中計算量的負(fù)載均衡。

圖3 按計算強度劃分?jǐn)?shù)據(jù)方式示例Fig.3 Data partition based on computational intensity

綜上，按計算強度進(jìn)行數(shù)據(jù)劃分的FCM并行算法流程如圖4所示。

4 基于MPI的FCM并行算法實現(xiàn)

目前，并行編程模型主要有基于消息傳遞（如Message Passing Interface，MPI）和基于共享存儲（如OpenMP）兩種。當(dāng)計算節(jié)點較多時，基于共享存儲編程模型的并行性能遠(yuǎn)不如消息傳遞編程模型［15］，因此本文選擇基于MPI實現(xiàn)FCM并行算法。上述的FCM并行算法也可基于OpenMP實現(xiàn)。

圖4 按計算強度劃分?jǐn)?shù)據(jù)的并行FCM算法流程Fig.4 Flowchart of parallel FCM algorithm

5 應(yīng)用評價

5.1 研究區(qū)數(shù)據(jù)

本文研究區(qū)域位于黑龍江省黑河市嫩江縣鶴山農(nóng)場，面積約60.2km2。數(shù)據(jù)是由該研究區(qū)0.5m分辨率的DEM數(shù)據(jù)生成的坡度、沿等高線曲率、沿坡面曲率以及地形濕度指數(shù)4個環(huán)境因子圖層數(shù)據(jù)，每個圖層由22 260×18 640個柵格組成，數(shù)據(jù)格式為.tif，總數(shù)據(jù)大小為6.4GB。圖5中黑色虛線內(nèi)是研究區(qū)。

圖5 研究區(qū)（虛線為研究區(qū)邊界）Fig.5 Study area

5.2 測試環(huán)境及評價指標(biāo)

本文測試是高性能集群環(huán)境，包括4個計算節(jié)點和1個存儲節(jié)點，每個計算節(jié)點有兩顆Intel（R）Quad Core E5645Xeon（R）CPU，共12核，內(nèi)存大小為32GB，集群操作系統(tǒng)為Linux CentOS 64位，節(jié)點間文件系統(tǒng)為Network File System（NFS）。

本文選取運行時間、加速比、并行效率3種常見的測度指標(biāo)來表征FCM并行的效果。其中，運行時間僅含計算時間，不包括I／O時間；加速比S即串行算法運行時間T串除以并行算法運行時間T并，如式（11）；并行效率E為加速比S比進(jìn)程數(shù)N，如式（12）。

FCM算法每次都隨機選擇初始聚類中心，為了具有可靠的對比性，測試時指定算法從相同的初始聚類中心開始迭代，且所有測試迭代次數(shù)均相同（本實驗為40次）。

5.3 評價結(jié)果

FCM的串行算法運行時間為1 130min（約19 h）。如圖6所示，兩種并行算法均大大縮短了計算時間，按計算強度均勻劃分的并行算法比按區(qū)域大小均勻劃分的并行算法的效率高。隨著進(jìn)程數(shù)增加，無論是按區(qū)域大小均勻劃分還是按計算強度均勻劃分，并行算法的運行時間開始都明顯減少，隨后趨于穩(wěn)定。當(dāng)進(jìn)程數(shù)為48h，按區(qū)域大小均勻劃分和按計算強度均勻劃分的并行算法運行時間分別減少到41min（約為串行算法的1／28）和28min（約為串行算法的1／40），表明并行算法實現(xiàn)了該問題的快速、高效求解。

按區(qū)域大小均勻劃分的并行算法性能明顯不如按計算強度均勻劃分的并行算法（圖7）。從加速比（圖7a）看，按區(qū)域大小均勻劃分和按計算強度均勻劃分的并行算法加速比都呈線性增長，然而前者加速比明顯低于后者，當(dāng)進(jìn)程數(shù)為48h，前者加速比只有28左右，后者達(dá)到40左右。另外，從并行效率（圖7b）看，按區(qū)域大小均勻劃分和按計算強度均勻劃分的并行算法加速效率都呈遞減趨勢，而后者的并行效率穩(wěn)定地高于前者（20%以上）。

圖6 算法運行時間對比Fig.6 The runtime of algorithm

圖7 加速比與并行效率Fig.7 The speedup ratio and parallel efficiency

6 結(jié)論

本文研究了地理柵格數(shù)據(jù)FCM算法的并行化，針對并行化時按區(qū)域大小均勻劃分?jǐn)?shù)據(jù)的方式導(dǎo)致的任務(wù)（計算量）負(fù)載不均衡問題，引入了按計算強度均勻劃分的方法，針對全局最優(yōu)解算法復(fù)雜且相對次優(yōu)解加速效率提高不明顯的特點，設(shè)計了一種簡單求次優(yōu)解的方法來實現(xiàn)數(shù)據(jù)劃分，并利用MPI實現(xiàn)了并行算法。實驗結(jié)果表明，所設(shè)計的FCM并行方法大大縮短了計算時間，加速性能明顯，按計算強度劃分的并行算法相較于按區(qū)域大小劃分的并行算法，加速效率更加明顯。本文實現(xiàn)的按計算強度均勻劃分的方法適用于研究區(qū)域不規(guī)則以及數(shù)據(jù)中存在大量無值區(qū)的情況。

［1］ 胡姝婧，胡德勇，趙文吉.基于LSMM和改進(jìn)的FCM提取城市植被覆蓋度——以北京市海淀區(qū)為例［J］.生態(tài)學(xué)報，2010，30（4）：1018－1024.

［2］ 楊琳，朱阿興，李寶林，等.應(yīng)用模糊c均值聚類獲取土壤制圖所需土壤－環(huán)境關(guān)系知識的方法研究［J］.土壤學(xué)報，2007（5）：784－791.

［3］ 邱超.模糊聚類分析在水文預(yù)報中的研究及應(yīng)用［D］.浙江大學(xué)，2007.

［4］ 蔣衛(wèi)國，陳強，郭驥，等.基于HPSO和FCM的多光譜遙感圖像濕地分類［J］.光譜學(xué)與光譜分析，2010，30（12）：3329－3333.

［5］ KWOK T，SMITH K，LOZAN S，et al.Parallel fuzzy c－means clustering for large data sets［A］.Euro－Par 2002Parallel Processing Proceedings［C］.2002，2400：365－374.

［6］ GONG X J，CI L L，YAO K Z.A FCM algorithm for remotesensing image classification considering spatial relationship and its parallel implementation［C］.2007International Conference on Wavelet Analysis and Pattern Recognition，2007.994－998.

［7］ PETCU D，ZAHARIE D，PANICA S，et al.Fuzzy clustering of large satellite images using high performance computing［A］.High－Performance Computing in Remote Sensing［C］.2011.

［8］ LIU G，LIANG X，HE X.Graphics processing unit based Fuzzy C－Means clustering segmentation［J］.Computer Science，2012，39（1）：285.

［9］ 王結(jié)臣，王豹，胡瑋，等.并行空間分析算法研究進(jìn)展及評述［J］.地理與地理信息科學(xué)，2011，27（6）：1－5.

［10］ WANG S W，ARMSTRONG M.A theoretical approach to the use of cyberinfrastructure in geographical analysis［J］.International Journal of Geographical Information Science，2009，23（2）：169－193.

［11］ BEZDEK J C.Pattern Recognition with Fuzzy Objective Function Algorithms［M］.Plenum Press，1981.7－10.

［12］ KIRKPATRICK S，GELATT C D，VECCHI M P.Optimization by simulated annealing［J］.Science，1983，220（4598）：671－680.

［13］ LORIES G.Toward a practice of autonomous systems：Proceedings of the first European conference on Artificial life［J］.Behavioural Processes，1996，37（2－3）：257－258.

［14］ DORIGO M，MANIEZZO V，COLORNI A.Ant system：Optimization by a colony of cooperating agents［J］.Ieee Transactions on Systems Man and Cybernetics Part B－Cybernetics，1996，26（1）：29－41.

［15］ 張林波，遲學(xué)斌.并行計算導(dǎo)論［M］.北京：清華大學(xué)出版社，2006.3－6.