一種基於SGOG瓦塊的大區域真三維地理場景自適應構建方法與流程

2023-06-11 22:21:16

本發明屬於地理信息科學技術領域，具體涉及一種基於SGOG瓦塊的大區域真三維地理場景自適應構建方法。

背景技術：

1.研究背景

地理信息科學中，三維地理場景的構建方法一般有兩種：一種是面元模型加紋理的方法，往往可以達到較好的地形和建築可視化效果，在3DGIS領域應用廣泛；另一種是基於體素模型的建模，常用來表達地質體垂向的多個界面，並方便水平與垂直剖面的製作，在三維地學模擬信息系統中應用廣泛。面模型不能很好表達洞穴、拱橋和懸崖等複雜的地理實體，難以實現地形體數據表達和實時編輯，基本上不可能實現地理空間過程的模擬；而基於體素的模型可以實現虛擬地理場景的快速構建與多解析度實時拆分，實現體素模型、表面模型的一體化表達、多維地理對象及地理過程的一體化表達，地理空間數據可定製及自適應表達，面向複雜地學分析，能夠較好地表示空間粒子之間的相互作用關係，適用於表達地理空間過程的動力學機制，進而實現對多尺度地理現象與過程的模擬、分析與預測。但上述傳統方法，無論是面元模型還是體素模型，都是面向局部的、小範圍的歐氏投影模型。數字地球背景下，傳統的投影模型在面向大尺度和全球性問題的應用時，存在著嚴重的局限。

在大區域或全球尺度，地球曲率不容忽略。因而，投影模型不再適用，必須基於顧及地球曲率的全球離散網格模型(Discrete Global Grid Systems，DGGs)構建三維場景,利用退化四叉樹網格實現了全球多尺度地形無縫可視化。其實，以Google為代表的50餘種第一代數字地球平臺，都是利用球面剖分網格(瓦片或面片)構建全球DEM框架。但以上所有基於DGGs的三維模型歸根結底仍屬於表面模型的範疇，不涉及地表上下的空間。近些年來，在地理信息科學領域，球體剖分網格，即地球系統空間網格(Earth System Spatial Grid，ESSG)模型被提出,利用球體退化八叉樹網格(Spheriod Degenerated Octree Grid,SDOG)建立了青藏高原巖石圈的三維建模。基於ESSG瓦塊的地理空間真三維表達方興未艾，亟待發展。需要說明的是，地球系統空間網格與傳統的地球物理和地球系統網格有很大區別。而地球系統網絡是一種科學計算的工具，基本不涉及或不強調地球場域的統一定位框架及一體化空間數據的組織、管理與應用。

2.SGOG剖分方案

球體大圓弧QTM八叉樹網格(Sphere Geodesic Octree Grid，SGOG)是球體三維網格剖分理論。QTM(Quaternary Triangular Mesh)指四元三角網。簡單的理解就是把一個三角形三邊的中點用直線連接，就可以把一個三角形劃分為四個子三角形，同理，四個子三角形又可以劃分為16個孫三角形……，以此類推。設定剖分層次，就可以構成一張三角形格網。顯然，四元三角網是一種面剖分方式，可以是平面或曲面。球體QTM八叉樹網格剖分方法可分為兩步：把地球看做一個球體，首先進行球面離散剖分，然後再進行球體徑向剖分。具體為：第一步，球面剖分。首先，以0°～180°首子午圈和與之垂直的東西經90°子午圈和赤道把球面分為等積的8個球面三角形。接著，對每個三角形分別取三邊弧的中點，並用大圓弧連接，形成四個二級的球面三角形；然後依據此方法依次對此後各級的球面三角形進行遞歸剖分，直到滿足應用為止(從八分之一球體起算，遞歸的層次數為n，下同)。第二步，徑向剖分；設球的半徑為R，從網格的每個節點向球心引直線，然後以R/2n將半徑等分，在各分層處以球面進行剖分(以球心為原點，以k*R/2n為半徑畫球，k＝1,2,…2n-1)。這樣，就可以把球分割為整齊一致的球面三稜台(上下底面為球面，側面為平面)，其中球心處為球面三稜錐。以稜台/錐的幾何中心為網格參考點。把球理解為包括大氣層在內的地球，可實現對整個重力場的剖分，其三維可視化效果如圖1、2所示。

球體SGOG剖分還包括徑向不等長八叉樹剖分，由於本發明不涉及不等長剖分，所以這裡不再介紹。

3.SGOG編碼模型

已有研究表明：在重力場中，從地心到磁層，16倍的地球半徑就可以滿足包含地球系統所有圈層、基本覆蓋人類活動空間的要求。考慮到網格定位及計算的方便，SGOG網格採用以下編碼模型：圈層碼(十六進位)+八分體標識碼(八進位)+球面位置編碼(四進位)+徑向深度碼(二進位)。這裡的圈層，不是嚴格意義上的地球圈層，而是以整數倍地球半徑表示的從地心和磁層的距離。八分體標識碼表示網格所在的八分體的位置，從0°開始，沿經度增加的方向，按每90°為一個卦限，北半球分為0～3、南半球分為4～7共八個卦限。球面QTM四叉樹網格已有很多種成熟的編碼方案，如固定方向編碼、ZOT編碼和LS編碼等。本發明採用固定方向編碼，徑向深度採用二叉樹編碼(靠近球心的位置取0，如圖4所示)。

球面QTM四叉樹固定方向編碼的原理：

如圖3所示，該法將八面體按等邊三角形投影(Equal-Triangles Projection，ETP投影)，首先把經緯度通過ETP投影轉換為x，y值；然後比較位置點(x，y)到四個三角形中心點的距離；選擇距離最小的三角形，記錄其地址碼；依此類推，直到劃分到一定層次或位置點到三角形中心的距離小於一定的值。

徑向二叉樹編碼原理：

如前所述，球體八叉樹剖分是球面四叉樹和球徑二叉樹剖分的組合，用球面四叉樹編碼加上徑向二叉樹編碼即可確定任一八叉樹網格的空間位置(一個完整的八叉樹網格編碼還包括首位0,1,2，……7的八分體識別碼)。二叉樹的基本編碼方法如下：把每次剖分的格按從球心到球面的次序轉化為與層次相對應的二進位位編碼即可。如圖4所示：第一次剖分的編碼為球心處0、1，第二次剖分為00、01、10、11，依此類推。

4.SGOG解碼方案——編碼到經緯度和空間直角坐標的相互換算

1)編碼到經緯度坐標的相互換算

所謂離散網格編碼到經緯度的相互轉換是指：已知網格的編碼，求網格體元中心點的經緯度；或者反過來，給定某點的經緯度求該點所在某層次(或滿足某種精度)網格的體元編碼。本發明以一倍的地球半徑為例，並將網格的定位分為球面四叉樹經緯度求解和徑向二叉樹深度(離地心的距離)求解兩個步驟。

依據相關研究實現方向四叉樹編碼到經緯度相互轉換。該法的核心是將八面體按等邊三角形投影，把經緯度轉換為笛卡爾直角坐標，並以其為橋梁，實現從地址碼到經緯度的相互變換。其主要數學模型如式(1)和(2)所示。

式中，λ、φ為經緯度；x、y為ETP投影坐標。

徑向二叉樹編碼的實質就是以某層次的網格稜長去度量地球的半徑。徑向二叉樹的編解碼原理：設某一點到球心的距離為d，剖分的層次為n，用d除以地球的平均半徑得商e，把e轉化為以2n為分母的分數，取其分子的整數部分加上1(如果為整數則不加1)，然後再將其化為n位二進位碼，即得到其二叉樹碼；反之，已知某二叉樹碼(對於一個完整的八叉樹碼，需根據其規則截取二叉樹碼)，設其位數為n，把其轉化為十進位數，再乘以地球的平均半徑與2n的商就得到該編碼所對應的點到球心的距離。由此可知，二叉樹的編碼沒有誤差，而解碼的精度在一個網格範圍內。變長八叉樹的編解碼方法與八叉樹原理上是一樣的，其精度也相仿，但剖分的比例不同。將變長比例代入相應的解碼公式即可。

2)編碼到空間直角坐標的相互換算

這裡網格編碼到空間直角坐標的轉換是指由網格的編碼求網格單元各頂點的空間直角坐標，進而可以求得SGOG瓦塊內外三角形的重心坐標，將二者取平均值，可以得到SGOG瓦塊的幾何中心的坐標；而空間直角坐標到網格編碼的轉換是指已知某點的空間直角坐標和剖分層次，求其所在網格單元的編碼。由於SGOG採用了大圓弧中分的剖分規則，所以整個網格體系與空間直角坐標對應十分整齊。下層網格瓦塊新增頂點與上層瓦塊的頂點之間存在著簡單的中分關係。通過求出弦的中點，將其投影到相應的剖分球面上即可。

本發明在算法設計時，對上面的編碼模型做了一些變動：將整個碼用標示符A和B分為三個碼段。A之前的為圈層碼，由0～n個1構成。若沒有1，則表示是固體地球本身，即1倍的地球半徑；每多一個1，球的半徑就乘以2，以此類推。A和B之間的二叉樹碼表示網格體元瓦塊的徑向位置。B之後的第一位八進位碼為瓦塊所在的八分體標識碼，其餘的四進位碼表示瓦塊在球面的橫向位置。根據以上編碼和剖分規則，給定初始條件，設計八叉樹編碼與空間直角坐標轉換的算法思路如圖5所示。圖5為編碼到空間直角坐標，圖6為空間直角坐標到編碼。

5.SGOG網格可視化繪製技術

SGOG網格體系包括球面三稜台和球面三稜錐兩種瓦塊。瓦塊是離散網格的基本單位，也是地理空間建模的基礎。因而，對瓦塊結構體的設計與實現是首要的關鍵技術。本發明中的瓦塊繪製程序基於C++平臺的應用程式接口OSG開發，其結構體如下：

由結構體繪製線框模型共需3步：

(1)外層3個頂點兩兩之間繪製大圓弧；

(2)內層3個頂點兩兩之間繪製大圓弧；

(3)內外層對應頂點之間繪製直線，繪製時用多段的弦來逼近弧。

技術實現要素：

本發明所要解決的技術問題是大區域真三維地理場景構建問題，利用IDL語言處理DEM模型源數據，以OSG作為圖形引擎，以VS2010為平臺，以標準C++作為開發語言，構建實驗系統。基於球體大圓弧QTM八叉樹網格剖分理論，利用全球共享數據，建立了中國大陸地區的真三維地形框架，實現了大尺度地理場景無縫集成、自適應建模與可視化，為真三維地理空間計算打下堅實基礎。

本發明所採用的技術方案為：一種基於SGOG瓦塊的大區域真三維地理場景自適應構建方法，包括如下步驟：

S1，SGOG瓦塊網格點坐標計算與存儲；

據SGOG剖分的大圓弧中分規則，由大圓弧的兩個端點坐標求其弦的中點坐標，然後將其沿直線延伸到球面上，得到該網格點的坐標；

S1a，每層網格的坐標存儲三個TXT文件，分別存放X、Y、Z坐標；

S1b，按方向編碼順序存儲第一層網格中N個三角形的頂點坐標；

S1c，按照第一層父三角形存儲的順序，存儲第二層網格中子三角形，每四個子三角形仍按方形編碼排序，依次類推，存儲下一層自三角形，直至存儲完成；

S2，進行DEM高程匹配與計算；

S2a，在IDL平臺下讀取存儲的SGOG最外層網格點X、Y坐標，並輸入地球半徑，換算為經緯度；

S2b，遍歷DEM影像數據，根據影像的最大、最小行列號，換算出影像的經緯範圍，判斷網格點是否落在當前影像中；

S2c，網格點落在當前影像則讀入影像高程值，並把落在影像內部的網格點經緯度換算為所處當前影像的行列號，根據行列號讀出像元值賦給網格點高程，如果網格點不在當前影像範圍內則判斷下一幅影像；

S3，繪製瓦塊構建真三維地形框架；

基於VS和OSG平臺，在獲取最外層網格點高程的基礎上，按照所有可視範圍的SGOG網格點坐標，按照網格可視化繪製方法繪製網格，將高程值適當放大，增強視覺效果，進而實現大區域真三維地形框架的建模與可視化；

S4，自適應建模；

根據地形起伏度大小和建模需求，設定相應的閾值，靈活採用不同剖分層次的網格進行建模可視化。以某一層為基礎層次，起伏度大於閾值的地方，進一步細分；起伏度小於閾值的地方，不再細分，直到滿足精度要求為止。這樣既有效減少了數據量，又保留了重要的地形特徵，進而達到較好的可視化效果。

S5，進行分層設色和紋理渲染；

根據網格點的不同高程值為每個網格點賦上相應的顏色值，利用OSG渲染引擎自帶的函數，隨著瓦塊編碼的遍歷而繪製瓦塊的各個表面並進行顏色渲染，可得到DEM場景暈渲圖。當然，也可以在網格表面粘貼遙感影像紋理，構建虛擬地理場景。

進一步，所述的閾值可以採用某層網格的起伏度平均值或其相應的比例。

進一步，所述的高程插值可以採用反距離加權法、線性內插、雙線性內插或最鄰近像元法。

本發明產生的有益效果是：

1、基於地表大數據的一體化無縫組織，實現了考慮地球曲率的、大區域地殼真三維可視化建模，克服了基於投影空間模型的變形和裂縫等缺陷，恢復了地理空間的自然流形性質。

2、提出了基於地表起伏度閾值的自適應可視化建模方法，既減少了數據量，又保留了重要地形特徵細節，保證了建模與可視化的效果。

3、本發明的方法為建立真三維數字地球平臺打下了堅實基礎，並可以應用於天地一體化的空間數據的組織、管理與應用。

附圖說明

圖1為以八分之一球體為例的球體QTM八叉樹剖分效果圖；

圖2為整球剖分下的球體QTM八叉樹第2層剖分效果圖；

圖3為ETP投影算法原理與編碼方案；

圖4為徑向二叉樹編碼方案；

圖5為SGOG編碼轉換空間直角坐標的流程圖；

圖6為空間直角坐標轉換SGOG編碼的流程圖；

圖7為真三維地理場景構建技術路線；

圖8為高程匹配方法流程圖；

圖9為大區域真三維地形框架的第4層瓦塊和第9層瓦片；

圖10為大區域真三維地形框架中第6、7、8、9層的自適應可視化；

圖11為圖10的側視圖；

圖12為自適應真三維渲染圖；

圖13為圖12的局部放大圖；

圖14為真三維地理場景效果圖；

圖15為圖14的局部放大圖。

具體實施方式

以下為本發明的一種實施例。

考慮到單機的計算與存儲能力，本實施例選擇中國大陸地球為研究區域，以SGOG第9層剖分瓦塊為建模基礎，該網格球面邊長約為19.5km，徑向稜長約12.4km。從共享網站下載覆蓋中國大陸地區的1144幅90m解析度的DEM數據。基本參數：投影UTM/WGS-84，GeoTIF格式，3601×3601像元，高程為相對於WGS-84橢球的大地高。

本發明利用IDL語言處理DEM模型源數據，以OSG作為圖形引擎，以VS2010為平臺，以標準C++作為開發語言，構建實驗系統。

如圖7所示，首先進行SGOG瓦塊網格點坐標計算與換算，通過經緯度坐標與DEM網格匹配，接著進行DEM高程匹配與計算，然後繪製瓦塊構建真三維地形框架，最後進行分層設色和紋理渲染。

對SGOG瓦塊網格點坐標計算與存儲

SGOG採用大圓弧中分遞歸剖分，網格上下層之間存在著整齊的對應關係。網格點坐標計算的方法是：首先由大圓弧的兩個端點坐標，求其弦的中點坐標，然後將其投影即沿徑向直線延伸到球面上即可。由於在現有的研究中已經實現網格編碼、網格頂點包括幾何中心點的直角坐標及其經緯度三者之間的相互換算。

不同層次的大圓弧QTM網格之間有著明確的包含與被包含關係，按照一定的順序存儲各層次的網格點坐標可以方便的查找與某一網格有著包含或被包含關係的網格。由SGOG剖分思想可知，每個三角網遞歸剖分1次形成4個三角形，所以可推出第N層網格的一個三角形網格包含的第N+I層網格的三角形網格數量為4(I-1)，網格點數為三角形數量的3倍。第1層三角形無父三角形，末層三角形無子三角形。這樣除第1層和末層之外，每個父三角形有4個子三角形。從第1個子三角形開始，每4個子三角形有1個公共父三角形。從坐標來看，父三角形的3個頂點坐標與4個子三角形12個頂點相關聯，存儲時保留這種對應關係，方便後續數據操作時對父子三角形的統一處理。

最外層網格具體存儲方法為：

1)每層網格的坐標存儲3個TXT文件，分別存放X、Y、Z坐標；

2)按方向編碼順序存儲第1層網格中N個三角形的頂點坐標；

3)按照第1層父三角形存儲的順序，存儲第2層網格中子三角形，每4個子三角形仍按方向編碼排序。然後依次類推，存儲下一層子三角形，直至存儲完成。

在進行真三維DEM可視化和自適應可視化操作時，徑向各層網格的存儲與操作與最外層網格一致，以保證嚴格的對應關係。

SGOG網格高程匹配

要建立基於ESSG的大區域地形模型，必須將傳統的基於投影的DEM網格高程「映射」到離散網格上。這裡藉助IDL語言實現，以經緯度坐標為橋梁，將WGS-84橢球近似為球體後採用最近鄰像元法匹配高程。其高程匹配方法如圖8所示。

具體步驟為：

1、在IDL平臺下讀入存儲的SGOG最外層網格點X、Y坐標，並輸入地球半徑，換算為經緯度；

2、遍歷DEM影像數據，根據影像的最大、最小行列號，換算出影像的經緯範圍，判斷網格點是否落在當前影像中；

3、網格點落在當前影像則讀入影像高程值，並把落在影像內部的網格點經緯度換算為所處當前影像的行列號，根據行列號讀出像元值賦給網格點高程，如果網格點不在當前影像範圍內則判斷下一幅影像。

在高程匹配過程中，各種坐標的變換是關鍵環節。DEM影像中像元值按行列存儲，而網格點位置是以右手坐標系下的三維坐標形式表達，所以需要把網格點的三維坐標和影像像元的行列號都換算為經緯度，才能進行網格點定位，匹配像元值。在匹配完所有DEM影像數據後，由於在OSG平臺可視化時使用的是右手坐標系下的三維坐標，則需要再將所有網格點的經緯度和匹配好的高程值轉換成相應的三維坐標。

真三維場景可視化

DEM可視化基於VS2010和OSG平臺實現。在進行最外層網格點高程的基礎上，按照所有可視範圍的SGOG網格點坐標，利用網格可視化技術繪製網格，實現大區域真三維地理可視化。為突出DEM高低起伏的效果，將高程擴大30倍。由於機器性能的限制，為達到較好的顯示效果，靈活設定模型的網格層次。圖9是利用第9層的近7萬個面網格疊加第4層體網格建立的中國大陸地區的三維地形框架。

大區域真三維地理場景自適應可視化

根據中國大陸地區的地形起伏度和建模需求，設定相應閾值，靈活採用多層網格體系，例如採用第6、7、8、9層網格進行自適應可視化。具體的方法是：選取第6層次作為基層網格，然後依次遞增網格層級進行自適應剖分建模。對於第6層三角形網格而言，取每個三角形三個頂點的最大高程差代表該三角形的地勢起伏，取該層所有三角形地勢起伏的平均值作為該層的地勢起伏度，設置為閾值。大於等於此值，則認為該三角形內地勢起伏較大，則繼續對其進行第7層剖分處理；小於此值，則認為地勢平坦，不再對三角形進行剖分處理。對於第7、8層按照上述方法處理，以第9層為最終分層網格。但是在實驗過程中，發現採用上述的閾值，得到的自適應效果並不好，網格過於稀疏，因此我們將該閾值修改為原來的三分之一，即為取該層所有三角形最大高程差的平均值的三分之一作為閾值。效果如圖10、圖11所示。

模型高程表達精度分析

源數據、圖9和圖10模型的地面高程表達的特徵值統計如表1所示：

表1各種DEM模型高程表達的特徵值

從表中可以看出，第9層網格和自適應網格的點數都遠遠小於源數據，自適應網格是第9層網格的1/4。兩者的高程最小值都遠大於源數據的最小值，自適應網格比第9層高28m。二者高程最大值是一樣的，比源數據低2262米。從平均值來看，第9層網格比源數據低830m，而自適應網格要比源數據高530m，二者相差1383m。從最值來衡量，二者的精度相當。從均值來看，二者相差較大。90mDEM網格與SGOG第9層網格屬於不同的剖分類型，二者網格點在進行「匹配」時，具有一定的偶然性，它們之間的差異主要由偶然因素引起；而從自適應網格的剖分方法可知，高差越大，剖分越細，因而它與前述兩種方法相比，總是「趨高」。

場景渲染

根據網格點的不同高程值為每個網格點賦上相應的顏色值，利用OSG渲染引擎自帶的函數，隨著瓦塊編碼的遍歷而繪製瓦塊的各個表面並進行顏色渲染，可得到中國大陸地區的DEM場景暈渲圖，如圖12和圖13所示。青藏高原側面的瀑布狀渲染，並不代表積雪的覆蓋範圍，而是由於渲染技術細節不完善造成，有待進一步優化。

選擇四川盆地局部地區，以地質分層示意的SGOG第8、9、10層構建三維地殼框架，表面粘貼遙感影像紋理，得到地理場景效果如圖14和圖15所示。