鮑軍鵬， 楊科， 周靜

(1. 西安交通大學(xué)電子與信息工程學(xué)院, 西安 710049； 2. 寧夏軍區(qū), 銀川 750021)

時序數(shù)據(jù)是一種帶有時間標記的常見大數(shù)據(jù)類型。航天、氣象、交通、電力、工業(yè)、金融、科研[1]等眾多領(lǐng)域日復(fù)一日地不斷產(chǎn)生此類型數(shù)據(jù)。在大數(shù)據(jù)和云計算時代，人們需要運用數(shù)據(jù)挖掘方法分析研究這些時序大數(shù)據(jù)，以獲得蘊藏在數(shù)據(jù)背后的異常、關(guān)聯(lián)、模式、趨勢等知識，并利用所得知識進行異常檢測、關(guān)聯(lián)分析、故障診斷、監(jiān)測預(yù)警等應(yīng)用，從而體現(xiàn)出數(shù)據(jù)的價值。然而，時序大數(shù)據(jù)應(yīng)用卻面臨著大數(shù)據(jù)量處理和實時性要求的挑戰(zhàn)。如何利用各種優(yōu)化手段有效提升海量數(shù)據(jù)挖掘效率，縮短運行時間，是當前時序數(shù)據(jù)應(yīng)用領(lǐng)域重點研究的一個問題。

針對衛(wèi)星時序數(shù)據(jù)的挖掘方法有很多，包括異變過程多類型特征分析、異變過程及特征變化規(guī)律獲取[2-3]、多參數(shù)關(guān)聯(lián)知識挖掘、狀態(tài)關(guān)聯(lián)知識挖掘[4-5]、狀態(tài)異常形態(tài)挖掘[6-7]、異常和故障模式挖掘[8-9]、多誘因復(fù)合狀態(tài)異變檢測[10-11]、故障診斷等。這些方法在執(zhí)行過程中涵蓋了數(shù)據(jù)預(yù)處理、表示、分割、相似性度量、分類、聚類等多種算法。采取有效措施優(yōu)化這些算法，縮短運行時間對于提高系統(tǒng)性能至關(guān)重要。

在更短時間內(nèi)處理更多數(shù)據(jù)，一直是計算機科學(xué)追求的目標之一。該目標有多種實現(xiàn)策略，包括粗粒度的分布式計算及節(jié)點級較細粒度的并行計算。本文僅討論用Python語言處理航天大數(shù)據(jù)的節(jié)點級并行與優(yōu)化策略，并不涉及分布式計算。本文工作能夠提高航天大數(shù)據(jù)挖掘效率，有利于提高衛(wèi)星運行狀態(tài)實時監(jiān)控效能，對于保障衛(wèi)星長期穩(wěn)定在軌運行具有重要意義。

1 程序并行與優(yōu)化方法

單個計算節(jié)點上的并行計算包括CPU+GPU等4種典型架構(gòu)[12]。本文用到的主要方法有基于多核CPU的并行計算、基于CPU+GPU的異構(gòu)計算、串行算法優(yōu)化與向量化等。

Python語言存在全局解釋鎖(Global Interpreter Lock，GIL)問題，故在多核CPU上采用多進程而不是多線程進行并行計算。使用Python的multiprocessing并行庫，可以并行開啟多個進程，沒有數(shù)量限制。每個進程都運作各自的GIL。

GPU(Graphics Processing Unit)以單指令多線程流(SIMT)方式運行，其顯著特點有[13]：更高的內(nèi)存帶寬；更多的計算單元；更具性價比優(yōu)勢。GPU在浮點運算及并行計算性能方面可優(yōu)于CPU達數(shù)十倍乃至上百倍。因此，基于GPU的并行計算已經(jīng)成為很多領(lǐng)域應(yīng)用中的主流，如天文觀測[14]、航空計算[15]、氣象預(yù)報[16]、水文模擬[17]、軍事仿真[18]、生物研究[19]等。另外還有一大部分應(yīng)用采用CPU+GPU異構(gòu)平臺。

除了硬件加速，對串行算法本身進行向量化改進也是一種廣泛采用的代碼優(yōu)化方法。在不使用GPU的情況下，這種方法也能大幅提升矩陣或向量運算速度。

衛(wèi)星數(shù)據(jù)是一種典型的時序大數(shù)據(jù)，其基本優(yōu)化思路是：分析數(shù)據(jù)挖掘過程中各種算法的性能瓶頸，找出耗時大、時間復(fù)雜度、具有邏輯或數(shù)據(jù)獨立性、可并行化的部分，在多核CPU或GPU上執(zhí)行；對于耗時小或不可并行部分則仍在CPU上串行執(zhí)行[20]，并結(jié)合向量化等方法做進一步優(yōu)化，以獲得最優(yōu)效果。

2 異變過程多類型特征分析

時序數(shù)據(jù)分析一般是在訓(xùn)練樣本上找出數(shù)據(jù)統(tǒng)計特性和發(fā)展規(guī)律性，構(gòu)建時序數(shù)據(jù)模型，然后進行樣本外預(yù)測。同理，在衛(wèi)星時序數(shù)據(jù)挖掘中，一些方法是為了挖掘頻繁出現(xiàn)的模式，期望發(fā)現(xiàn)某種規(guī)律，異常數(shù)據(jù)被作為噪聲忽略；而另外一些方法則更重視異常數(shù)據(jù)背后可能隱藏的重要信息，如挖掘異常檢測、故障預(yù)測等模式在衛(wèi)星應(yīng)用中有很高價值。

如圖1所示，異變過程多類型特征(Multiple Feature of Anomaly Process，MFAP)分析是一個典型的衛(wèi)星時序數(shù)據(jù)挖掘過程，包括了數(shù)據(jù)預(yù)處理、特征提取、機器學(xué)習算法等基本步驟。本文以該過程為代表，針對其中求周期、提特征、聚類等3個關(guān)鍵步驟，討論其并行與優(yōu)化方法。本文方法完全可以推廣到其他類似過程。

圖1 異變過程多類型特征分析流程Fig.1 MFAP analysis flowchart

3 自適應(yīng)獲取周期優(yōu)化

3.1 自適應(yīng)獲取周期算法

對于明顯周期性變化的數(shù)據(jù)，其最小正周期是一個非常重要的信息。在數(shù)據(jù)挖掘過程中，數(shù)據(jù)周期可以人為設(shè)定，但是人為指定法缺乏適應(yīng)能力。一旦衛(wèi)星參數(shù)發(fā)生變化，難以及時響應(yīng)到正確周期，這不但增加了系統(tǒng)維護負擔，而且增加了設(shè)定錯誤的風險。而自適應(yīng)獲取周期算法根據(jù)數(shù)據(jù)序列自動識別出該參數(shù)的最小完整周期，不必人工逐一測算。若直接對原始數(shù)據(jù)用傅里葉變換求周期容易受到噪聲干擾，因此本文首先求取數(shù)據(jù)窗口的相關(guān)性，然后在相關(guān)性(相似度)序列的基礎(chǔ)上再求周期。算法步驟如下：

步驟1設(shè)置觀察向量長度為L，選擇從第1條數(shù)據(jù)開始的L條數(shù)據(jù)作為基準向量v0。

步驟2從起點依次往后移動kΔt，生成一系列等長度的偏移向量vk，k=1,2,…，直到向量的長度小于L。

步驟3計算基準向量v0與偏移向量vk間的相似度(內(nèi)積)，得到相似度列表。

步驟4利用傅里葉變換獲取相似度列表的周期。

步驟5根據(jù)數(shù)據(jù)量、最大能量的頻率、采樣間隔Δt及傅里葉變換公式，得到原始數(shù)據(jù)周期。

3.2 串行算法優(yōu)化與向量化

3.2.1 偏移向量的局部更新

自適應(yīng)獲取周期算法的步驟2生成了一系列滑動窗口。一般滑動窗口都有重疊，即滑動偏移量通常小于劃分窗口大小。其實，每次窗口偏移后并不需要更新窗口內(nèi)的所有元素(見圖2)，只需向前次偏移向量中添加δ個新元素，再舍棄δ個舊元素即可。其中δ為偏移量大小。這稱為局部更新。顯然元素重疊部分占比越高，計算性能提升越明顯。

圖2 偏移向量局部更新示意圖Fig.2 Partial update diagram of offset vector

3.2.2 numpy向量化

Python的numpy包提供了經(jīng)過優(yōu)化的向量運算庫。自行編寫循環(huán)逐點計算的效率要遠遠低于使用numpy包向量化運算的效率。所以在程序中應(yīng)當盡量消除循環(huán)，而使用Python的numpy相應(yīng)向量化方法，同時更新向量中的多個元素，從而提高窗口向量整體的更新速度。在自適應(yīng)獲取周期算法的步驟3中，向量間內(nèi)積也應(yīng)用向量化計算方法numpy.dot實現(xiàn)。

3.2.3 計算相似度序列的改進

根據(jù)局部更新思想，計算相似度列表的方法也可進行優(yōu)化。窗口重疊部分數(shù)據(jù)不需要重復(fù)計算，只需計算一次。也就是說，在自適應(yīng)獲取周期算法的步驟3中無需每次偏移后都在全部維度上重新計算偏移向量與基準向量的點積，而只需在前一次偏移時先記錄好點積結(jié)果，后一次偏移時加上更新部分的元素點積，同時減去舍棄部分的元素點積。這樣便省去了重復(fù)計算重疊部分數(shù)據(jù)點積的開銷。算法復(fù)雜度由此大為降低，并且窗口重疊部分占比越高性能提升越大。

3.2.4 實驗結(jié)果

本文中的所有實驗均采用以下平臺。GPU：Tesla K40c(2880 cores)；CPU：Xeon(R) CPU 2.30 GHz×20；Memory：15.6 GB；Disk：5.9 TB；OS：ubuntu 16.04.3 desktop x64；Python：v3.5；Cuda：v8.0。

表1對比了自適應(yīng)獲取周期算法的串行代碼在改進優(yōu)化前后的運行時間，其中數(shù)據(jù)大小是經(jīng)過預(yù)處理之后數(shù)據(jù)序列所含數(shù)據(jù)點的個數(shù)?？梢钥闯?，當窗口偏移量相同時，數(shù)據(jù)量越大，加速比越高。

表1 自適應(yīng)獲取周期算法串行代碼優(yōu)化前后耗時對比Table 1 Comparison of adaptive cycle achieving algorithm’s time consumption before and after serial optimization

3.3 基于多核CPU的多窗口向量相似度獲取

3.3.1 優(yōu)化方法

觀察窗口依次平移并獲取偏移向量，目的都是為了計算每個偏移向量與基準向量間的相似度，再匯總為相似度列表并據(jù)此算出周期。這在原始串行程序中體現(xiàn)為一個高次循環(huán)內(nèi)反復(fù)調(diào)用偏移向量獲取函數(shù)及相似度計算函數(shù)的過程，是影響周期獲取快慢的主要因素，應(yīng)做重點優(yōu)化。

優(yōu)化后的程序流程如圖3所示。通過使用Python的multiprocessing庫實現(xiàn)并行獲取多窗口向量相似度。進程池中的任務(wù)分批執(zhí)行，每一批可同時執(zhí)行cores個任務(wù)，從而能夠在整體上加快相似度列表的獲取速度。不過需注意，由于這種方法忽略了依次獲取窗口向量時的前后依賴關(guān)系，所以不能使用3.2.1節(jié)局部更新的優(yōu)化方法，而應(yīng)多個進程同時讀取含有重疊部分的初始向量。

圖3 基于多核CPU并行的相似度列表獲取流程Fig.3 Acquisition flowchart of similarity list based on multi-core CPU parallelization

3.3.2 實驗結(jié)果

針對同一數(shù)據(jù)用例，實驗中設(shè)置了不同偏移量δ(決定了計算量)，考察征用不同數(shù)目CPU核心完成并行計算的耗時,結(jié)果如圖4所示。

表2對比了多核CPU并行的最優(yōu)耗時與串行代碼優(yōu)化后的耗時?？梢钥吹?，并行優(yōu)化的效果并非總是優(yōu)于串行優(yōu)化，而與偏移量δ有關(guān)。

原因在于Python以多進程方式實現(xiàn)多核CPU并行計算過程，其最大限制是非內(nèi)存共享，即創(chuàng)建幾個函數(shù)進程就要復(fù)制并傳入幾份相同的參數(shù)，這會消耗大量內(nèi)存和時間。故應(yīng)根據(jù)實際計算量設(shè)置合理進程數(shù)，使得程序并行化后“節(jié)省的計算時間超過復(fù)制數(shù)據(jù)時間”，才能得到更優(yōu)的結(jié)果。經(jīng)測試，最優(yōu)條件為“活動核心數(shù)最少且每個核心的利用率接近滿載”，亦即進程內(nèi)計算量與進程間通信量的比值越高越好。因此在該實驗中看到，偏移量值越小時生成的向量越多，實際計算量越高，此時用CPU多核并行效果更好，否則有可能不如串行優(yōu)化。

圖4 多窗口向量相似度獲取耗時Fig.4 Time consumption of similarity acquisition of multiple window vectors

表2 不同方法優(yōu)化前后耗時結(jié)果對比Table 2 Comparison of time consuming results before and after different optimization methods

3.4 基于GPU的單一窗口向量相似度計算

3.4.1 優(yōu)化方法

兩向量(基準向量與偏移向量)間的相似度計算屬于數(shù)值密集型計算，特別是向量長度很大時，串行計算點積的時間就會很長。本文利用PyCuda將Cuda代碼直接嵌入Python程序中，實現(xiàn)基于GPU的多線程并行，從而大幅縮短計算耗時。設(shè)備端代碼流程如圖5所示。

求和中使用原子操作和共享內(nèi)存的組合可以避免內(nèi)核二次調(diào)用。如果硬件支持操作符⊕的原子操作，那么就可以簡單地使用單遍規(guī)約操作。如對于加法操作，只需調(diào)用atomicAdd()將塊中部分結(jié)果加到全局內(nèi)存中即可。如果硬件不支持，則采取兩遍規(guī)約，啟動兩次GPU內(nèi)核。這主要是針對Cuda線程塊無法同步問題的解決方法。實際上，其異步執(zhí)行效率并不低。

圖5 基于GPU并行的相似度計算流程Fig.5 Flowchart of similarity calculation based on GPU parallelization

3.4.2 實驗結(jié)果

實驗中選取不同長度的向量作為數(shù)據(jù)用例，并分別采用CPU(包含pure python及numpy兩種模式)和GPU計算兩向量間的相似度，耗時結(jié)果如圖6所示。

對于相似度計算來說，經(jīng)過向量化的代碼在多核CPU上的執(zhí)行時間與經(jīng)過Cuda并行的代碼在GPU上的執(zhí)行時間不相上下。不過該結(jié)果也并非一成不變，而與主機和設(shè)備端的硬件性能、數(shù)據(jù)量、計算復(fù)雜度等因素均有關(guān)，實際中應(yīng)根據(jù)具體情況選擇最優(yōu)策略。

圖6 單一窗口向量相似度計算耗時Fig.6 Time consumption of single window vector similarity calculation

4 特征向量提取與計算優(yōu)化

4.1 特征向量提取與計算方法

分析時序數(shù)據(jù)特征時，一般要把數(shù)據(jù)劃分成連續(xù)觀察窗口。對于周期數(shù)據(jù)，窗口大小一般為周期的整倍數(shù)。而對于非周期數(shù)據(jù)，則需要人工指定窗口大小，如128。如果窗口中數(shù)據(jù)不完整，有缺失，則可以采用刪除策略，直接刪除該段窗口，如數(shù)據(jù)盲區(qū)前后的數(shù)據(jù)，以及最后一個不完整的窗口等。

窗口劃分完之后，需要獲取每個觀察窗口的多種不同特征，主要包括統(tǒng)計特征、傅里葉特征、PCA特征和小波特征等。通過提取觀察窗口的多類型特征，可以從不同角度獲取數(shù)據(jù)特征信息，更有利于發(fā)現(xiàn)隱藏在數(shù)據(jù)中的規(guī)律或者模式。

4.2 基于多核CPU的多窗口向量特征提取

4.2.1 優(yōu)化方法

4.1節(jié)所述多種特征可以各自構(gòu)成單一特征向量，也可以經(jīng)過組合之后構(gòu)成合成特征向量。各種特征的計算過程相互獨立，互不依賴，完全可以并行執(zhí)行，而不必串行循環(huán)。即可以用cores個進程同時獲取多個窗口的特征，實現(xiàn)較大粒度并行，其中cores為CPU核心數(shù)。

本文利用Python的multiprocessing庫將4.1節(jié)串行程序替換為多核CPU并行程序。與圖3中計算多個向量間相似度過程不同，不用人工設(shè)定每批執(zhí)行的進程數(shù)，直接自動獲取CPU最大邏輯核心數(shù)作為最優(yōu)進程數(shù)即可。因為點積運算相對簡單，CPU會很快完成計算過程。如果進程數(shù)過多，單個進程的數(shù)據(jù)量偏少，就有可能造成在進程間切換，復(fù)制傳遞數(shù)據(jù)的消耗反而大于計算數(shù)據(jù)的消耗。但是在特征提取過程中，傅里葉變換、小波變換、PCA計算等復(fù)雜度比較高，故無需擔心進程數(shù)設(shè)置過大導(dǎo)致單個CPU核心負載過低的問題。

4.2.2 實驗結(jié)果

實驗中選取2種不同長度的劃分窗口，針對不同大小的測試用例分別開展單進程和多進程運算，耗時結(jié)果如圖7所示。

可以看出，在CPU未滿載的情況下，對大窗口數(shù)據(jù)用多進程提取特征的加速比高于多進程對小窗口數(shù)據(jù)的加速比。因為當數(shù)據(jù)量較大時，每個CPU核心的利用率較高。算法花在CPU計算上的時間會遠遠多于花在進程切換和復(fù)制數(shù)據(jù)的時間。當窗口較小時，更多時間被浪費在進程切換和復(fù)制數(shù)據(jù)上，所以導(dǎo)致多進程加速比減小。

圖7 不同大小窗口的特征提取耗時Fig.7 Time consumption of feature extraction from different sizes of window

4.3 基于GPU的單一窗口向量特征計算

4.3.1 優(yōu)化方法

通過傅里葉變換提取頻譜特征是處理時序數(shù)據(jù)的常見操作。4.2節(jié)利用多進程在多核CPU上實現(xiàn)較大粒度的任務(wù)級并行。在單個進程內(nèi)，還可以利用GPU多線程實現(xiàn)較小粒度的線程級并行，進一步提升并行性能。

盡管Python的numpy庫中已經(jīng)提供了向量化的numpy.fft快速傅里葉變換函數(shù)，底層用C語言實現(xiàn)，其運算效率并不低，但是采用基于GPU的多線程并行還能帶來更高的加速比。

需要注意的是，3.4節(jié)簡單的計算是在Python中嵌入Cuda代碼來使用GPU，而此處則通過調(diào)用Cuda庫中已有的cufftExec方法來使用GPU。此時需要編寫基于Cuda語句的 .cu文件，導(dǎo)入cufft.h頭文件，調(diào)用cufftExec方法計算傅里葉變換，然后將Cuda文件編譯為動態(tài)鏈接庫 .so文件，再用Python的ctypes實現(xiàn)調(diào)用。另外 .so文件應(yīng)采用混合編譯，才能正常調(diào)用相關(guān)庫。

調(diào)用過程為：主機端先將原始程序運行過程中待處理的Python類型向量轉(zhuǎn)換為C類型，再將主機內(nèi)存中的C類型數(shù)據(jù)拷貝至設(shè)備端(GPU)全局內(nèi)存中。設(shè)備端讀取全局內(nèi)存中的數(shù)據(jù)后開始核心傅里葉變換計算，并將結(jié)果保存在全局內(nèi)存中。之后再將C類型結(jié)果數(shù)據(jù)從設(shè)備端拷貝回主機端，并轉(zhuǎn)換成Python類型。這個過程中應(yīng)特別注意Python與C類型數(shù)據(jù)間相互轉(zhuǎn)換的效率問題。

4.3.2 實驗結(jié)果

本實驗對比了針對不同長度窗口數(shù)據(jù)，在CPU上運行numpy.fft獲得傅里葉變換結(jié)果和在GPU上運行cufftExec方法獲得傅里葉變換結(jié)果的運行時間。耗時結(jié)果和加速比如圖8所示。其中，基于CPU的numpy.fft方法中調(diào)用一次fft所花費時間記為該方法的核心耗時。基于GPU的cufftExec方法中將不考慮數(shù)據(jù)類型轉(zhuǎn)換與數(shù)據(jù)拷貝過程的fft計算過程耗時記為核心耗時。算法整個完成時間記為總耗時?？梢钥闯?，數(shù)據(jù)量越大，GPU加速效果越明顯。

圖8 不同方法優(yōu)化后的特征計算效率Fig.8 Characteristic calculation efficiency after optimization by different methods

5 聚類過程優(yōu)化

5.1 TK-Means聚類算法

聚類是一種非常重要的數(shù)據(jù)挖掘方法，特別是對于時序大數(shù)據(jù)，絕大部分數(shù)據(jù)都沒有人工標簽，無法使用分類方法，而只能用聚類進行初步挖掘。

傳統(tǒng)的K-means聚類算法需要人工指定簇的數(shù)目，但是對于實際數(shù)據(jù)，最優(yōu)簇數(shù)往往很難事先確定。在TK-Means聚類算法中，簇的數(shù)不是一個固定值，而是一個范圍[min,max]。TK-Means聚類算法首先根據(jù)隨機初始閾值t得到初次聚類結(jié)果。如果初次聚類結(jié)果的簇數(shù)目不在控制范圍內(nèi)，則調(diào)整閾值t，重新進行初次聚類。若初次聚類結(jié)果的簇數(shù)目滿足要求，執(zhí)行正常的K-Means聚類過程，直到聚類結(jié)果穩(wěn)定。

TK-Means聚類算法步驟如下：

步驟1初次聚類。

步驟2如果初次聚類得到的簇數(shù)目在[min,max]范圍中，則執(zhí)行步驟3；否則調(diào)整閾值t并重復(fù)執(zhí)行步驟1。

步驟3對聚類結(jié)果進行調(diào)整，迭代次數(shù)遞增。

步驟4判斷聚類結(jié)果是否穩(wěn)定。如果聚類結(jié)果不穩(wěn)定，且迭代次數(shù)小于最大迭代閾值，則重復(fù)執(zhí)行步驟3；否則，聚類結(jié)束。

本文通過TK-Means聚類算法分別對第4節(jié)各特征向量進行聚類，并將得到的聚類結(jié)果表示成特征字符。對于聚類結(jié)果中點個數(shù)最多的簇，用a特征表示，即該簇中所有觀察窗口都是a特征。該特征的支持度，也就是該簇的支持度，定義為簇中窗口數(shù)目與該數(shù)據(jù)所有窗口數(shù)目之比。

5.2 基于GPU的特征向量聚類

5.2.1 優(yōu)化方法

通常的并行K-Means聚類算法只針對“計算所有點到簇心間距離”進行并行，而本文將算法迭代過程所有步驟都進行了并行優(yōu)化，并綜合運用GPU中的共享內(nèi)存、常量內(nèi)存、二維線程等措施優(yōu)化提升并行聚類過程時間性能。本文實現(xiàn)的并行聚類算法步驟如下：

步驟1初始化每個樣本的簇標簽。

步驟2迭代過程：①統(tǒng)計每一簇的樣本和；②統(tǒng)計每一簇的樣本個數(shù)；③計算每一簇的中心：樣本和/樣本個數(shù)；④計算每個樣本與每個簇中心的歐氏距離；⑤根據(jù)歐氏距離更新樣本的簇標簽。

上述過程從流程上看似乎與串行K-Means聚類算法沒有區(qū)別，但實際上其迭代過程是在GPU中并行完成的。本文以310 000個樣本(維度為10)聚成80類為例。在GPU中各個函數(shù)分解方式如下：

1) 初始化樣本簇標簽函數(shù)。由于只執(zhí)行一次，且計算量非常小，故不作為優(yōu)化的重點。線程塊維度：256；線程格維度：(310 000+256-1)/256=1 211。每個線程負責初始化一個樣本的簇標簽。

2) 統(tǒng)計每一簇的樣本和函數(shù)。如果每個簇的樣本都放在連續(xù)空間，那么此函數(shù)可使用類似規(guī)約求和的方式實現(xiàn)，且效率很高。但是此處每個簇是分散的，所以換一種方式。線程塊維度：(16,16)；線程格維度：((10+16-1)/16,(310 000+16-1)/16)=(1,19 375)。每個線程負責一個樣本中的一個數(shù)據(jù)。此處要使用原子操作，因為多個線程可能同時寫一個聚類中心的數(shù)據(jù)。也可以按下列方式劃分線程塊和線程格。線程塊維度：256；線程格維度：(310 000+256-1)/256=1 211。每個線程負責更新一個樣本，但是此時的效率通常不如之前高。

3) 統(tǒng)計每一簇的樣本個數(shù)函數(shù)。與初始化簇函數(shù)類似，計算量很少，每個負責處理一個樣本的計數(shù)，只有一個原子操作的加法，所以依舊采用一維的線程方式。線程塊維度：256；線程格維度：(310 000+256-1)/256=1 211。由于會有310 000個線程寫 80個位置，所以會存在許多沖突，此處對其進行優(yōu)化，開辟線程數(shù)變?yōu)椋壕€程塊維度：1 024；線程格維度：(310 000+1 024-1)/1 024=303。與此同時，在每個塊內(nèi)申請80個整數(shù)大小的共享內(nèi)存，先在塊內(nèi)進行統(tǒng)計，再寫到全局內(nèi)存，而不是像第一種方式那樣，直接寫全局內(nèi)存，這樣能避免很多沖突，獲得更高加速比。

4) 計算每一簇的聚類中心函數(shù)。此函數(shù)用于對每個樣本求和之后取平均操作，計算量極其少，不是優(yōu)化的重點，開辟二維線程塊。線程塊維度：(16,16)；線程格維度：((10+16-1)/16,(80+16-1)/16)=(1,5)。每個線程負責更新簇中心中的一個數(shù)。

5) 計算每個樣本與每個簇中心的歐氏距離函數(shù)。優(yōu)化重點，因 K-Means 聚類算法絕大部分計算量都集中在求每個樣本與每個簇中心的歐氏距離，故此處也開辟二維線程塊。線程塊維度：(16,16)；線程格維度：(80+16-1)/16,(310 000+16-1)/16)=(5,19 375)。每個線程負責計算一個樣本與一個簇中心的歐氏距離。

因函數(shù)計算過程中簇中心是不變的，可以考慮使用常量內(nèi)存將內(nèi)存訪問合并，進而隱藏內(nèi)存訪問延遲，在一定程度上提高讀取效率，以此加快對簇中心的訪問速度。

因此函數(shù)計算方式與矩陣乘法類似，可以考慮使用共享內(nèi)存，每個線程塊的任務(wù)是計算16個樣本與16個簇中心的距離，將16個樣本與16個簇中心的數(shù)據(jù)存到共享內(nèi)存中。

實驗中發(fā)現(xiàn)使用共享內(nèi)存的計算時間遠遠小于未使用共享內(nèi)存的計算時間。但是使用常量內(nèi)存的計算時間卻不一定比不使用常量內(nèi)存的時間還長，而根據(jù)簇心數(shù)據(jù)大小和設(shè)備硬件具體分析才行。不僅要看GPU最大可用常量內(nèi)存，還要看每個SM可供常量內(nèi)存使用的高速緩存是多大。

6) 樣本簇標簽更新函數(shù)。用于尋找每個樣本最近的簇中心，將當前樣本劃歸到該簇，可以開辟一維線程。線程塊維度：256；線程格維度：(310 000+256-1)/256=1 211。每個線程用于查找當前樣本對應(yīng)的與簇心間的最短距離(共80個)。因為求最小值與規(guī)約類似，而上述方式的每個樣本卻是完全串行的方式，所以對其進行優(yōu)化，開辟二維線程。線程塊維度：(16，16)；線程格維度：((1,(3 100+16-1)/16)=(1,19 375)。每個線程塊用于計算16個樣本的最小距離，也就是說用16個線程來完成原先1個線程的工作。最后，將規(guī)約后剩余2個元素的較小值作為最短距離。

5.2.2 實驗結(jié)果

實驗中選取含有不同樣本數(shù)目的用例，分別用原始串行程序和經(jīng)過并行與優(yōu)化后的程序進行聚類，得到兩者的耗時及加速比，如圖9所示?？梢钥吹?，在保證系統(tǒng)可用資源充足的情況下，聚類樣本數(shù)越大，并行化之后的加速比越高。

5.3 基于Cuda的自適應(yīng)聚類

5.3.1 優(yōu)化方法

TK-Means聚類算法的自適應(yīng)聚類依賴于距離閾值的調(diào)整，雖然可以采取二分法加快閾值收斂的過程，但卻不能被有效的并行化。

實驗中根據(jù)“簇大則分裂簇小則合并”的思想，提出一種K值自適應(yīng)調(diào)整算法，將K值估計問題轉(zhuǎn)化為單個簇中的樣本數(shù)估計問題，從而能夠在Cuda并行聚類程序中實現(xiàn)簇數(shù)的自適應(yīng)調(diào)整。算法流程如圖10所示。

分裂與合并過程為：設(shè)定一個簇中最小樣本數(shù)作為參考閾值，當前聚類迭代至最后一次時會將樣本個數(shù)小于該值的簇樣本單列，不參與之后的簇號更新，并可根據(jù)這些簇的有無或多少判定最大簇數(shù)是設(shè)置大了還是小了，根據(jù)設(shè)大了則合并設(shè)小了則分裂的原則，調(diào)整K值后開始下一次聚類。Cuda自適應(yīng)聚類算法則將上述過程簡化為只合并不分裂，即K從雙向收斂改為由最大值向最小值收斂。

圖9 并行優(yōu)化前后的聚類效率Fig.9 Clustering efficiency before and after parallel optimization

圖10 Cuda自適應(yīng)聚類流程Fig.10 Cuda adaptive clustering flowchart

5.3.2 實驗結(jié)果

實驗中選取的數(shù)據(jù)用例包含31個樣本，由5.1節(jié)的算法可知，最大聚類數(shù)Kmax應(yīng)設(shè)為6。此外，將最小聚類數(shù)設(shè)為3，類中的最少樣本數(shù)設(shè)為2，每次調(diào)K后重新聚類時的最大迭代次數(shù)設(shè)為100。圖11示例了改進后的Cuda聚類程序執(zhí)行過程和結(jié)果。

在示例中，自適應(yīng)聚類過程的初始K值等于6，經(jīng)過2次調(diào)整后變?yōu)?，最終的聚類結(jié)果中含有3個類，其中含有最小樣本個數(shù)的類為[2,2]，此時若K再調(diào)整為3，簇會進一步合并，那么所有簇中的樣本數(shù)都將超過參考閾值2，這與題設(shè)不符，故K調(diào)整至4便結(jié)束。由此可見，自適應(yīng)聚類中的整個調(diào)整過程完全符合預(yù)期。

圖11 Cuda自適應(yīng)聚類過程示例Fig.11 An example of Cuda adaptive clustering process

6 結(jié) 束 語

衛(wèi)星智能化應(yīng)用的不斷拓展對衛(wèi)星數(shù)據(jù)挖掘效率提出了更高要求。針對常規(guī)數(shù)據(jù)挖掘代碼進行并行優(yōu)化改進具有重要意義。本文以異變過程多類型特征分析過程為典型代表，針對窗口劃分與向量相似度計算、特征提取、傅里葉變換、聚類等常見數(shù)據(jù)挖掘操作，提出了多種優(yōu)化策略，大幅提高了算法執(zhí)行效率。綜合使用多種優(yōu)化策略說明對衛(wèi)星時序數(shù)據(jù)挖掘過程進行并行優(yōu)化改進是一個復(fù)雜的過程，需要全面考慮多個算法運算效率的影響。本文方法和思路有較好的通用性，可應(yīng)用于各種時序數(shù)據(jù)挖掘過程。當然，本文方法也并非完美無瑕，在CPU多進程+GPU多線程混合并行方面，以及基于Cuda流的多任務(wù)并行等方面，還有進一步改進空間。