奮進號歷險記（中）

01-28

2013年冬，我搬家到亞省，是一個 NASA 項目，全球高解析度耕地製圖（GFSAD30），這件事情有別於此前工作的地方在於，區域性研究往往基於特定的先驗知識，就比如說，你想研究北京大興的西瓜地，那麼跑到大興去吃兩天西瓜，四下里看一看，也就可以了，甚至你再去保定跑兩趟，說不定就可以了解到華北的西瓜是個什麼情況，至少是大致上。但你肯定不會以為你手頭的數據可以把新疆的西瓜都代表進來。事實上，全球各地，還有你沒聽說過的很多地方都種著西瓜，這樣你不知道那裡的西瓜長成什麼樣了。

這是種繁瑣的說法，簡單的說法就是，為了了解數據的總體分布或特徵（全球的西瓜），我們需要一些數據探索手段。否則的話，就只能守著一些局部知識（大興的西瓜）過日子，而局部知識有多大的適用範圍，能代表總體分布里的百分之多少，也是需要了解的內容。

Landsat遙感衛星日夜監測地表長達四十年，已經建立起一套地表觀測記錄，這套數據去掉雲，留下晴好的部分，再使用一些數據處理技術進行填補，可以得到每個月的全球多波段影像。現在，假設我們已經知道了全球耕地的分布（以前 NASA有一些地表覆蓋類別的數據，並不準確但可以作為參考），怎麼樣組織我們的採樣方案，在複雜的地區多采些樣，在簡單的地區少採些樣，使得採樣的結果更接近於全球耕地的總體分布？

當然，實際上的採樣方案是一個多目標的的妥協結果，要考慮到具體採樣的難度、時間和經濟費用，但是，首先需要一個對全球耕地的分布刻畫：哪裡種類多，地塊破碎、分布；哪裡種類單一，地塊集中。這樣的話，對全球遙感數據做一次聚類分析，就是首要處理的問題，對於採樣，數據挖掘，分類器設計，能夠提供給我們很多的信息。

聚類分析是遙感圖像處理里的基本方法，問題在於如何處理如此大的數據（全球單景影像總量達 10TB）。一個簡單的辦法就是抽稀了後再去聚類，把 30m 的數據變成 1km再聚類，數據量少了，也能提供一些數據。但對於數據探索而言，這樣就失去了意義，因為你不能假定抽稀了以後總體分布不變。真實地表不具有這種連續性。比如，很多耕地都分布在河流沿岸，三角洲處，你一抽就可能把它們全跳過去了。

那麼一般的分布化，把這麼多數據切片，比如說切成 2GB，聚完類後再拼接起來，可不可以呢？

這是不行的。我先簡單介紹一下聚類演算法。其原理簡單來說就是，想像你有一堆星星需要聚成不同名字的星團，想要的聚類效果就是「星團內的星星都足夠近，星團間的星星都足夠遠」。首先你要確定所有的星星最後聚成幾類，然後挑選幾個點作為初始中心點，再然後依據預先定好的啟發式演算法（heuristic+algorithms）給數據點做迭代重置（iterative+relocation），直到最後到達「類內的點都足夠近，類間的點都足夠遠」的目標效果。也正是根據所謂的「啟發式演算法」，形成了k-means演算法及其變體包括k-medoids、k-modes、k-medians、kernel+k-means等演算法。另外，很多教程都告訴我們Partition-based+methods聚類多適用於中等體量的數據集，但我們也不知道「中等」到底有多「中」，所以不妨理解成，數據集越大，越有可能陷入局部最小。所以，我們還需要試探性地選擇初始類中心，即所謂的 kmean++。

具體說來，首先隨機選取k個宇宙中的點（或者k個星星）作為k個星團的質心，然後第一步對於每一個星星計算其到k個質心中每一個的距離，然後選取距離最近的那個星團作為所屬類，這樣經過第一步每一個星星都有了所屬的星團；第二步對於每一個星團，重新計算它的質心（對裡面所有的星星坐標求平均）。重複迭代第一步和第二步直到質心不變或者變化很小。也就是說，在這個過程中，所有的星星坐標（每個像素的多波段信息）都必須是全局可以見的，無法分布化，除非你設計出極為複雜的索引系統和聚類合併演算法。

奮進號，買賣上門了！

有了奮進號這種神器，心中自然狂野。為了把TB 級數據送入奮進號聚類，我需要把計算和I/O都分布化，否則數據的讀寫會成為新的瓶頸。我先是找到了parallel-netcdf庫，由西北大學和阿貢國家實驗室聯合開發的，是一種高速並行讀寫庫。pnetcdf 的開發者之一，Wei-keng Liao，在 pnetcdf 上實現了一個 k-means。他的學生 Jing Zhang 用 MPI 實際了並行 K-means，但是沒有跑過TB級數據。MPI 看起來不錯，在奮進號上也有優化後的編譯器。我下載了並行 K-means 的代碼，第一步先把地址長度從 32 改成 64（必須的嘛），然後——

然後程序就崩了。

西北大學的Wei-keng Liao教授溫文爾雅，給我很多適時的支援。在他的幫助下，我在奮進號上布署了一套海量數據的聚類器。用200個 CPU把2TB 的數據送入奮進號，2小時後，奮進號吐出了結果。