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

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

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

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

1 填補算法技術(shù)路線

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

圖1 地質(zhì)真三維模型構(gòu)建時產(chǎn)生的“漏洞”

1.1 算法總體思路

本文以SGOG網(wǎng)格為例，基于網(wǎng)格編碼的拓?fù)渫评硖岢鼍暥葞茠咛钛a算法。在SGOG網(wǎng)格體系中，球面四叉樹編碼確定了網(wǎng)格的球面位置，徑向二叉樹編碼確定了網(wǎng)格的徑向位置。對位于同一球面位置的徑向內(nèi)部漏洞，可根據(jù)某位置上的最上、最下層網(wǎng)格直接進(jìn)行漏洞編碼填補。

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

1.2 算法關(guān)鍵技術(shù)

1.2.1 網(wǎng)格編碼緯度帶的劃分

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

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

根據(jù)上述方法，取出每個位于同一緯度帶上的所有內(nèi)部和邊界的SGOG球面四叉樹編碼，從邊界球面四叉樹編碼出發(fā)自東向西的順序?qū)ν痪暥葞系木幋a進(jìn)行鄰近計算(如圖4所示，方向m至n即從東向西)，若發(fā)現(xiàn)鄰近網(wǎng)格不存在，則對該網(wǎng)格進(jìn)行填補。修正方向編碼具有子三角形相對于其二級父三角形位置固定的特點，據(jù)此設(shè)計其鄰近搜索算法(另文著述)。

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

圖中虛線部分為模型區(qū)域擬定的邊界，藍(lán)色部分為存在的網(wǎng)格，則共有5種編碼填補情況：

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

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

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

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

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

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

遞歸上述方法完成對所有緯度帶內(nèi)編碼的填補即完成模型橫向模型填補，對橫向網(wǎng)格還需進(jìn)行高度位置的內(nèi)插，進(jìn)而完成整個模型的空間重構(gòu)。

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

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

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

2 實驗與分析

2.1 填補實例

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

圖6 航空港區(qū)模型填補前后比較圖

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

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

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

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

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

2.2 算法的時間復(fù)雜度測算

實驗環(huán)境：Intel(R) Xeon(R) Silver 4110 CPU @ 2.10GHz，RAM 64.0GB，GPU NVIDIA Quadro P2000，Win10專業(yè)版，使用MySQL數(shù)據(jù)庫對編碼數(shù)據(jù)的儲存和調(diào)用。實驗的基本原理為：選取12個存在漏洞的地質(zhì)模型，對這些模型進(jìn)行漏洞填補，記錄每個模型填補的漏洞個數(shù)與填補耗費的時間，繪制填補數(shù)量與時間關(guān)系曲線，如圖9所示?？梢钥闯觯W(wǎng)格填補數(shù)量(個)與工作時間(s)呈線性正相關(guān)，本文填補算法的時間復(fù)雜度為T(n) =O(f(n))。

圖9 填補算法效率趨勢圖

3 結(jié) 語

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