顧及計算強度的地質矢量大數據計算域 均衡分解方法研究

朱月琴，郭明強，黃 穎，吳 亮,5，謝 忠,5

(1.自然資源部地質信息技術重點實驗室，北京 100037； 2.中國地質調查局發展研究中心，北京 100037； 3.中國地質大學(武漢)信息工程學院，湖北 武漢 430074； 4.武漢中地數碼科技有限公司，湖北 武漢 430074； 5.國家地理信息系統工程技術研究中心，湖北 武漢 430074)

0 引 言

隨著科學研究方式第四范式的誕生，即數據密集型的知識發現，地質學領域具有多源、異構等復合型數據已經被列入大數據的范疇，被稱之為“地質大數據”[1-2]。大數據量地理信息數據在線服務應用提供的矢量服務，其響應速度和數據量、服務器性能，網絡環境相關性大，目前矢量大數據服務應用存在的瓶頸有以下四點：數據發布準備工作費時費力[3-4]、時效性難以保證[5-6]、無法在線編輯處理、難以支持大數據量應用[7-9]，迫切需要研究適應矢量大數據實時渲染的并行處理方法。

1 地質矢量大數據實時渲染系統總體框架

面向地質大數據實時渲染的計算域均衡分解系統層次架構的設計包括：數據準備層、前端應用層、業務邏輯層、GIS服務層。數據準備層為獨立的桌面軟件，與后面層級沒有直接的關系，但后面層級在運行之前必須首先通過數據準備層生成必要的準備數據。前端應用層通過ajax向業務邏輯層發送http請求，業務邏輯層通過多線程并行的向多個MapGIS IGServer服務節點發送REST服務請求獲得矢量地圖圖像，如圖1所示。

圖1 系統架構圖

數據準備層：該層為一個獨立的桌面軟件，即桌面自動采集工具，主要的功能是對矢量大數據進行樣本自動采集，并通過對采集的樣本數據進行統計分析，得到矢量大數據渲染時間的計算公式。同時實現了對矢量大數據自動網格化，為后面矢量大數據優化做數據準備工作。

前端應用層：該層為面向地質大數據實時渲染的計算域均衡分解系統的前端部分，主要實現了地圖顯示、地圖基本操作，包括放大，縮小，復位等功能，以及動態的監視并顯示每一個服務器節點渲染矢量數據所耗費的時間。

業務邏輯層：該層為本系統的核心部分，主要實現了最佳網格級別查找、可視域范圍內的網格單元檢索、可視域范圍內矢量要素渲染時間統計、可視域范圍的均衡劃分、通過多線程對矢量數據并行渲染、無縫拼接渲染結果等一系列功能。

GIS服務層：即MapGIS IGServer，該部分主要負責地質矢量大數據服務的管理與服務發布，提供了豐富的GIS功能服務以及完善的二次開發接口。

2 地質矢量大數據實時渲染關鍵技術

2.1 矢量地圖服務時間樣本自動采集

矢量地圖服務時間樣本的自動采集主要實現了在矢量數據的空間范圍內隨機生成大小、位置各異的矩形范圍，利用每次生成的矩形范圍對矢量數據進行空間查詢，記錄當前矩形范圍內要素的個數。并且利用以當前隨機矩形范圍作為矢量地圖的渲染范圍，從服務器獲取圖像，同時記錄獲取圖像所耗費的時間。將每一個隨機矩形范圍對應的要素總數及矢量地圖的渲染時間記錄到樣本文件并保存至指定的文件夾下。實現流程圖如2所示。

圖2 矢量地圖服務時間樣本自動采集實現流程圖

2.2 矢量地圖服務時間公式自動統計

矢量地圖服務時間公式自動統計的功能是在矢量地圖服務時間樣本自動采集完成的基礎上進行的，其目的是對采集的樣本數據通過一元線性回歸的方法，計算出在矢量要素渲染的過程中要素數量與渲染時間之間的線性關系CI=a×count+b中的系數a與截距b，以及R2。其主要的實現思路為：首先從采集的樣本數據中獲取多組觀察值(count1，CI1)，(count2，CI2)，…，然后調用已封裝好的最小二乘法計算接口，求出系數a、截距b以及R2。

實現效果圖如3所示。

圖3 矢量地圖服務時間公式自動統計實現效果圖

2.3 矢量圖層時間網格自動生成

矢量圖層時間網格自動生成主要實現了將矢量圖層使用四叉樹的方式進行網格劃分，劃分的實現思路為：網格組織方式按照四叉樹(4n)組織，第0級一個網格，第1級4個網格，依次類推。算法：首先獲取圖層的外包矩形范圍[Xmin，Ymin，Xmax，Ymax]，然后將其擴展為一個正方形范圍[XRmin，YRmin，XRmax，YRmax]，擴展分為兩種情況：

第一種，外包矩形寬度(Xmax-Xmin)大于高度(Ymax-Ymin)時，擴展計算公式見式(1)～(4)。

XRmin=Xmin

(1)

(2)

XRmax=Xmax

(3)

(4)

第二種情況，外包矩形寬度(Xmax-Xmin)小于高度(Ymax-Ymin)時，擴展計算公式見式(5)～(8)。

(5)

YRmin=Ymin

(6)

(7)

YRmax=Ymax

(8)

計算第0～n級網格的分辨率，計算公式見式(9)～(11)。

(9)

(10)

…

(11)

計算指定級別l下第r行第c列網格單元的矩形范圍[XCmin，YCmin，XCmax，YCmax],規定左下角(XRmin，YRmin)為計算原點，具體計算公式見式(12)～(15)。

XCmin=XRmin+c.Rl

(12)

YCmin=YRmin+r.Rl

(13)

XCmax=XCmin+Rl

(14)

YCmax=YCmin+Rl

(15)

在MapGISLocal數據源指定數據庫中創建區圖層，將計算的每個網格單元的矩形范圍作為幾何參數，該網格單元的要素個數、渲染時間、行號、列號等作為屬性參數添加到上述新建的區圖層中。

2.4 矢量地圖服務均勻劃分

矢量地圖服務均勻劃分首先需要在均衡分解模塊的全局配置文件中配置結點數量c和一個均衡參數p。如現有部署方式為：1+2模式，1臺機器為示例網站部署機器，2臺MapGIS IGServer服務器，則c=2。可視域范圍采用網格劃分方法，則可視域到底在哪個分辨率的網格下劃分最均衡至關重要。為計算最佳網格分辨率，設置了均衡參數p。設現有可視域范圍為[XVmin，YVmin，XVmax，YVmax],最佳網格分辨率Rp計算公式見式(16)。

(16)

計算得到最佳分辨率Rp之后，通過遍歷所有級別的網格的分辨率，得到分辨率最接近Rp的網格的級別數lp。然后利用當前可視域范圍對級別數為lp的網格區圖層進行空間查詢，查詢該范圍內的所有計算時間網格單元矩形要素。若可視域范圍不能完全覆蓋正方形范圍，則邊界求占比，修改可視域范圍覆蓋不完整的網格單元的渲染時間并且修改網格單元的矩形范圍。經過研究分析，將視域范圍不能完全覆蓋的網格單元按照其所在不同的位置分別計算其邊占比，參與邊占比計算的網格單元的分布情況如圖4所示，具體的算法如下所述。

圖4 參與邊占比計算的網格單元分布圖

為了便于閱讀下列公式，此處給出公式相關參數含義：網格單元矩形范圍為[XCmin，YCmin，XCmax，YCmax],可視域矩形范圍為[XVmin，YVmin，XVmax，YVmax],修改后單元格范圍為[XNmin，YNmin，XNmax，YNmax],邊占比為Ratio，修改后網格單元的渲染時間為CInew。

左上角單元格計算公式為式(17)。

(17)

修改XCmin為XVmin，YCmax為YVmax。

左邊單元格計算公式式(18)。

(18)

修改XCmin為XVmin。

左下角單元格計算公式為式(19)。

(19)

修改XCmin為XVmin，YCmin為YVmin。

下邊單元格計算公式為式(20)。

(20)

修改YCmin為YVmin。

右下角單元格計算公式為式(21)。

(21)

修改XCmax為XVmax，YCmin為YVmin。

右邊單元格計算公式為式(22)。

(22)

修改XCmax為XVmax。

右上角單元格計算公式為式(23)。

(23)

修改XCmax為XVmax，YCmax為YVmax。

上邊單元格計算公式為式(24)。

(24)

修改YCmax為YVmax。

修改渲染時間的公式為式(25)。

CInew=CI×Ratio

(25)

網格數據檢索完成之后，遍歷所有的網格單元，計算當前可視域范圍內網格單元內矢量要素的渲染時間總和。根據服務器節點數量求出平均時間。利用可視域范圍內的網格單元進行可視域劃分。劃分思路可簡述為：將可視域范圍內的所有的單元格按照其行號進行分行，遍歷單元格行，累加單元格渲染時間同時累加每一個單元格的空間范圍，直至渲染時間總和大于等于平均渲染時間為止，此處分為兩種情況：第一種，渲染時間總和等于平均渲染時間，則直接將空間范圍的累加結果分配給第一個服務器；第二種，渲染時間總和大于平均渲染時間，則計算超出平均時間部分與平均時間的占比，根據該占比調整空間范圍，然后再將該范圍分配給第一個服務器節點，多出的空間范圍則累計至下一個服務器節點的空間范圍，依此類推。實現流程圖如圖5所示。

3 結 語

本文在地質大數據爆發式增長的背景下，結合目前在線大數據地圖服務在實際應用中的瓶頸。面向地質大數據實時渲染的計算域均衡分解系統采用了一系列的優化方案，包括矢量地圖服務時間樣本自動采集、矢量地圖服務時間工時自動統計、矢量地圖網格化、基于網格的可視域劃分、多線程并行渲染，渲染結果合并等，在矢量大數據在線服務的響應速率提升方面獲得了比較明顯的效果。其中基于網格的可視域劃分是該系統的核心部分，結合矢量大數據的特征，采用網格單元將矢量要素的渲染時間拆分，再根據服務器節點的數量按照平均渲染時間進行網格單元合并的劃分方案，保證了每個服務器節點的渲染時間基本相等，負載均衡效率提升明顯。

圖5 矢量地圖服務均勻劃分實現流程圖