柵格地理數(shù)據(jù)模糊C均值聚類算法的并行化研究

王維一，裴韜，秦承志

（1.中國(guó)科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100101；2.中國(guó)科學(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算法屬于計(jì)算密集型算法，計(jì)算耗時(shí)長(zhǎng)，現(xiàn)今隨著地學(xué)中高分辨率數(shù)據(jù)的應(yīng)用、輸入圖層數(shù)的增加以及研究尺度的擴(kuò)展，地學(xué)計(jì)算的數(shù)據(jù)量急劇增加，面對(duì)如此海量數(shù)據(jù)，傳統(tǒng)的串行算法在計(jì)算效率和性能方面越來(lái)越不能滿足需要，因此急需將串行FCM并行化。

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

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

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

1 模糊C均值聚類算法

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

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

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

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

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

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

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

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

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

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

FCM算法具有計(jì)算密集型特點(diǎn)。在地學(xué)應(yīng)用中多用于多圖層?xùn)鸥駭?shù)據(jù)，采用數(shù)據(jù)并行模式［9］，將數(shù)據(jù)根據(jù)地理范圍劃分為數(shù)據(jù)塊，分配給不同的進(jìn)程，各個(gè)進(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）計(jì)算本地模糊隸屬度矩陣U，uij為模糊隸屬矩陣U中表示第j個(gè)元素屬于第i個(gè)中心的隸屬度值。

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

（5）利用全局規(guī)約函數(shù)AllReduce（），對(duì)各進(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ù)載均衡問(wèn)題

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

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

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

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

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

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

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

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

在此基礎(chǔ)上，選擇一種相對(duì)簡(jiǎn)單的求次優(yōu)解的方法劃分柵格數(shù)據(jù)。假設(shè)將數(shù)據(jù)劃分成K分，則只需要進(jìn)行K－1次劃分，每一次劃分都盡量保證是對(duì)所要?jiǎng)澐值臄?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時(shí)，記tu－1為0。

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

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

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

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

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

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

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

5 應(yīng)用評(píng)價(jià)

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

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

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

5.2 測(cè)試環(huán)境及評(píng)價(jià)指標(biāo)

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

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

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

5.3 評(píng)價(jià)結(jié)果

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

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

圖6 算法運(yùn)行時(shí)間對(duì)比Fig.6 The runtime of algorithm

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

6 結(jié)論

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

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

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

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

［4］ 蔣衛(wèi)國(guó)，陳強(qiáng)，郭驥，等.基于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)展及評(píng)述［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é)斌.并行計(jì)算導(dǎo)論［M］.北京：清華大學(xué)出版社，2006.3－6.