球體網格真三維地質模型漏洞的緯度帶推掃填補算法*

王金鑫，曹澤寧，陳藝航，秦子龍，石 焱

(1.鄭州大學 地球科學與技術學院，河南 鄭州 450001；2.鄭州大學 水利科學與工程學院，河南 鄭州 450001)

地球剖分網格(Earth Tessellation Grid, ETG)是一種基于(橢)球體的對整個地球重力場進行遞歸剖分而形成的多分辨率的離散網格地球系統空間參考模型[1]，包括傳統面向地球系統大型科學計算的地球物理和地球系統網格[2]以及面向時空大數據整合、集成、建模與應用的數字地球平臺網格[3]兩種。數字地球平臺網格又包括全球離散網格(Discrete Global Grid System, DGGS)[4]和地球系統空間網格(Earth System Spatial Grid，ESSG)[5]兩種。球體測地線八叉樹網格(Sphere Geodesic Octree Grid，SGOG)[6-7]是ESSG模型的一種，首先進行球面大圓弧QTM四叉樹剖分，然后進行徑向二叉樹剖分，再完成球體的三維網格分割。具有剖分規則簡明、網格形體簡單、排列規整、幾何特征明晰的特點[6-7]。

與二維圖件相比較，三維地質成果圖件具有諸多優點，如基礎資料利用的充分性與靈活性、成果的綜合性、應用的直觀性等[8]，并可動態地顯示和查詢地層模型的幾何特征[9]，對地學工作者和地質工程師進行地質分析極為重要[10-14]。如何根據已有的離散地層數據，構建出完整的三維柵格地質模型，是采用鉆孔數據建立三維地質模型的關鍵和難點。本文針對球體網格真三維地質建模存在的漏洞(橫向與徑向漏洞)修補問題，基于多邊形填充原理和SGOG編碼鄰近搜索方法，設計了修補算法，取得了良好的效果。

1 填補算法技術路線

傳統利用鉆孔數據和體素(元)構建三維地質模型的方法，一般是找到鉆孔數據的空間外包圍盒(常為規則的幾何形體)，然后進行空間剖分，邏輯上屬于由外向內的方法。這種方法雖然不存在漏洞問題，但不利于復雜地質體的模型建立。球體網格概念的提出，基于網格的規則剖分和編碼機制，可直接利用(或經適當加密的)離散采樣點數據進行三維構模，邏輯上屬于由內向外的方法。該方法十分有利于復雜地質體的三維建模，但由于采樣點數量的稀疏性，常常會出現由于數據密度低而產生模型“漏洞”問題，包括橫向表面漏洞(如圖1(a)所示)和徑向內部漏洞(如圖1(b)所示)。

圖1 地質真三維模型構建時產生的“漏洞”

1.1 算法總體思路

本文以SGOG網格為例，基于網格編碼的拓撲推理提出緯度帶推掃填補算法。在SGOG網格體系中，球面四叉樹編碼確定了網格的球面位置，徑向二叉樹編碼確定了網格的徑向位置。對位于同一球面位置的徑向內部漏洞，可根據某位置上的最上、最下層網格直接進行漏洞編碼填補。

橫向上下表面漏洞的編碼填補是本算法的關鍵。其主要思路為：首先，根據實際的應用需求，確定地質體三維建模的球面剖分層次(包括徑向剖分層次)；其次，依據SGOG網格球面四叉樹編碼，進行網格的緯度帶劃分；再次，根據空間網格的球面拓撲關系以及填補規則，以地層鉆孔數據圈定的空間范圍作邊界約束條件，遍歷進行橫向漏洞的編碼填補；然后，對新增的橫向球面網格進行徑向上下面的高度位置內插并對其進行徑向填充，進而完成整個地層模型的完整構建。對于地質體中多地層數據，重復進行上述流程，完成每層模型的構建，基于實際地理坐標進行集成即可。

1.2 算法關鍵技術

1.2.1 網格編碼緯度帶的劃分

在一定的球面剖分精度下，模型上的球面網格所在的緯度位置(并非真正意義上的緯度)是該球面三角形編碼的函數，而緯度位置可以用緯度帶作為行數據表示。需要說明的是，SGOG網格的邊為測地線大圓弧，與緯線小圓弧不重合。這里的緯度帶僅作為填補的順序導向，不是嚴格意義上的網格范圍，填補通過編碼實現，所以SGOG網格的邊不影響算法的實現。

對于給定層級SGOG剖分，4個球面三角形(屬于某個上層父三角形)被劃分到2條緯度帶中，分別記為第‘0’緯度帶和第‘1’緯度帶，如圖2所示。那么當該級子剖分單元的編碼(這里采用修正方向編碼方法[7])為‘1’時，返回當前層級的緯度帶標記值index=1；當該級子剖分單元的編碼為‘0’、‘2’、‘3’時，返回當前層級的緯度帶標記值index=0，即

圖2 上、下三角形的緯度帶劃分

(1)

式中，index為計算當前層級的緯度帶標記值；code[n]為SGOG球面編碼序列；i為當前計算編碼的剖分子級(從0開始計)。

這里的緯度帶是根據球面剖分層次從低緯度到高緯度按照升序的方式劃分的，記低緯度帶的編碼為‘0’，高維度帶編碼為‘1’(北半球的情況)。當上一層級對應的球面三角為上三角時(即上一層級編碼為‘1’)，編碼為‘1’的子三角處于高緯度位置(如圖2左圖所示)；當上一層級對應的球面三角為下三角時(即上一層級編碼為‘0’)，編碼為‘1’的子三角處于低緯度位置(如圖2右圖所示)；故此時存在緯度帶劃分的正反序問題，必須進行相應的調整。當計算的上一層級球面三角為下三角時，記編碼為‘1’時返回的標記值為‘0’，編碼為‘0’、‘2’、‘3’時返回的標記值為‘1’，故對式(1)修改得：

(2)

通過計算SGOG網格每一層級編碼的緯度帶標記值，可得到一條長度(位數)與網格編碼相同的緯度帶標記值序列index[n]，對該序列每一位數按照如下公式計算，即可得到該SGOG球面編碼所在的緯度帶值，即

(3)

式中，Nlat為計算得到的緯度帶編號；index[n]為緯度帶標記值序列；n為返回的標記值長度；i為當前標記值位置(從0開始計)。

當計算位于不同半球的SGOG球面四叉樹編碼的緯度帶時，SGOG球面四叉樹編碼所處的位置只會因為南北半球的不同而導致東西位置的改變，并不會改變編碼所在的南北位置，故不會影響緯度帶編號的計算，如圖3所示，所以此方法適用于全球SGOG球面四叉樹編碼。

圖3 南北半球在2級剖分層次下的各編碼緯度帶位置示意圖

1.2.2 基于緯度帶的SGOG編碼填補

根據上述方法，取出每個位于同一緯度帶上的所有內部和邊界的SGOG球面四叉樹編碼，從邊界球面四叉樹編碼出發自東向西的順序對同一緯度帶上的編碼進行鄰近計算(如圖4所示，方向m至n即從東向西)，若發現鄰近網格不存在，則對該網格進行填補。修正方向編碼具有子三角形相對于其二級父三角形位置固定的特點，據此設計其鄰近搜索算法(另文著述)。

圖4 基于編碼的填補情況示意圖

圖中虛線部分為模型區域擬定的邊界，藍色部分為存在的網格，則共有5種編碼填補情況：

情況1：編碼a左、右鄰近編碼均存在，則不需要進行編碼填補；

情況2：編碼b僅存在左鄰近或僅存在右鄰近或不存在左右鄰近，對缺少的編碼進行編碼填補，并對填補后的編碼繼續做左(右)鄰近判斷，直到觸碰到邊界為止；

情況3：編碼c位于邊界上且存在另一邊的鄰近編碼，則不需要進行編碼填補；

情況4：編碼d位于邊界上但不存在另一邊的鄰近編碼，對缺少的編碼進行編碼填補，并對填補后的編碼繼續做左(右)鄰近判斷，直到觸碰到邊界為止；

情況5：編碼e已超出邊界范圍，刪除超出邊界范圍的編碼。

根據以上規則，對同一緯度帶的編碼填補工作流程如下：(1)記在同一緯度帶下所有存在的網格編碼集合為D，根據SGOG網格鄰近關系將集合D中各元素Di進行位置順序的排列；(2)獲取在該緯度帶上邊界的位置信息，規定填補方向(由東向西或由西向東)，以由東向西為例，在填補時每個緯度帶上可能存在多個邊界(總數為偶數)，則根據從東向西的順序分別對邊界賦予“進”與“出”的標簽；(3)對“進-出”標簽內的網格編碼Di根據圖4中的前4種情況進行填補，對“出-進”標簽內的網格編碼Di根據圖4中的第五種情況，認為其已超出給定的邊界范圍外則刪去。

遞歸上述方法完成對所有緯度帶內編碼的填補即完成模型橫向模型填補，對橫向網格還需進行高度位置的內插，進而完成整個模型的空間重構。

另外，對于存在空洞(真實存在的洞)的復雜地質體，必須首先給出空洞的數據邊界(空洞部分的地層數據)，然后將其與球體網格匹配(柵格化)，得出內邊界后再進行填補。

綜上，填補算法流程圖如圖5所示。

圖5 基于SGOG方向修正編碼填補算法流程圖

2 實驗與分析

2.1 填補實例

為了驗證本文填補算法的可行性，選取鄭州市航空港地區某地質層三維構模作為研究對象，構建橫向17層剖分層次、徑向25層SGOG剖分層次的地質模型(建模前鉆孔數據經過適當加密)，模型填補效果如圖6所示。共填補球面橫向漏洞5 078個，其約占研究區域總面積的2.5%(研究區域總面積約643.5 km2，漏洞面積約為16.4 km2)。

圖6 航空港區模型填補前后比較圖

由于上述模型的邊界較為簡單，故截取地形邊界較為復雜的長沙市長沙縣(山地丘陵地形)的DEM數據(來自于共享網站：http：//srtm.csi.cgiar.org，將DEM數據隨機刪去20%來模擬漏洞)作為地層數據進行地表真三維建模來驗證算法有效性(模型剖分參數為：SGOG橫向16層，徑向25層)。研究區域表層橫向網格總數量129 830個，其中人為刪除(模擬漏洞)數量21 638個，使用本文算法共填補漏洞網格21 312個，填補產生誤差小于2%，滿足三維模型的可視化及空間分析等需求，填補效果如圖7所示。

圖7 長沙市長沙縣模型填補前后比較圖

根據圖6、圖7對比看出，無論地形復雜與否，本算法均達到了良好的修補效果。

對航空港地區的地質模型，給定2個空洞的邊界，其修補效果如圖8所示，可以看出本文算法在避開固有空洞的前提下對模型中的漏洞填補效果良好，并未出現對固有空洞造成錯誤填補的情況。

圖8 存在固有空洞地質模型的漏洞填補示意圖

2.2 算法的時間復雜度測算

實驗環境：Intel(R) Xeon(R) Silver 4110 CPU @ 2.10GHz，RAM 64.0GB，GPU NVIDIA Quadro P2000，Win10專業版，使用MySQL數據庫對編碼數據的儲存和調用。實驗的基本原理為：選取12個存在漏洞的地質模型，對這些模型進行漏洞填補，記錄每個模型填補的漏洞個數與填補耗費的時間，繪制填補數量與時間關系曲線，如圖9所示。可以看出，網格填補數量(個)與工作時間(s)呈線性正相關，本文填補算法的時間復雜度為T(n) =O(f(n))。

圖9 填補算法效率趨勢圖

3 結 語

本文以SGOG網格為例，就利用鉆孔數據生成真三維地層模型的漏洞填補問題，提出了緯度帶推掃填補算法。該算法充分利用球體網格的剖分規則及其空間拓撲關系，以網格編碼為核心，實現地質實體的遍歷填補。研究表明：(1)利用規則球體網格建立真三維地質實體模型是一種可行的方法，尤其是對大區域和全球尺度地質空間，具有較大的精度優勢(考慮地球曲率，無投影變形)；(2)SGOG網格兼具TIN與Grid的特點，在真三維地表(地殼)可視化表達與空間分析中具有較大應用潛力；(3)本文提出的方法對所有球體網格具有普適性，填補效果良好。本研究對新一代數字地球平臺建設及時空大數據管理與應用具有參考價值。