三調(diào)橢球面積計(jì)算工具設(shè)計(jì)與實(shí)現(xiàn)

閭海洋，湯 玲，王 剛

(江蘇省地質(zhì)測(cè)繪院，江蘇 南京 211102)

0 引 言

國務(wù)院于2017年10月16日下發(fā)《國務(wù)院關(guān)于開展第三次全國土地調(diào)查的通知》(國發(fā)〔2017〕48號(hào))，決定自2017年起開展第三次全國土地調(diào)查(以下簡稱“三調(diào)”)。為統(tǒng)一三調(diào)技術(shù)標(biāo)準(zhǔn)，方便各個(gè)縣區(qū)的數(shù)據(jù)庫能夠融合到一起，自然資源部于2019年1月28日發(fā)布了《第三次全國國土調(diào)查技術(shù)規(guī)程》(TD/T 1055-2019，以下簡稱《三調(diào)技術(shù)規(guī)程》)。三調(diào)以縣(區(qū))為單位進(jìn)行調(diào)查，每個(gè)縣都有一個(gè)控制面積，每個(gè)縣內(nèi)所有圖斑的面積和等于縣控制面積，所有縣的面積和等于國家的控制面積。通過這種方式進(jìn)行控制，能夠保證全國的面積是一個(gè)固定值，不會(huì)因?yàn)橛?jì)算方式的不同而產(chǎn)生不同的總面積。三調(diào)圖斑沒有采用高斯投影面上圖形的幾何面積，而是使用橢球面積作為圖斑的面積[1]。計(jì)算圖斑橢球面積這一功能相對(duì)復(fù)雜，在主流的GIS建庫軟件中都沒有現(xiàn)成的工具，因此需要使用GIS軟件提供的開發(fā)接口進(jìn)行二次開發(fā)。

1 面積計(jì)算公式解析

《三調(diào)技術(shù)規(guī)程》中對(duì)圖斑的橢球面積計(jì)算方式做出了明確的規(guī)定：對(duì)任意一個(gè)封閉區(qū)域(圖斑)進(jìn)行分割，將封閉區(qū)域分割成多個(gè)梯形圖塊，先對(duì)每個(gè)梯形圖塊進(jìn)行計(jì)算，然后再對(duì)所有梯形的面積求和從而得到封閉區(qū)域的面積。對(duì)于橢球面上的任意一個(gè)梯形圖塊的應(yīng)采用公式(1)進(jìn)行計(jì)算：

(1)

其中，a表示橢球的長半軸，單位為m;b表示橢球的短半軸，單位為m;ΔL表示梯形圖塊的經(jīng)度差值，單位為弧度；BΔ表示梯形圖塊的緯度差值，單位為弧度；Bm表示梯形圖塊的南北緯度的平均值，單位為弧度。

《三調(diào)技術(shù)規(guī)程》規(guī)定平面坐標(biāo)使用2000國家大地坐標(biāo)系，因此參數(shù)α、b、e2、A、B、C、D、E都為常數(shù)，在計(jì)算圖斑面積之前可以固化在代碼中，具體數(shù)值如表1所示。

表1 基本常數(shù)表

圖斑是土地實(shí)際的利用狀況在數(shù)據(jù)庫的表現(xiàn)，都是不規(guī)則的多邊形，因此需要對(duì)圖斑進(jìn)行分割處理，將不規(guī)則的圖斑分割成若干個(gè)規(guī)則的梯形圖塊進(jìn)行計(jì)算，分割方式如圖1所示。

圖1 任意圖斑梯形分割示意圖

因此圖斑ABCDE的總面積公式可以表示為公式(2)：

SABCDE=SABB1A1+SBCC1B1+SCDD1C11+SDEE1D1+SEAA1E1

(2)

式中的每個(gè)面積子項(xiàng)都可以通過梯形圖塊面積計(jì)算公式得到。計(jì)算梯形圖塊面積時(shí)，使用的坐標(biāo)是地理坐標(biāo)，而實(shí)際生產(chǎn)過程中使用的坐標(biāo)是高斯平面坐標(biāo)，在計(jì)算時(shí)還需要通過高斯投影反算公式將高斯平面坐標(biāo)反算為地理坐標(biāo)。

2 長邊問題處理

在計(jì)算圖斑的橢球面積過程中，為提高面積計(jì)算精度，如果圖斑節(jié)點(diǎn)上任意相鄰節(jié)點(diǎn)之間的線段長度大于70 m時(shí)，需要進(jìn)行插值計(jì)算。內(nèi)插點(diǎn)僅作為橢球面積計(jì)算使用，不需要對(duì)線段進(jìn)行分割。進(jìn)行內(nèi)插時(shí)，內(nèi)插點(diǎn)應(yīng)分布均勻，插值計(jì)算流程如圖2所示。

圖2 插值計(jì)算流程

3 面積平差方法

通過面積計(jì)算公式計(jì)算出每一個(gè)圖斑的橢球面積，但由于計(jì)算過程中會(huì)產(chǎn)生一些誤差，最終導(dǎo)致面積匯總時(shí)總面積會(huì)出現(xiàn)不確定性。為解決這一問題，三調(diào)采用了控制面積制度。國土調(diào)查數(shù)據(jù)庫中的縣級(jí)及縣級(jí)以上行政區(qū)域界線采用全國陸地行政區(qū)域勘界成果確定的界線，鄉(xiāng)鎮(zhèn)級(jí)行政區(qū)域界線采用各縣(市、區(qū))最新確定的界線。為保證全國各個(gè)縣(市、區(qū))的國土調(diào)查數(shù)據(jù)庫拼到一起后的總面積的確定性，在計(jì)算面積時(shí)實(shí)行了控制面積制度。在開展三調(diào)建庫時(shí)，國家統(tǒng)一下發(fā)了各個(gè)縣(市、區(qū))的控制面積，在一個(gè)縣(市、區(qū))范圍內(nèi)，所有圖斑的總面積應(yīng)等于下發(fā)控制面積。面積計(jì)算以平方米為單位，保留2位有效小數(shù)。

在三調(diào)中，圖斑號(hào)以村為單位進(jìn)行編號(hào)，控制面積實(shí)行三級(jí)控制，即縣控制各個(gè)鄉(xiāng)鎮(zhèn)的面積，鄉(xiāng)鎮(zhèn)控制各個(gè)村(村級(jí)調(diào)查區(qū))，村再控制內(nèi)部的各個(gè)圖斑。面積平差的總體流程為：① 計(jì)算各個(gè)鄉(xiāng)鎮(zhèn)的橢球面積；② 匯總各個(gè)鄉(xiāng)鎮(zhèn)面積，得到匯總面積與縣級(jí)控制面積的差值，再將差值以鄉(xiāng)鎮(zhèn)面積大小為權(quán)分配到各個(gè)鄉(xiāng)鎮(zhèn)中；③ 計(jì)算各個(gè)村(村級(jí)調(diào)查區(qū))的橢球面積；④ 匯總各個(gè)村的面積，得到匯總面積與鎮(zhèn)控制面積的差值，再將差值以村面積大小為權(quán)分配到各個(gè)村中；⑤ 計(jì)算各個(gè)圖斑的橢球面積；⑥ 以村為單位進(jìn)行圖斑橢球面積匯總，得到匯總面積與村控制面積的差值，再將差值以圖斑面積大小為權(quán)分配到各個(gè)圖斑中。

4 工具設(shè)計(jì)

橢球面積計(jì)算工具分為3個(gè)模塊來實(shí)現(xiàn)，分別為數(shù)據(jù)管理模塊、橢球面積計(jì)算模塊和面積平差模塊。3個(gè)模塊相互調(diào)用，共同實(shí)現(xiàn)圖斑橢球面積計(jì)算的全部功能。

4.1 數(shù)據(jù)管理模塊

三調(diào)數(shù)據(jù)庫是一個(gè)空間數(shù)據(jù)庫中，行政區(qū)(鎮(zhèn))、村級(jí)調(diào)查區(qū)與地類圖斑分別保存于XZQ、CJDCQ、DLTB層中，數(shù)據(jù)結(jié)構(gòu)如表2-表4所示[2]。

表2 XZQ屬性結(jié)構(gòu)表(部分)

表3 CJDCQ屬性結(jié)構(gòu)表(部分)

表4 DLTB屬性結(jié)構(gòu)表(部分)

此模塊主要功能有：行政區(qū)、村級(jí)調(diào)查區(qū)與地類圖斑數(shù)據(jù)的加載、讀取、賦值，要素類屬性結(jié)構(gòu)檢查，要素的遍歷、空間參考獲取等操作。使用IWorkspaceFactory.Open方法將三調(diào)數(shù)據(jù)庫加載到內(nèi)存，得到Workspace對(duì)象，Workspace對(duì)象同時(shí)實(shí)現(xiàn)了IFeatureWorkspace。使用IFeatureWorkspace. OpenFeatureClass方法代碼行政區(qū)或地類圖斑的名稱，得到FeatureClass對(duì)象。然后通過IFeatureClass.Search方法獲取IFeatureCursor對(duì)象后，就可以利用IFeatureCursor.NextFeature方法實(shí)現(xiàn)對(duì)行政區(qū)、地類圖斑層中各個(gè)要素的遍歷[3]。

4.2 橢球面積計(jì)算模塊

行政區(qū)、村級(jí)調(diào)查區(qū)與地類圖斑在空間數(shù)據(jù)庫中都是以多邊形(Polygon)的形式存在，可以使用統(tǒng)一的方法計(jì)算其橢球面積。在GeoDatabase數(shù)據(jù)模型中，多邊形是由一個(gè)或多個(gè)環(huán)(Ring)組成的。環(huán)可以分為外環(huán)(ExteriorRing)和內(nèi)環(huán)(InteriorRing)，一個(gè)外環(huán)內(nèi)可以有0個(gè)或多個(gè)內(nèi)環(huán)，內(nèi)環(huán)不能獨(dú)立存在。外環(huán)的面積為正數(shù)，內(nèi)環(huán)的面積為負(fù)數(shù)。對(duì)于沒有島洞的多邊形而言，其只有一個(gè)環(huán)[4-5]。

橢球面積計(jì)算模塊主要功能是計(jì)算環(huán)的橢球面積及多邊形內(nèi)部環(huán)的遍歷。計(jì)算環(huán)的橢球面積是此模塊中的重點(diǎn)，也是整個(gè)工具中最難實(shí)現(xiàn)的部分。

對(duì)于一個(gè)環(huán)而言，其主要接口為IRing,但環(huán)對(duì)象也實(shí)現(xiàn)了IPointCollection接口，通過對(duì)環(huán)中每個(gè)節(jié)點(diǎn)的遍歷，每兩個(gè)連續(xù)點(diǎn)組成一個(gè)線段，實(shí)現(xiàn)橢球面積的計(jì)算，主要實(shí)現(xiàn)方法為：① 初始化節(jié)點(diǎn)索引N=0 ；② 通過IPointCollection.Point方法，代入索引號(hào)得到邊 Edeg(N,N+1)；③ 計(jì)算Edeg(N,N+1)的長度，如果小于或等于70 m，則直接使用，如果大于70 m，則進(jìn)行內(nèi)插，形成一個(gè)或多個(gè)子邊；④ 對(duì)于每個(gè)子邊而言，將首尾結(jié)點(diǎn)通過IPoint.Project方法計(jì)算點(diǎn)的經(jīng)緯度坐標(biāo)；⑤ 通過梯形面積公式(1)計(jì)算每個(gè)梯形的面積；⑥ 使用N=N+1，回到步驟2再進(jìn)行計(jì)算，直到計(jì)算到最后一個(gè)點(diǎn)；⑦ 使用公式(2)中提供的方法，匯總所有梯形的面積，得到環(huán)的面積。

XZQ層的橢球面積保存于JSMJ字段中，CJDCQ層的橢球面積保存于JSMJ字段中，DLTB層的橢球面積保存于TBMJ字段中。

4.3 面積平差模塊

面積平差關(guān)系到整個(gè)數(shù)據(jù)庫的面積匯總值能否與國家下發(fā)的控制面積保持一致。受計(jì)算策略、計(jì)算精度及小位取位的影響，直接計(jì)算得到的面積很難與控制面積保持一致，因此得到圖斑的橢球面積后，必須要進(jìn)行面積平差。進(jìn)行平差時(shí)，分配多出面積的基本策略就以橢球面積值為權(quán)進(jìn)行分配，橢球面積值越大，分配到的誤差值就越大。關(guān)鍵算法如下：

(1)匯總XZQ層的JSMJ字段，得到總控制面積與匯總面積的差值，將差值以JSMJ值為權(quán)進(jìn)行分配，平差后的面積保存于DCMJ字段。

(2)一個(gè)XZQ內(nèi)部有一個(gè)或多個(gè)CJDCQ，以XZQ為單位，分組匯總CJDCQ的JSMJ字段，得到每組的匯總值與相應(yīng)XZQ的DCMJ值之間的差值，再將差值以CJDCQ的JSMJ值為權(quán)進(jìn)行分配，平差后的面積值保存于CJDCQ的DCMJ字段。

(3)一個(gè)CJDCQ內(nèi)部有一個(gè)或多個(gè)DLTB，以CJDCQ為單位，分組匯總DLTB的TBMJ，得到每組的匯總值與相應(yīng)CJDCQ的DCMJ值之間的差值，再將差值以DLTB的TBMJ值為權(quán)進(jìn)行分配，平差后的面積值保存于DLTB的TBMJ字段。

在進(jìn)行面積平差時(shí)，受小數(shù)取位的影響，可能一次分配不能達(dá)到匯總值與控制值完全相同，還需要進(jìn)行第二次分配，再次分配還如果還有差值，那么就將這個(gè)差值直接分配給面積最大的那個(gè)圖形。

5 功能實(shí)現(xiàn)

橢球面積計(jì)算工具的開發(fā)基于ArcGIS Add-in框架，IDE選擇Visual Studio2013，使用VisualBasic.NET作為編程語言實(shí)現(xiàn)以上的全部功能。

所有代碼編寫、調(diào)試完成后，即可對(duì)此項(xiàng)目進(jìn)行編譯。編譯通過后，會(huì)生成一個(gè)后綴為esriAddIn的安裝包文件。將此文件拷到具有相應(yīng)ArcGIS版本的電腦上，雙擊此文件后在彈出的窗口中點(diǎn)擊“安裝”完成工具部署。打開ArcMap軟件，在工具欄中將此工具條顯示出來即可使用。圖3為橢球面積計(jì)算工具的操作界面。

第三次國土調(diào)查已經(jīng)在全國范圍內(nèi)開展，橢球面積計(jì)算工具已經(jīng)應(yīng)用到了我院承接了多個(gè)區(qū)縣的三調(diào)項(xiàng)目中，使用此工具能夠快速對(duì)縣級(jí)數(shù)據(jù)庫進(jìn)行面積計(jì)算，計(jì)算結(jié)果能夠通過國家質(zhì)檢軟件的檢查。

圖3 橢球面積計(jì)算

6 結(jié) 語

本文通過橢球面積計(jì)算方法的分析，并以ArcMap AddIn為開發(fā)框架，開發(fā)了專用的橢球面積計(jì)算工具，能夠在ArcMap中快速進(jìn)行橢球面積計(jì)算，方便內(nèi)業(yè)人員的建庫工作。文中關(guān)于ArcMap Add-in開發(fā)的案例對(duì)開發(fā)其他基于ArcMap的應(yīng)用也有一定的借鑒意義。