閭海洋，湯 玲，王 剛

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

0 引 言

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

1 面積計算公式解析

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

(1)

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

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

表1 基本常數(shù)表

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

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

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

SABCDE=SABB1A1+SBCC1B1+SCDD1C11+SDEE1D1+SEAA1E1

(2)

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

2 長邊問題處理

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

圖2 插值計算流程

3 面積平差方法

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

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

4 工具設(shè)計

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

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

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

4.2 橢球面積計算模塊

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

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

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

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

4.3 面積平差模塊

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

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

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

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

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

5 功能實現(xiàn)

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

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

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

圖3 橢球面積計算

6 結(jié) 語

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