一種基于地球剖分網(wǎng)格的區(qū)域面積計算方法

齊澄宇,程承旗,濮國梁,陳 波

(1. 北京大學(xué) 遙感與地理信息系統(tǒng)研究所,北京 100871;2. 北京大學(xué) 工學(xué)院,北京 100871)

0 引 言

地球剖分網(wǎng)格是一種科學(xué)簡明的空間參考系統(tǒng),是地球信息數(shù)據(jù)的離散化、準確化、系統(tǒng)化的表現(xiàn)[1].在地球剖分的思路上,產(chǎn)生了多種全球離散格網(wǎng)系統(tǒng)(Discrete Global Grids System, DGGS)[2].Goodchild、郭華東等認為,DGGS采用預(yù)先計算的瓦片結(jié)構(gòu)優(yōu)化本地計算,提高使用中的細節(jié)管理水平,這是"虛擬地球最吸引科學(xué)工作者的特點之一"[3].特別是隨著對地觀測能力的增強,空間數(shù)據(jù)體量迅猛增長.利用網(wǎng)格進行數(shù)據(jù)組織、用帶有位置信息的網(wǎng)格編碼替代經(jīng)緯度來記錄數(shù)據(jù)的位置屬性,成為多源異構(gòu)空間大數(shù)據(jù)解決方案.其思路在于,在適用于大數(shù)據(jù)的新型列數(shù)據(jù)庫中,使用網(wǎng)格編碼作為數(shù)據(jù)索引主鍵,實現(xiàn)數(shù)據(jù)與位置的自動關(guān)聯(lián).相關(guān)研究包括真三維數(shù)據(jù)表達[4]、警用大數(shù)據(jù)[5]、遙感大數(shù)據(jù)[6]、空間對象快速可視化[7]等.

地表區(qū)域面積計算一直是測繪科學(xué)的一個重要問題,在國土面積測量、土地資源管理、民政、遙感、減災(zāi)等方面都有實際應(yīng)用.在計算機輔助軟件產(chǎn)生之前,基于地圖數(shù)據(jù)的面積測量常用的方法是網(wǎng)格法、求積儀法.地理信息系統(tǒng)產(chǎn)生之后,在常用的矢量和柵格兩類數(shù)據(jù)模型下,借助計算機,可以方便地進行面積計算.如矢量模型以點為最基本的實體組成.計算矢量多邊形的面積,一種可行的方法是根據(jù)點的坐標借助格林公式[8]精確計算.

網(wǎng)格型地理信息系統(tǒng)有望憑借其在處理大體量多源異構(gòu)數(shù)據(jù)中的優(yōu)勢,成為大數(shù)據(jù)地理信息系統(tǒng)的解決方案,基于網(wǎng)格的空間量算變得尤為重要.在剖分框架下,實體和數(shù)據(jù)以剖分網(wǎng)格組織,實體的面積就是其對應(yīng)的網(wǎng)格集的面積之和.問題在于,球面空間與歐氏空間不同胚,導(dǎo)致同一層級剖分網(wǎng)格的面積不盡相同.以等經(jīng)緯度剖分網(wǎng)格為例,由于不同緯度下相同經(jīng)度差對應(yīng)的地表長度不同,因此同一層級的網(wǎng)格由赤道到兩極,面積變小.計算剖分網(wǎng)格集對應(yīng)的地表范圍面積,常常需要先將網(wǎng)格對應(yīng)回經(jīng)緯度范圍,再借助積分方法求解,地表區(qū)域面積計算效率低.

針對以上問題,本文提出了一種新的地表面積計算方法,以GeoSOT地球空間剖分參考框架為基礎(chǔ),根據(jù)GeoSOT剖分網(wǎng)格幾何特性,結(jié)合地表區(qū)域面積計算的需求,構(gòu)建全球網(wǎng)格弧長數(shù)據(jù)庫,以查表求和方法取代經(jīng)緯度積分方法,實現(xiàn)在該框架下地表區(qū)域面積快速、精確計算,有望應(yīng)用于測繪、遙感、減災(zāi)、國土資源管理等諸多領(lǐng)域.

1 GeoSOT地球剖分方案

基于空間信息剖分組織理論,北京大學(xué)程承旗團隊[9]提出"2n一維整型數(shù)組的全球經(jīng)緯度剖分參考網(wǎng)格"(Geographic Coordinate Subdivision Grid with One Dimension Integral Coding on 2n-Tree,簡稱GeoSOT),以探尋建立一種空間信息或數(shù)據(jù)組織的專用網(wǎng)格,并能夠與現(xiàn)有的主要空間數(shù)據(jù)組織網(wǎng)格有很好的區(qū)位關(guān)聯(lián)和尺度聚合關(guān)系;為空間信息的統(tǒng)一區(qū)位組織和區(qū)域關(guān)聯(lián)調(diào)度,構(gòu)建一種更適合空間信息組織的區(qū)位標識體系.GeoSOT全球剖分方案提出之后,很多學(xué)者針對多源數(shù)據(jù)在GeoSOT剖分參考框架下的表達、組織、管理、分發(fā)等技術(shù)與應(yīng)用層面展開了深入研究,如氣象水文[10]、減災(zāi)[11-12]、地理國情監(jiān)測[13]、遙感[14-15]、導(dǎo)航[16]、無人機[17]等.并且由于GeoSOT剖分參考框架的科學(xué)性與實用性,相關(guān)部分成果已在國內(nèi)高分專項、北斗專項等重大項目應(yīng)用,并正在形成國防及公安、民政、新型智慧城市等重要行業(yè)部門的國家應(yīng)用標準及草案.

GeoSOT網(wǎng)格剖分以國家大地坐標系CGCS2000為空間基準;沿用大地測量理論與技術(shù)體系,能夠最大限度地繼承現(xiàn)有歷史空間數(shù)據(jù);包含現(xiàn)有不同網(wǎng)格體系公共的4°、2°、1°、2′、1′、2″、1″、0.5″八個基本網(wǎng)格;能夠聚合生成現(xiàn)有主要基于經(jīng)緯度的標準網(wǎng)格.其剖分方案如下:

通過三次地球擴展(將地球表面真實空間擴展為512°X512°、將1°擴展為64′、將1′擴展為64″)實現(xiàn)整度、整分的四叉樹劃分,從而形成一個上至全球(0級)、下至厘米級邊長(32級)的多尺度網(wǎng)格體系[18],如圖1所示.

圖1 GeoSOT地球剖分參考框架[9]Fig.1 GeoSOT reference frame

在GeoSOT格網(wǎng)劃分的基礎(chǔ)上,按照"Z"序為每個格網(wǎng)賦予層次性編碼,編碼形式主要有四進制一維數(shù)組(如圖2所示)、二進制一維數(shù)組(四進制編碼的二進制形式)和二進制二維數(shù)組(二進制一維數(shù)組的交叉分組)3種,均由度級、分級、秒及秒小數(shù)三部分組成.例如,四進制一維形式的編碼12,其二進制一維形式的編碼為0110,二進制二維形式的編碼為(01,10).

圖2 GeoSOT格網(wǎng)編碼示意圖[9]Fig.2 GeoSOT coding rule

2 基于GeoSOT的地表區(qū)域面積計算方法

2.1 計算原理

本文基于構(gòu)建全球網(wǎng)格弧長數(shù)據(jù)庫的思路提出一種新的計算方法,構(gòu)建全球網(wǎng)格弧長數(shù)據(jù)庫,采用預(yù)先計算、查表求和的方式,取代網(wǎng)格轉(zhuǎn)換為經(jīng)緯度范圍再進行積分的計算方式,以期獲得更高的計算效率.方案如下:

首先,采用數(shù)學(xué)方法計算基層級每一網(wǎng)格的邊長.基層級需在19級4″及以下層級選定,以確保不包含不規(guī)則層級.隨后,根據(jù)網(wǎng)格擴展方法,設(shè)計數(shù)據(jù)查找規(guī)則.對于輸入的網(wǎng)格編碼,確定其對應(yīng)的數(shù)據(jù)塊,進行求和操作,計算網(wǎng)格的面積,流程如圖3所示.

圖3 流程圖Fig.3 Flow chart

2.2 經(jīng)緯線弧長計算

經(jīng)緯線弧長的計算方法可以根據(jù)幾何方法求解[17].在地球橢球面上,緯線是圓.在任意一條緯線φ上求兩點之間經(jīng)差為μ的緯線弧長,可以利用下式:

地球橢球面上,經(jīng)線圈是橢圓,曲率半徑隨著緯度變化.求兩條緯線φ1φ2之間經(jīng)線長度,需要進行積分,公式如下:

兩式中,ae表示地球橢球赤道半徑,e1表示第一偏心率.按照公式可以精確計算經(jīng)緯度弧長.

GeoSOT采用CGCS2000大地坐標系,因此,本文地球橢球參數(shù)采用此大地基準給出的相關(guān)參數(shù),與GeoSOT保持一致.CGCS2000大地基準中,ae=6 378137 m,扁率f=1/298.257 222101 ,利用扁率計算第一偏心率,公式為:

在計算中,取e12=0.006694380022 901.

GeoSOT進行經(jīng)緯度劃分時,在南北緯度88°以內(nèi)采用四等分,在南北極處采用不同的劃分策略.考慮到對稱性,本研究只進行北緯0°~88°的計算,南半球的數(shù)據(jù)與北半球相同.

2.3 弧長查詢方法設(shè)計

GeoSOT網(wǎng)格上至全球,下至1/2 048″,共分為32個層級,每一層級的網(wǎng)格經(jīng)緯度跨度類似但不盡相同,主要原因在于GeoSOT的第二次和第三次擴展.在第9級向第10級劃分時,首先將1°網(wǎng)格擴展到64′,隨后再進行劃分.1°網(wǎng)格向下劃分到第10級,得到的4個網(wǎng)格跨度分別是32′X32′,32′X28′,28′X32′以及28′X28′,如圖4a所示.對這4個網(wǎng)格進一步劃分,將得到16個11級網(wǎng)格,其中16′X16′的9個,16′X12′的3個,12′X16′的3個,12′X12′的1個.進一步劃分將得到64個網(wǎng)格,其中大部分是規(guī)則的8′網(wǎng)格,但仍有不規(guī)則網(wǎng)格,直到再下一次劃分,才得到統(tǒng)一的共計225個4′網(wǎng)格.

圖4 不規(guī)則網(wǎng)格的產(chǎn)生Fig.4 Irregular grids generation

在1′網(wǎng)格向下劃分時,會出現(xiàn)相似的過程.不規(guī)則網(wǎng)格的出現(xiàn)給網(wǎng)格面積計算帶來了困難.如28′X28′的網(wǎng)格與32′X32′網(wǎng)格層級一致,但面積比其小約1/4.

由球面對稱性,可僅計算北緯0°~88°、東經(jīng)0°~1°的單位緯線、經(jīng)線長度,其余由旋轉(zhuǎn)對稱性和軸對稱性可以求得.為方便計算,采用第22層級1/4″網(wǎng)格作為最小跨度,原因如下:1/4″網(wǎng)格精確度到米級,可以滿足較高精度的使用需求;并且,1/4″網(wǎng)格精度高于由于擴展而產(chǎn)生的不規(guī)則網(wǎng)格的最小尺度,不存在殘缺的網(wǎng)格,在進行更為細小的劃分時,只需將長度進行均分即可.使用Matlab,根據(jù)式(2)以22層級網(wǎng)格1/4″為跨度進行積分,將結(jié)果保存數(shù)組中.數(shù)組的每一項表示北緯0°~88°,每條1/4″經(jīng)線的長度,共有4X60X60X88=1 267200項.隨后對其進行如下變換.

對第一個數(shù)組各項進行求和,可以得到0°緯線和88°緯線之間經(jīng)線段的長度.將其變換為(4X60)X(60X88)的矩陣,則此時矩陣的每一列之和代表1矩經(jīng)線段的長度.在這一矩陣下鄰接一個(4X4)X(60X88)的0矩陣,更新為(4X64)X(60X88)的新矩陣,每列之和仍代表1°經(jīng)線段的長度,但行數(shù)由4X60變到了4X64.對其進行一次劃分,得到(4X32)的兩個矩陣,第一個矩陣每一項都為正數(shù),而第二個矩陣的最后4X4行都由0組成.這樣的一次矩陣變換實現(xiàn)了1°向64′擴展操作在經(jīng)線方向上的投影.

將矩陣鄰接之后得到的(4X64)X(60X88)矩陣變換為(4X64X60)X88的新矩陣,則矩陣的每一列之和代表1°經(jīng)線的長度,88列代表共88條1°經(jīng)線.在其下方鄰接(4X64X4)X88的0矩陣,則同上一次矩陣鄰接操作類似,實現(xiàn)了1′向64″擴展操作在經(jīng)線方向上的投影.最后,把矩陣整理為4X64X64X88項的一個向量,記為Lon,則向量Lon代表了從0°到88°的GeoSOT網(wǎng)格在經(jīng)線方向上的投影.變換過程示意如圖5所示.產(chǎn)生的結(jié)果向量_Lon、_Lat1以及_Lat2.隨后,求取_Lon、_Lat1兩個向量的內(nèi)積,并數(shù)乘_Lat2向量的各項之和,結(jié)果即為該網(wǎng)格編碼對應(yīng)的地表面積.計算過程可以用下式表示:

對網(wǎng)格編碼集中各面積進行求和,就可以得到所表示的地表區(qū)域面積.其中,由GeoSOT四進制編碼求取對應(yīng)劃分的方法如下:設(shè)網(wǎng)格A的四進制一維編碼Codei=a1a2…an,這里假設(shè)A在東北半球,因此,a1=1對于其他情況,由對稱性可以方便求取.令i=2,3,…,n,對于ai%4=0或1,取Lon、Lat1的右區(qū)間,否則取左區(qū)間;對于ai%2=0,取Lat2左區(qū)間,否則取右區(qū)間,直到窮盡四進制一維編碼的每一位,最終得到_Lon、_Lat1以及_Lat2.

圖5 變換過程示意Fig.5 Matrix transformation process

3 實驗結(jié)果與分析

經(jīng)過矩陣變換,原來共計4X60X60X88項的數(shù)組,變換成為4X64X64X88的向量.任意一個1至22層級網(wǎng)格都能在向量Lon上找到一個區(qū)間來對應(yīng),區(qū)間包含的所有項之和構(gòu)成該網(wǎng)格的經(jīng)向長度.

緯向長度與經(jīng)向長度略有不同.沿經(jīng)線方向,不同緯度的同一跨度緯線段長度不相同,并且會受到擴展的影響;沿緯線方向,不同緯度的同一跨度緯線段長度相同,但也會受到擴展的影響.首先,沿經(jīng)線方向計算不同緯度的同一跨度緯線段長.根據(jù)式(1),以1/4″為跨度進行計算,得到一個4X60X60X88項的數(shù)組.接下來的處理方法類似于經(jīng)線長度,最終得到4X64X64X88項的一個向量,記為Lat1,代表0~88表緯線每隔1/4隔跨度的1/4隔緯線段長度.緯線方向的處理較為簡單,不需要進行長度的實質(zhì)計算.對于1線的跨度范圍,構(gòu)造一個4X60X60長度、各項為1的向量.隨后,在向量的每隔4X60項插入4X4個0元素,并且在最后插入4X64X4個0元素,得到一個4X64X64的新向量,記為Lat2,其所有項由1和0組成.

在以上3個向量的基礎(chǔ)上,以GeoSOT網(wǎng)格編碼集表示的地表區(qū)域面積計算方法為:

對于網(wǎng)格編碼集C中每一網(wǎng)格編碼Codei,首先,根據(jù)網(wǎng)格編碼,在向量Lon、Lat1和Lat2上找到對應(yīng)劃分所

根據(jù)文章第2部分所述的方法建立全球網(wǎng)格弧長數(shù)據(jù)庫,數(shù)據(jù)文件總大小50 Mb左右.在此基礎(chǔ)上,進行可行性實驗及效率實驗.實驗環(huán)境:4 GB內(nèi)存,Intel(R)Core(TM) i5-3230M 2.60GHz CPU,Win10 64-bit,軟件環(huán)境為Matlab 2015b,實驗過程及結(jié)果如下.

3.1 可行性實驗與分析

可行性實驗的主要思路是利用文中方法和傳統(tǒng)借助經(jīng)緯度的方法分別計算GeoSOT編碼集代表的地表區(qū)域面積,對結(jié)果進行比較,通過相對誤差大小來檢驗算法的可行性.

其中,經(jīng)緯度法求地球表面規(guī)則區(qū)域面積方法的計算公式為:

式中,φ1、φ2代表始、終緯度,μ代表經(jīng)度跨度.

采用隨機數(shù)生成器,模擬生成東北半球不同層級的GeoSOT網(wǎng)格四進制一維編碼,作為算法輸入.實驗進行500次,生成隨機長度的四進制一維網(wǎng)格編碼,使用兩種方法分別進行計算并對結(jié)果進行比較.其中,無效網(wǎng)格編碼78個.剩余422次結(jié)果的相對誤差大多分布在±10-5之內(nèi),少量突破了10-5,但也在10-4以內(nèi).

實驗結(jié)果表明,本文提出的計算方法在不同層級的計算中均能得到較準確的結(jié)果,計算相對誤差在0.000 1以內(nèi).本文提出的基于GeoSOT的面積計算方法在不同層級的面積計算中具有可行性.

3.2 效率實驗與分析

對于隨機生成的東北半球的GeoSOT網(wǎng)格四進制一維編碼,使用Matlab分別利用本文提出的方法以及積分方法進行多次計算,記錄并比較程序總運行時間.實際測試時,每次進行2 000次面積計算,記錄總用時.50次比較結(jié)果如圖6所示,其中橫軸表示實驗次數(shù),左縱軸表示運行時間,右縱軸表示經(jīng)緯度積分方法與本文方法運行時間之比.

圖6 時間記錄Fig.6 Time record

實驗結(jié)果顯示,本文提出的方法在運行時間上優(yōu)于傳統(tǒng)方法,效率約為傳統(tǒng)方法的7.5倍左右.傳統(tǒng)方法使用復(fù)雜的積分進行面積計算.隨著待求區(qū)域?qū)蛹壧岣?、面積變大,積分時間顯著提高.本文提出的方法則采用事先計算的方案,從存儲器中讀出結(jié)果進行求和,使得計算效率大大提高.

4 結(jié)束語

本文提出了一種基于GeoSOT的區(qū)域面積計算方法,該算法以GeoSOT地球空間剖分參考網(wǎng)格為基礎(chǔ),根據(jù)GeoSOT剖分方案,設(shè)計編碼查詢規(guī)則,構(gòu)建全球網(wǎng)格弧長數(shù)據(jù)庫,根據(jù)編碼查表計算網(wǎng)格面積.該方法基于全球剖分組織框架設(shè)計,具有全球統(tǒng)一性.通過應(yīng)用實驗,驗證了該算法計算不同層級網(wǎng)格編碼的可行性,并且相比傳統(tǒng)方法,時間效率有大幅提高.該方法為地表空間面積計算提供了一種新思路.全球離散格網(wǎng)系統(tǒng)有望在大數(shù)據(jù)時代成為空間數(shù)據(jù)組織體系的標準框架,離散格網(wǎng)系統(tǒng)下的立體空間量算,是下一步研究中應(yīng)重點考慮的問題.