近鄰搜索在多孔材料格點模型建模中的應用

劉任濤，陳　衛

LIU Rentao1，2,CHEN Wei1

1.中國科學院 計算機網絡信息中心 高性能計算技術與應用發展部，北京 100190

2.中國科學院大學，北京 100049

1.Department of High Performance Computing Technology and Application Development,Computer Network Information Center,ChineseAcademy of Sciences,Beijing 100190,China

2.University of ChineseAcademy of Sciences,Beijing 100049,China

1　引言

分治算法是計算機領域中一種廣為人知的算法，它的功能體現在于兩個方面：其一，通過多次的遞歸迭代自動地將復雜的問題分解成許許多多小問題；其二，將每個小問題的解沿著分解該問題的逆方向合并，直到得到原問題的解[1]。分治算法實現這兩個功能的具體方式在于將程序段多次遞歸迭代。近年來，由于能夠以相對簡單的程序段解決復雜的問題，分治算法的分而治之的思想被運用在諸多領域：例如在計算機科學領域利用分治算法研究聚類問題和群體智能算法[2-3]；在計算機圖形學領域利用分治策略劃分Delaunay三角網格[4-5]；在數值計算領域利用分治算法求矩陣特征值[6-7]；在系統科學領域利用分治算法研究種群的進化[8]。然而分治策略很少被應用于與計算物理和化學物理相關的計算機模擬當中。

多孔材料是一種具有非均勻內部結構的系統，在物理學[9-12]、化學[13]、材料科學[14]等領域引發了越來越多研究者的關注。將多孔材料體系離散為格點系統來進行計算機模擬是一種高效的計算機模擬方法，在基礎研究領域不斷涌現出新成果[15-16]。多孔材料的內部由兩部分組成：一部分是母體，被固體物質填充；另一部分是孔隙，可以容納流體（例如氣體和液體）。多孔材料格點模型（Porous Material with Lattice Model，PMLM）的計算機模擬分為初始化建模和體系演化兩個步驟，其中的初始化建模步驟需要把全體空間中的格點區分為兩類進行標記：一類格點標記為母體填充物（例如固體）；另一類格點標記為材料當中的孔隙。然而在實際研究過程中，對格點進行分類的過程并不是完全隨機的。為了逼近真實物理體系，通常母體區塊采用具有硬球勢或Lennard-Jones勢的球體[9，15]。原始建模方法不得不在固定每一個母體球的條件下遍歷地搜索全體空間的每一個格點，其計算量與系統體積的平方成正比。隨著模擬體系規模的增大，計算機模擬過程中的初始化建模步驟的計算量急劇增加，這是許多關于大規模多孔材料的研究無法承受的。

k近鄰搜索是一種用來取代全空間遍歷搜索的高效的搜索方式，在模式識別、機器學習、數據挖掘、時間序列分析等方面有諸多應用。利用k近鄰搜索可以預測數據樣本的分類，或以較小的代價尋找目標樣本[17-20]。受到k近鄰分類算法的啟發，本文把基于分治策略的近鄰搜索（Divide-Conquer based Nearest-Neighbor searching，DCNNs）方法應用到PMLM的建模當中，詳細介紹了在PMLM系統中應用DCNNs建模的方法。本文還在理論和實驗上對原始建模方法和DCNNs建模方法的運行效率作了比較。

2　基于分治策略的近鄰搜索方法理論

2.1　PMLM系統的基于分治策略近鄰搜索的優化建模方法

在一個三維格點模型中，格點的編號遵循：首先沿x方向布滿再向y方向延拓，布滿每一個xy平面后再向z方向延拓。因此，只要得到一個格點的序號就可以知道該格點在三維體系中的坐標(x,y,z)；反之，根據任意一個格點的坐標就可以知道該格點的序號。

多孔材料的顯著特征就是構型當中的母體物質與孔隙分散的占據全體空間。最常采用的硬球勢體系須將母體區域視為由一個個小球堆疊而成，孔隙區域占據剩余的空間。隨機選擇一個格點作為一個母體小球的球心A(x0,y0,z0)。為了標記出被球形母體占據的格點，最簡單的篩選辦法就是首先固定母體球的球心坐標，然后計算各個格點與球心之間的幾何距離：

球形母體區域的限制條件：如果某個格點到該球心之間的幾何距離di小于設定的母體球半徑R，則認為此格點在該母體球內部。

本文提出的DCNNs建模方法是以PMLM的一個母體區域的中心A作為根節點，通過分治的辦法將近鄰區域的格點劃分為不相交的區域，遞歸地向外部空間擴展以實現近鄰搜索。利用該DCNNs可以快速地標記出PMLM當中應當被母體占據的格點。

圖1顯示了用于PMLM建模過程中的DCNNs在不同類型格點位置上向近鄰格點擴展的狀態空間圖，圖中的每一個箭頭對應于一個方向上的擴展（體現在計算機程序中就是一個子程序）。算法對應的偽代碼如算法1～7所示。

圖1　近鄰搜索各個格點的狀態空間圖

偽代碼算法1表述了向6個方向擴展根節點的算法，狀態空間如圖1（a）所示。

偽代碼算法2和偽代碼算法3分別表述了向5個方向擴展z+節點或z-節點的算法，狀態空間如圖1（b）所示。

算法 2 UP（z，queue）

偽代碼算法4和偽代碼算法5分別表述了向3個方向擴展y+節點或y-節點的算法，狀態空間如圖1（c）所示。

DCNNs過程中的數據結構可以描繪成樹狀圖，如圖2所示，除了葉節點以外：

（1）根節點有6個子節點；

（2）每個z+或z-節點有5個子節點；

（3）每個y+或y-節點有3個子節點；

（4）每個x+或x-節點有1個子節點。

圖2　三維格點系統的DCNNs的搜索樹

從幾何空間上看，這樣的搜索機制類似于“碰壁”：設想在母體球外包裹了一層球殼，當計算機沿著一條搜索路徑向前搜索直到訪問到球殼上的格點時返回尋找下一條路徑繼續搜索。該DCNNs方式在搜索順序方面類似于深度優先搜索，搜索方向沿著搜索樹直到葉節點才返回，一言以蔽之就是“一路到底”。但區別在于：傳統的深度優先搜索的目的只是為了找到一個目標節點，因此一旦找到目標節點就會終止搜索；而本文介紹的搜索方法并沒有特定的目標節點，只有把搜索樹的所有葉節點全部搜索了一遍以后才會終止。

該基于分治策略的近鄰搜索方法吸收了數據分割的思想[20]，通過將近鄰空間格點劃分成合理的數據結構，保證了在同一次近鄰搜索當中既不會遺漏任何一個應當被搜索到的格點，又不會重復搜索任何一個格點。采用該DCNNs方法對PMLM系統進行建模，能夠準確地篩選出初始構型中應當被母體物質占據的格點。此外，只要改變程序中對于球形母體區域的限制條件，就可以將該DCNNs方法推廣用來對具有更復雜的母體形狀的PMLM系統進行建模。

2.2　建模時間復雜度的理論分析與比較

2.2.1原始方法建模的時間復雜度

對于PMLM系統，原始的算法是：首先隨機生成每一個母體區塊的球心的坐標，然后在固定該坐標的條件下遍歷搜索空間中的所有格點，計算各個格點的坐標與母體球心坐標的空間距離，從而篩選出位于該母體球內部的格點。

不妨設整個PMLM體系的格點總數為N；體系中被母體物質填充的格點個數與格點總數的比例是固定的，設為η(0＜η＜1)；同一體系中的各個母體球的半徑大小是相同的，設每個母體球占據n個格點。對于不同的體系，當母體球的半徑大小固定不變時，每個母體球中包含的格點個數n也是固定的常數（對于多孔材料n?N）。格點系統建模過程中的時間復雜度：

2.2.2DCNNs方法建模的時間復雜度

同樣，不妨設整個PMLM體系的格點總數為N；體系中被母體物質填充的格點個數與總的格點個數的比例為η(0＜η＜1)；每個母體球占據n個格點（n固定，且n?N）。對于同樣的PMLM系統，基于分治策略的近鄰搜索首先也是隨機生成每一個母體球的球心的坐標，然后在固定該坐標的條件下利用2.1節當中介紹的規則進行分類，在每一個母體球球心附近搜索到的格點個數，僅僅是母體球理論上包含格點個數n加上母體球外的“壁”包含的格點個數（其數量級一般不高于n具有的數量級）。格點系統建模過程中的時間雜度：

很明顯，在理論上，對于n?N的情形（例如多孔材料），DCNNs方法建模的時間復雜度更低，能夠節省大量計算量。

3PMLM系統氣液相變的計算機模擬

3.1　Gibbs系綜下基于密度泛函的格點模型平均場理論的熱力學函數

假定空間被離散化為三維的格點網絡，共有N個格點，第i個格點的粒子數密度為ρi，設第i個格點上的外場為Vi，第i、j格點之間的相互作用強度為εij，則整個體系的總勢能：

對于經典體系的熱力學正則系綜（固定的溫度、體積、總粒子數），全體粒子所處的密度狀態服從熱力學統計規律，任意一個給定的粒子數分布{ρi}出現的熱力學概率為：

其中，正則系綜配分函數：

因此，體系的Helmhotz自由能：

構造一個無外場無相互作用(Vi=0,εij=0)的參考體系，設參考體系的總能量為E?，則實際體系的能量E可以寫為參考體系能量與附加能量ΔE的加和：

在參考體系中，各個熱力學構型出現的幾率相等：

此時配分函數Z的物理意義可以形象地表述為M個粒子隨機分配給N個格點分配方式總個數，即服從超幾何分布：

對于M和N的值很大的情形（例如宏觀熱力學體系），運用Stirling近似ln(N!)≈NlnN-N ，則該參考體系的Helmhotz自由能：

其中，kB是Boltzmann常數，一個關于溫度和能量的物理常數。

由于該參考體系的粒子數分布可以視為宏觀上均勻的，設參考體系的粒子數密度ρ=M/N，因而該參考體系的Helmhotz自由能：

其中的密度ρ是無量綱的約化密度(0＜ρ＜1)。

引入平均場近似，設外場對每一個格點的作用強度為?i，則附加能量可以表示成：

3.2　密度迭代演化的計算機模擬方法

對于實際的物理系統，系統達到平衡時的密度分布使式（11）取得極小值，即滿足微分方程：

將格點之間的相互作用當作短程力處理，僅僅考慮最近鄰和次近鄰格點之間的相互作用[14]，并將表達式（11）代入式（12）即可得到達到熱力學平衡狀態下系統關于密度的代數方程[14]：

對于每一個格點，每次把舊的密度值代入上式右端，計算結果賦值給等號左端的新的密度值ρrhsi。為了防止密度取值跳動過大以致溢出（必須滿足ρ≤1），密度值與舊密度值進行線性組合，使得新得到的密度值僅僅在原有值的基礎上小范圍改動。因此，整個迭代的完整過程用等式（14）、（15）、（16）表示[14]：

等式（15）中的系數α(α?1)用來限制每圈迭代的速率；等式（16）用來標定體系的總粒子數守恒，得到每一輪迭代被更新的密度值。

反復迭代等式（14）、（15）、（16），使得系統的密度分布最終收斂。當某一圈迭代過后系統中密度值的改變量大小的平均值小于設定的閾值時，認為體系已經達到熱力學近平衡，終止迭代。

3.3　計算機模擬結果的比較

本文對采用周期性邊界條件的邊界長為50、65、80、100的立方體形格點模型分別進行計算機模擬。使用原始方法和DCNNs方法分別進行建模能夠構造相同的初始構型（流體的初始密度均勻ρ=0.28），并經過迭代演化達到同樣的熱力學平衡態。一個包含106個格點的三維體系的模擬結果，在x=0的截面的密度分布如圖3所示。圖右側標出了各種顏色代表的流體密度值（無量綱約化密度0＜ρ＜1），圖中空白區域是建模時建立的母體（流體密度為零），圖中密度稠密的團狀聚集區域是液相凝聚區，證明該體系已經發生了氣液相變。

圖3　計算機模擬多孔材料格點模型氣液相變達到熱力學近平衡時在x=0截面上的密度分布

盡管使用原始的建模方法與DCNNs建模方法能夠得到相同的物理結果，但是兩種方法的運行時間差異巨大。將程序在初始化和迭代過程中占用的運行時間進行比較，如圖4所示。隨著格點系統體積的增加，用DCNNs方法進行建模初始化占用的時間是線性增加的，系統遵循熱力學規律的迭代計算時間也近似于線性遞增，然而用原始的建模方法進行初始化占用的運行時間隨著系統的格點總數的二次方急劇增加。

圖4　不同規模體系模擬的各步驟的運行時間

當系統體積達到106個格點時，應用原始的建模方法與DCNNs建模方法進行計算機模擬的過程中各步驟占用時間的比例，如圖5所示。

圖5　N=106的體系在計算機模擬過程中運行時間的組成

不難發現，隨著體系規模的增大，若采用原始建模方法則初始化步驟會占據大部分的運行時間。加之系統演化過程可以很方便地并行加速，初始化步驟的運行時間在總運行時間當中占據的比例有可能會更大。因而對于采用原始方法建模的大規模PMLM系統，初始化步驟嚴重拖累了整個計算機模擬過程的運行時間，成為用計算機模擬大規模PMLM系統的技術瓶頸。應用DCNNs方法進行建模則只需占用幾乎可以忽略不計的運行時間來完成初始化，為更大規模PMLM系統的計算機模擬創造了條件。

4　結束語

本文主要研究了一種適用于PMLM系統進行計算機建模的優化方法。將分治思想與近鄰搜索策略相結合，提出的基于分治策略的近鄰搜索方法能將計算機建模PMLM系統的時間復雜度由O(N2)降低至O(N)。對PMLM系統氣液相變過程的計算機模擬，驗證了該方法在不影響物理結果的前提下的高效性。本文提出的優化方法突破了計算機模擬大規模PMLM系統的一個技術瓶頸。

參考文獻：

[1]Cormen T H，Leiserson C E，Rivest R L，et al.算法導論：Introduction to algorithms[M].殷建平，徐云，王剛，等譯.3版.北京：機械工業出版社，2013：16-17.

[2]王寶文，閻俊梅，劉文遠，等.基于分治法的高維大數據集模糊聚類算法[J].計算機工程，2007，33（24）：60-62.

7月中旬，我們形成了10個文件，經浦東開發領導小組、市政府常務會、市委常委會和市人大常委會相繼審議通過。

[3]李田來，劉方愛，王新華.基于分治策略的改進人工蜂群算法[J].控制與決策，2015，30（2）：316-320.

[4]余杰，呂品，鄭昌文.Delaunay三角網構建方法比較研究[J].中國圖象圖形學報，2010，15（8）：1158-1167.

[5]Cignonit P，Montanit C，Scopigno R.DeWall：A fast divide and conquer delaunay triangulation algorithm in Ed[J].Computer-Aided Design，1998，30（5）：333-341.

[6]Gansterer W N，Ward R C，Muller R P.An extension of the divide-and-conquer method for a cla-ss of symmetric block-tridiagonaleigenproblems[J].ACMTransactions on Mathematical Software，2002，28（1）：45-58.

[7]魏立峰，李曉梅.求解對稱帶狀廣義特征值問題的擴展分治算法[J].計算機研究與發展，2004，41（5）：861-867.

[8]王攀，萬君康，馮珊，等.一類基于分治原理的多種群協同進化算法[J].系統工程與電子技術，2004，26（11）：1687-1690.

[9]Brennan J K，Dong W.Phase transition of one-component fluids absorbed in random porous media：Monte Carlo simulations[J].The Journal of Chemical Physics，2002，116（20）：8948-8958.

[11]馬強，陳振乾.分形多孔材料雙尺度孔隙內氣體脫附擴散過程數值模擬[J].東南大學學報：自然科學版，2015，45（4）：728-733.

[12]劉高潔，郭照立，施保昌.多孔介質中流體流動及擴散的耦合格子 Boltzmann 模型[J].物理學報，2016，65（1）：282-290.

[13]馬強，陳俊，陳振乾.分形多孔介質傳熱傳質過程的格子Boltzmann模擬[J].化工學報，2014，65（S1）：180-187.

[14]崔靜潔，何文，廖世軍，等.多孔材料的孔結構表征及其分析[J].材料導報，2009，23（7）：82-86.

[15]Hughes A P，Thiele U，Archer A J.An introduction to inhomogeneousliquids，densityfunctionaltheory，and the wetting transition[J].American Journal of Physics，2014，82：1119-1129.

[16]Dickey J M，Paskin A.Computer simulation of the lattice dynamics of solids[J].Physical Review，1969，188（3）：1407-1418.

[17]祝繼華，尹俊，邗汶鋅，等.面向低維點集配準的高效最近鄰搜索法[J].模式識別與人工智能，2014，27（12）：1071-1077.

[18]譚駿祥，李少達，楊容浩.迭代最近點匹配算法的樹結構k近鄰搜索比較研究[J].測繪科學，2014，39（4）：152-155.

[19]衛煒，張麗艷，周來水.一種快速搜索海量數據集K-近鄰空間球算法[J].航空學報，2006，27（5）：944-948.

[20]肇瑩，劉紅星，王仲宇，等.最近鄰搜索用于分類問題的一種改進[J].南京大學學報：自然科學版，2009，45（4）：455-462.

[21]Beckmann N，Krigel H P，Schneider R，et al.The R*-tree：An efficient and robust access method for points and rectangles[C]//Hector G M，Jagadish H V.Proceedings of the 1990 ACM SIGMOD Conference on Management of Data.New York：ACM，1990：322-331.