球體網(wǎng)格真三維地質(zhì)模型漏洞的緯度帶推掃填補(bǔ)算法*

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

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

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

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

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

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

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

1.1 算法總體思路

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

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

當(dāng)計(jì)算位于不同半球的SGOG球面四叉樹編碼的緯度帶時(shí)，SGOG球面四叉樹編碼所處的位置只會(huì)因?yàn)槟媳卑肭虻牟煌鴮?dǎo)致東西位置的改變，并不會(huì)改變編碼所在的南北位置，故不會(huì)影響緯度帶編號(hào)的計(jì)算，如圖3所示，所以此方法適用于全球SGOG球面四叉樹編碼。

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

1.2.2 基于緯度帶的SGOG編碼填補(bǔ)

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

圖4 基于編碼的填補(bǔ)情況示意圖

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

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

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

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

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

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

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

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

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

綜上，填補(bǔ)算法流程圖如圖5所示。

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

2 實(shí)驗(yàn)與分析

2.1 填補(bǔ)實(shí)例

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

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

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

圖7 長(zhǎng)沙市長(zhǎng)沙縣模型填補(bǔ)前后比較圖

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

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

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

2.2 算法的時(shí)間復(fù)雜度測(cè)算

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

圖9 填補(bǔ)算法效率趨勢(shì)圖

3 結(jié) 語(yǔ)

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