橢球面三角形外心的地圖代數(shù)解法

蔣會平，譚樹東，胡　海

1. 武漢大學資源與環(huán)境科學學院，湖北 武漢 430079； 2. 國家海洋信息中心，天津 300171； 3. 中國科學院地理科學與資源研究所資源與環(huán)境信息系統(tǒng)國家重點實驗室，北京 100101； 4. 中國科學院大學，北京 100049

?

橢球面三角形外心的地圖代數(shù)解法

蔣會平1,3,4，譚樹東2，胡海1

1. 武漢大學資源與環(huán)境科學學院，湖北 武漢 430079； 2. 國家海洋信息中心，天津 300171； 3. 中國科學院地理科學與資源研究所資源與環(huán)境信息系統(tǒng)國家重點實驗室，北京 100101； 4. 中國科學院大學，北京 100049

Foundation support： The National Natural Science Foundation of China (Nos. 41271443；41471328); The National High-tech Research and Development Program of China (863 Program) (No. 2009AA12Z224)

摘要：橢球面三角形外心到3個相鄰頂點的大地線距離都相等。面向橢球面空間的外心大地坐標的求解對于橢球面Voronoi圖的生成和橢球面Delaunay三角網(wǎng)的構(gòu)造具有重要作用。利用基于地圖代數(shù)理論的矢柵結(jié)合方法，首先基于地圖代數(shù)測地變換建立高精度橢球面空間距離場，再通過邊界跟蹤配對確定外心所在的柵格范圍，最后通過數(shù)值計算內(nèi)插生成初始等距點并不斷逼近外心的精確大地坐標。試驗結(jié)果表明，采用本文方法求解的橢球面三角形外心大地坐標，在103～104km跨度內(nèi)其定位誤差小于0.001 m，且算法非常適用于海量空間數(shù)據(jù)的高精度快速計算。

關(guān)鍵詞：地圖代數(shù)；測地變換；橢球面；三角形外心；矢柵結(jié)合

從計算幾何[1]的觀點來看，橢球面三角形外心就是與3個相鄰頂點所對應(yīng)的橢球面Voronoi圖的公共頂點，而Voronoi圖具有許多優(yōu)秀的性質(zhì)，研究橢球面Voronoi圖對于全球數(shù)據(jù)的管理和空間關(guān)系的動態(tài)維護具有重要意義[2]。根據(jù)Voronoi圖和Delaunay三角網(wǎng)的對偶關(guān)系，采用間接Delaunay三角剖分算法由橢球面Voronoi圖能夠比較容易地構(gòu)造橢球面Delaunay三角網(wǎng)，而Delaunay三角網(wǎng)被廣泛應(yīng)用于數(shù)字高程模型(DEM)構(gòu)建、空中三角測量、大地控制網(wǎng)評估和連續(xù)運行參考系統(tǒng)(CORS)網(wǎng)布站優(yōu)化[3]。因此作為生成橢球面Voronoi圖并構(gòu)造相應(yīng)的Delaunay三角網(wǎng)的關(guān)鍵步驟，研究橢球面三角形外心的求解具有重要的理論和應(yīng)用價值。

在平面上由3個相鄰頂點計算三角形外心的坐標可采用確定的公式和成熟的算法[4]，然而大地線距離是由在橢球面某范圍內(nèi)足夠精確且不產(chǎn)生二義的大地線尺度函數(shù)所定義的，直接在橢球面上計算三角形外心大地坐標異常復(fù)雜[5]。因此在傳統(tǒng)的地理信息系統(tǒng)中大多數(shù)采用球面或投影平面對橢球面進行高或低精度的近似以降低實現(xiàn)難度[6-7]，例如在球面上103km跨度內(nèi)誤差達到了35 m，而在投影平面上103km跨度內(nèi)誤差竟達到了8 km[7]，可見兩者的誤差控制情況都不是很理想。

文獻[4]以構(gòu)建區(qū)域性橢球面數(shù)字地面模型(DTM)為目的，提出了一種在測地坐標系[8]中由橢球面三角形3個相鄰頂點求解其外心測地坐標的方法。對于平均邊長在1～2 km的橢球面三角形，計算所得外心到該三角形各頂點的大地線距離偏差小于0.1 mm[4]。因為測地坐標系是基于區(qū)域性橢球面建立起來的，適用范圍僅限于較大的局部區(qū)域[5]，所以該方法對于跨度遠遠超過2 km甚至達到103～104km的很多實際問題，如無人機航跡規(guī)劃[9]、巡航導(dǎo)彈突防路徑規(guī)劃[10]以及海域劃界中間線及其三接點的求解[11]等，無法保證其結(jié)果的高精度。

文獻[12]以求解海域劃界中間線上的轉(zhuǎn)折點位置為目的，提出了一種在球面和橢球面上求解三角形外心大地坐標的初步思路，但是未見該方法的有效試驗結(jié)果。對于橢球面三角形外心大地坐標的求解，首先要計算得到相應(yīng)的球面三角形外心的經(jīng)緯度坐標，然后采用求解線性方程組的一般迭代方法逼近外心的精確大地坐標，然而迭代過程中產(chǎn)生的向量序列并不一定逐步逼近線性方程組的解，且不同迭代方法的收斂效率差別較大[13]。

為了能夠在全球范圍內(nèi)高效計算橢球面三角形外心的精確大地坐標，解決上述方法各自存在的問題，本文基于地圖代數(shù)[14]中的測地變換[15]，提出了一種求解橢球面三角形外心的矢柵結(jié)合方法。算法可以直接在橢球面上快速計算大量離散點集的高精度Voronoi頂點，對于橢球面Voronoi圖的生成和橢球面Delaunay三角網(wǎng)的構(gòu)造具有重要作用，從而為建設(shè)動態(tài)全球性地理信息系統(tǒng)[16-17](Global GIS)并實現(xiàn)精確橢球幾何計算[18]與空間分析[19]打下基礎(chǔ)。

1方法原理

本文方法基于地圖代數(shù)理論，利用“矢柵結(jié)合使用長處”[20]的思想，通過柵格大地線尺度距離變換——測地變換，先精確求解橢球面空間中各柵格中心點到其最近地理實體上來源點的距離，后在臨界相鄰柵格范圍內(nèi)，以柵格中心點作為控制點，通過數(shù)值計算內(nèi)插初始等距點大地坐標并最終逼近橢球面三角形外心的精確位置。算法的誤差主要來自于在橢球面3～4個臨界相鄰柵格范圍內(nèi)由初始等距點逼近三角形外心所導(dǎo)致的精度損失，下面從橢球大地測量學[21]和球面三角學[22]的角度論證：在滿足一定條件下上述逼近過程與大地主題解算的精度[23]能夠保持一致。

對于橢球面三角形ABC，由橢球面正弦定理[21]可得

(1)

圖1　橢球面三角形ABC與球面三角形A0B0C0Fig.1　Triangle ABC on ellipsoidal surface and triangle A0B0C0on spherical surface

對于球面三角形A0B0C0，由球面正弦定理[28]可得

(2)

比較式(1)和式(2)可得

(3)

(4)

同理可得

(5)

(6)

式(6)為橢球面改正數(shù)公式，當橢球面三角形ABC的平均邊長小于200 km時，橢球面改正數(shù)小于0.001″，而當200 km線長的方向變化0.001″時，端點的位移小于0.001 m。因此在平均邊長小于200 km的條件下可直接把橢球面三角形當作球面三角形，而對于平均邊長在200～400 km之間的橢球面三角形需要經(jīng)過適當改正轉(zhuǎn)化為球面三角形[21]。

如圖2所示，橢球面三角形DEF和DEG共用同一邊DE，且大地線DF和DG長度相等，∠EDF和∠EDG也相等。對應(yīng)到本文方法中，D表示初始等距點，F(xiàn)和G分別表示相鄰兩個柵格中心點的來源點，E表示∠FDG平分線(大地線)上的一點。當DF和DG的長度不超過104km且DE的長度不大于1 km時，如果能夠論證EF和EG的長度偏差小于0.001 m，由初始等距點逼近三角形外心與大地主題解算的精度[23]相當?shù)慕Y(jié)論就是成立的。在密切球上做與橢球面三角形DEF和DEG對應(yīng)邊長度相等的三角形D0E0F0和D0E0G0，已知大圓弧D0F0和D0G0長度相等，假設(shè)∠E0D0F0和∠E0D0G0也相等，即兩邊及其夾角彼此相等。根據(jù)球面三角形全等定理[22]，球面三角形D0E0F0和D0E0G0全等，由全等球面三角形的性質(zhì)可知大圓弧E0F0和E0G0的長度也相等。實際上，將橢球面三角形轉(zhuǎn)化到對應(yīng)球面上以后，∠E0D0F0和∠E0D0G0并不相等，因此大圓弧E0F0和E0G0的長度不相等。

圖2　橢球面三角形DEF和DEG與球面三角形D0E0F0和D0E0G0Fig.2　Triangle DEF and DEG on ellipsoidal surface and triangle D0E0F0and D0E0G0on sphericalsurface

由式(3)可知，橢球面上的∠EDF與對應(yīng)球面上的∠E0D0F0僅相差一個橢球面改正數(shù)，其大小主要與橢球面上大地線DE和DF長度的乘積有關(guān)。因為大地線DF的長度不超過104km，且大地線DE的長度不大于1 km，所以橢球面改正數(shù)小于0.001″，由球面余弦公式[22]可知對應(yīng)球面上大圓弧E0F0端點的位移小于0.001 m，同理可知大圓弧E0G0端點的位移也小于0.001 m，即大地線EF與EG的長度偏差在0.001 m以內(nèi)可控，也即逼近過程與大地主題解算的精度[23]保持一致。

一般情況下大地線DF、DG和DE的長度越短，大地線EF與EG的長度偏差就越小，逼近的精度也就越高。在本文方法中DF、DG的大地線長度對應(yīng)于外心到橢球面三角形頂點的大地線距離，DE的大地線長度對應(yīng)于柵格單元的跨度。因此理論上算法的精度與大地線距離的長短呈負相關(guān)，與柵格分辨率的高低呈正相關(guān)。

2關(guān)鍵技術(shù)

本文方法基于地圖代數(shù)測地變換，充分發(fā)揮了矢柵緊密結(jié)合數(shù)據(jù)結(jié)構(gòu)的優(yōu)勢，在保證柵格計算效率的條件下，能夠計算定位誤差小于0.001 m的大跨度(103～104km范圍)內(nèi)橢球面三角形外心大地坐標。按照算法流程的先后順序，主要涉及以下幾個關(guān)鍵技術(shù)：測地變換、邊界跟蹤配對和求解外心大地坐標,其中邊界跟蹤配對用來確定外心所在的臨界相鄰柵格范圍。

2.1測地變換

所謂測地變換是指在橢球面上進行的柵格大地線尺度距離變換，其目的是建立橢球面空間的高精度距離場。地圖代數(shù)測地變換的時間復(fù)雜度為O(m)，其中m為計算規(guī)模，是全空間的柵格單元總數(shù)[15]，等于橫向分辨率(地理網(wǎng)格陣列的列數(shù))與縱向分辨率(地理網(wǎng)格陣列的行數(shù))的乘積。與利用窮舉算法進行距離變換相比，地圖代數(shù)測地變換完全不受生成元個數(shù)和生成元復(fù)雜程度的影響，不但能夠適用于復(fù)雜的自然圖形，而且可以大大提高計算效率，適用于本文問題的解決。作為大地線尺度下的空間度量，測地變換的精度最終取決于大地主題解算公式的精度，在20 000 km的范圍內(nèi)誤差小于數(shù)厘米[23]，因此可以認為柵格中心間的大地線距離是沒有偏差的，在全球范圍內(nèi)可用。

文獻[15]規(guī)定測地變換必須滿足3個基礎(chǔ)定義，其中定義3明確指出兩定點a、b間的柵格距離等同于a、b所屬柵格中心間的距離，顧及整個柵格和舍入誤差，有不大于0.5個柵格單元的長度[15]。雖然該算法的精度可以由柵格分辨率來控制，理論上能夠達到很高的精度，但是對于某些實際應(yīng)用，例如要在103～104km范圍內(nèi)計算誤差小于0.001 m的橢球面三角形外心，各柵格單元的跨度都不能超過0.001 m，相應(yīng)的橫向和縱向分辨率將達到109～1010，全空間的柵格單元總數(shù)就增加到了1018～1020，現(xiàn)階段普通的計算機顯然無法處理規(guī)模如此龐大的地理網(wǎng)格陣列。

2.2邊界跟蹤配對

文獻[24]提出了一種基于八鄰域邊界跟蹤的物體標號算法，一次掃描即可同時獲得物體的邊界點序列以及邊界鏈碼，為后續(xù)處理提供數(shù)據(jù)準備。在此基礎(chǔ)上，本文通過以下方法對測地變換所得分配場進行邊界跟蹤配對：首先找到最左上角的一個邊界柵格作為搜索起點，按順時針方向，自上而下、從左至右，搜索其八鄰域，并將搜索路徑上的四鄰域柵格作為上一個邊界柵格的配對柵格，直到找到下一個邊界柵格；然后以此邊界柵格為新的起點繼續(xù)搜索，這一搜索過程不斷重復(fù)下去，直至回到搜索起點。經(jīng)過邊界跟蹤配對以后，所有邊界柵格的鄰接匹配關(guān)系都可以被確定。遍歷所得邊界柵格對序列，按順序從前往后每兩對一組進行判別掃描，過濾出匹配關(guān)系發(fā)生變化的相鄰兩對邊界柵格對，即為橢球面三角形外心所在的柵格范圍。

2.3求解外心大地坐標

經(jīng)過邊界跟蹤配對確定橢球面三角形外心所在的柵格范圍以后，最終求解足夠精確的結(jié)果又可分為兩個步驟：初始等距點的內(nèi)插以及外心大地坐標的逼近，它們分別如圖3和圖4所示。其中O1和O2分別表示上述相鄰兩對邊界柵格對中的配對柵格中心點，B1和B2分別表示O1和O2的來源點，B1到O1的大地線距離和B2到O2的大地線距離分別等于S1和S2，O3和O4分別表示與O1和O2相互對應(yīng)的邊界柵格中心點，B3表示O3和O4的同一個來源點，B3到O3的大地線距離和B3到O4的大地線距離分別等于S3和S4。

2.3.1內(nèi)插初始等距點

經(jīng)過大地線尺度下的地圖代數(shù)測地距離變換，可以確保O1距離B1比B2更近，O2距離B2比B1更近。由大地線長度變化的連續(xù)性容易證明，在如圖3所示的O1與O2連接成的大地線段上必然存在一點O到B1和B2的大地線距離都等于S，即得到初始等距點O。

2.3.2逼近橢球面三角形外心

在103～104km范圍內(nèi)，當柵格單元的跨度不大于1 km時，在如圖4所示的3～4個臨界相鄰柵格范圍內(nèi)，∠α的平分線(大地線)上任意一點P到B1和B2的大地線距離偏差都小于0.001 m。由大地線長度變化的連續(xù)性同理可證，在該角平分線上必然存在一個到B1和B3的大地線距離相等或到B2和B3的大地線距離相等的點，不妨假設(shè)該點為P，則P到B1、B2和B3的大地線距離S′、S″和S?的偏差都小于0.001 m，即求解得到由B1、B2和B3構(gòu)成的橢球面三角形外心的精確位置。

圖3　內(nèi)插初始等距點Fig.3　Interpolation of the initial equidistant point

圖4　逼近橢球面三角形外心Fig.4　Approximation of the circumcenter of triangle on ellipsoidal surface

3具體實現(xiàn)及試驗結(jié)果分析

3.1具體實現(xiàn)

充分利用測地變換所生成的高精度距離場，通過數(shù)值計算先內(nèi)插初始等距點，再以同樣的方式不斷逼近橢球面三角形外心是本文方法的關(guān)鍵。原則上在103～104km范圍內(nèi)要使得各柵格單元的跨度都不超過1km，因此一般情況下將橫向分辨率范圍設(shè)定為103～104，縱向分辨率范圍也隨之確定，相應(yīng)的計算規(guī)模(全空間的柵格單元總數(shù))控制在106～108之間。本文方法的具體實現(xiàn)步驟如圖5所示，其中源點坐標場記錄的是各控制點的最近來源點經(jīng)緯度坐標信息，距離場記錄的是每個控制點到對應(yīng)來源點的大地線距離信息，分配場則記錄了所有控制點與其最近地理實體之間的對應(yīng)關(guān)系。

3.2試驗結(jié)果分析

為了保證方法在全球范圍內(nèi)的普遍適用性，本文采用了WGS-84坐標中的橢球基本參數(shù)：橢球長半軸a= 6 378 137.0m，橢球扁率f= 1.0/298.257 223 563，且采用了解算精度與距離無關(guān)的貝塞爾大地主題解算公式[25]。

圖5　橢球面三角形外心的求解流程Fig.5　Flow chart of determination of circumcenter of triangle on ellipsoidal surface

3.2.1驗證方法精度與大地線距離的關(guān)系

進行以下10組數(shù)據(jù)的計算，其中橫向分辨率都設(shè)定在5000以內(nèi)，各橢球面三角形3個頂點及其外心的位置如圖6所示，具體的大地坐標見表1。為了驗證所得結(jié)果的正確性，分別計算各外心到對應(yīng)橢球面三角形3個頂點的大地線距離及其標準差，詳細結(jié)果見表2。在表1和表2中，λ都表示經(jīng)度，φ則表示緯度。

由表2可知，隨著外心到對應(yīng)橢球面三角形各頂點的大地線距離變得越來越遠，標準差呈現(xiàn)不斷增大的趨勢，說明本文方法的計算精度與大地線距離的長短呈負相關(guān)。1—9號橢球面三角形外心到對應(yīng)各三角形3個頂點的大地線距離之差都小于0.001m，其中外心到對應(yīng)三角形頂點的最短大地線距離為307 648.857 7m，標準差不超過0.000 06m；最長大地線距離達到3 827 513.452 9m，標準差也不超過0.000 5m。10號橢球面三角形外心到對應(yīng)各三角形3個頂點的大地線距離分別是5 646 875.751 3m、5 646 875.751 8m和5 646 875.753 9m，雖然兩兩之差不都小于0.001m，但是標準差不超過0.002m。因此在中、長距離情況下，本文方法對于橢球面三角形外心的求解精度非常高。在超長距離情況下，計算精度雖有所降低，但仍然保持在0.001m左右。

圖6　橢球面上各三角形及其外心與外接大地圓Fig.6　Triangles on ellipsoidal surface and their geodetic circumcenters and circumcircles

3.2.2驗證方法精度與柵格分辨率的關(guān)系

進行以下5組數(shù)據(jù)的計算，其中構(gòu)成橢球面

三角形ABC各頂點的大地坐標分別是A(88.131 4°,55.252 6°)、B(100.333 3°,40.515 2°)和C(111.666 6 °,66.373 8 °)，分別計算不同分辨率下外心的大地坐標和相應(yīng)外心到各頂點的大地線距離及其標準差，詳細結(jié)果見表3，在表3中λ和φ分別表示經(jīng)度和緯度。

由表3可知，隨著分辨率變得越來越高，標準差呈現(xiàn)不斷減小的趨勢，說明本文方法的計算精度與分辨率的高低呈正相關(guān)。當橫向分辨率從1000逐漸增大到4000時，標準差從0.001 556 706 0m逐漸減小到0.000 208 166 6m，計算結(jié)果的精度在不斷提高。當橫向分辨率繼續(xù)增大到5000時，標準差反而增加到0.000 321 455 1m，說明當分辨率增加到一定程度后，對于算法精度的提高作用是有限的，由于大地主題解算的精度限制[23]，此時大地主題解算的誤差對結(jié)果的影響超過逼近誤差，精度發(fā)生異常波動屬于正?，F(xiàn)象。因此本文方法能夠在可計算的柵格分辨率范圍內(nèi)(橫向分辨率在103～104之間，全空間的柵格單元總數(shù)在106～108之間)高效計算定位誤差小于0.001m的大跨度(103～104km范圍內(nèi))橢球面三角形外心大地坐標。

表1　橢球面三角形各頂點及其外心的大地坐標

3.2.3分析橢球面三角形各頂點所處緯度的高低對算法精度的影響

由于本文方法在劃分網(wǎng)格單元的步驟中采用的是基于大地坐標系統(tǒng)的傳統(tǒng)平面格網(wǎng)系統(tǒng)構(gòu)建方法，不可避免地存在面積變形、形狀變形以及由赤道向南北極遞增的問題[26]，因此分析方法精度與橢球面三角形各頂點所處緯度的關(guān)系也很有必要。觀察表2中的6號和7號橢球面三角形，兩者的橫向分辨率都設(shè)定為4000，前者的外心到各頂點的大地線距離分別是1 680 030.998 3 m、1 680 030.998 3 m和1 680 030.999 1 m，標準差是0.000 461 880 1 m，后者的外心到各頂點的大地線距離分別是2 092 296.495 1 m、2 092 296.494 2 m和2 092 296.494 6 m，標準差是0.000 450 924 9 m，同屬中長距離情況，計算精度基本相當。由表1可知，6號橢球面三角形各頂點處于低緯度地區(qū)，而7號橢球面三角形各頂點處于中高緯度地區(qū)。因此在正確劃分網(wǎng)格單元的條件下，本文方法的精度與所處緯度無關(guān)。

表2　外心到各橢球面三角形3個頂點的大地線距離及其標準差

表3不同分辨率下橢球面三角形外心的大地坐標與外心到各頂點的大地線距離及其標準差

Tab.3Geodetic coordinates of the circumcenters of triangles on ellipsoidal surface and geodesic distances from the circumcenters to vertexes of triangles on ellipsoidal surface and related standard deviations in a variety of resolutions

橫向分辨率外心坐標/(°)外心到各頂點的大地線距離/mλφ到A的距離到B的距離到C的距離標準差/m1000111.469875441352.33384784731564069.06231564069.06101564069.05920.00155670602000111.469875407452.33384784121564069.06041564069.05911564069.05990.00065574393000111.469875378252.33384784721564069.05831564069.05841564069.05920.00049328834000111.469875388152.33384784731564069.05891564069.05881564069.05920.00020816665000111.469875385952.33384784651564069.05881564069.05871564069.05930.0003214551

3.2.4分析地理數(shù)據(jù)的復(fù)雜程度對算法效率的影響

選取如圖7所示的位于大西洋兩岸的100個城市結(jié)點進行試驗。在橫向分辨率設(shè)定為4000的條件下，采用配置為Intel Core 2 Duo E7500處理器和2 GB內(nèi)存的臺式電腦，本文方法計算100個城市的所有167個Voronoi頂點的整個過程耗時不超過20 s(包括測地變換過程在內(nèi))，與相同條件下計算單個橢球面三角形外心速度相當，說明算法效率基本不受對象個數(shù)的影響。作為一種面向橢球面空間且內(nèi)容無關(guān)的算法，本文方法更適合于海量空間數(shù)據(jù)的高精度快速計算。

綜上所述，本文方法精度與外心到對應(yīng)橢球面三角形各頂點的大地線距離的長短呈負相關(guān)，與柵格分辨率的高低呈正相關(guān)，與橢球面三角形各頂點所處緯度無關(guān)，能夠滿足在全球范圍內(nèi)以可計算的柵格分辨率(全空間的柵格單元總數(shù)在106～108之間)快速求解外心精確大地坐標(定位誤差小于0.001 m)的需求。

圖7　大西洋兩岸100個城市及其Voronoi頂點Fig.7　The hundred cities on both sides of the Atlanticand their Voronoi vertexes

4結(jié)論

本文方法采用矢柵緊密結(jié)合的數(shù)據(jù)結(jié)構(gòu)，基于地圖代數(shù)測地變換建立了橢球面空間的高精度距離場，在對分配場進行邊界跟蹤、配對、遍歷和過濾操作之后，通過數(shù)值計算內(nèi)插逼近的方法，最終計算出精確的大跨度橢球面三角形外心大地坐標。算法充分發(fā)揮了矢柵結(jié)合的優(yōu)勢，具有面向橢球面空間且內(nèi)容無關(guān)的特性，能夠同時兼顧計算精度和效率。該方法可以用來生成橢球面Voronoi圖和構(gòu)造橢球面Delaunay三角網(wǎng)，且已經(jīng)成功應(yīng)用于海域劃界三接點的高精度快速求解。后續(xù)工作將研究如何將本文方法拓展到橢球面緩沖區(qū)和等比例線的求解問題[27]，從而提高海域劃界相關(guān)計算的精度和效率。

參考文獻：

[1]周培德. 計算幾何[M]. 北京: 清華大學出版社, 2000.

ZHOU Peide. Computational Geometry[M]. Beijing: Tsinghua University Press, 2000.

[2]GOLD C M. The Global GIS[C]∥Proceeding of the International Workshop on Dynamic and Multi-Dimension GIS. Hong-Kong, China: [s.n.], 1997: 80-91.

[3]劉廣軍, 郭晶, 劉旭東. 基于橢球面大地線三角網(wǎng)的CORS布站研究與仿真[C]∥第13屆中國系統(tǒng)仿真技術(shù)及其應(yīng)用學術(shù)年會論文集. 黃山: 中國自動化學會系統(tǒng)仿真專業(yè)委員會, 2011.

LIU Guangjun, GUO Jing, LIU Xudong. Research and Simulation on CORS Station Distribution Based on Ellipsoid Geodetic-line Triangulation Network[C]∥Proceedings of 13th Chinese Conference on System Simulation Technology & Application. Huangshan: Professional Committee of China Automation Society System Simulation, 2011.

[4]施一民, 朱紫陽. 橢球面上Delaunay三角形的外接大地圓圓心的求解[J]. 同濟大學學報, 2004, 32(1): 78-81.

SHI Yimin, ZHU Ziyang. Determination of Center for Geodetic Circle Which Passes through Three Top Points of Triangle on Ellipsoidal Surface[J]. Journal of Tongji University, 2004, 32(1): 78-81.

[5]馮琰, 施一民. 基于區(qū)域性橢球面數(shù)字地面模型的研究[J]. 同濟大學學報, 2003, 31(8): 964-967.

FENG Yan, SHI Yimin. New DTM Based on Surface of Regional Ellipsoid[J]. Journal of Tongji University, 2003, 31(8): 964-967.

[6]趙學勝, 陳軍, 王金莊. 基于O-QTM的球面Voronoi圖的生成算法[J]. 測繪學報, 2002, 31(2): 157-163.

ZHAO Xuesheng, CHEN Jun, WANG Jinzhuang. QTM-based Algorithm for the Generating of Voronoi Diagram for Spherical Objects[J]. Acta Geodaetica et Cartographica Sinica, 2002, 31(2): 157-163.

[7]HOREMUZ M. Error Calculation in Maritime Delimitation between States with Opposite or Adjacent Coasts[J]. Marine Geodesy, 1999, 22(1): 1-17.

[8]施一民, 馮琰. 地球橢球面上另一種形式的測地坐標系的建立[J]. 同濟大學學報, 2001, 29(11): 1282-1285.

SHI Yimin, FENG Yan. Establishment of Another Form of Geodesic Coordinate System on Earth Ellipsoid[J]. Journal of Tongji University, 2001, 29(11): 1282-1285.

[9]趙文婷, 彭俊毅. 基于Voronoi圖的無人機航跡規(guī)劃[J]. 系統(tǒng)仿真學報, 2008(z2): 159-162, 165.

ZHAO Wenting, PENG Junyi. Voronoi Diagram-based Path Planning for UAVs[J]. Journal of System Simulation, 2008(z2): 159-162, 165.

[10]閻代維, 谷良賢, 王興治. 基于Voronoi圖的巡航導(dǎo)彈突防路徑規(guī)劃研究[J]. 彈箭與制導(dǎo)學報, 2005, 25(2): 11-13.

YAN Daiwei, GU Liangxian, WANG Xingzhi. The Study of Cruise Missile Path Planning with Voronoi Diagram[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2005, 25(2): 11-13.

[11]FERRERO S, PIEROZZI M, REPETTI L, et al. An Algorithm for the Unambiguous Determination of the Equidistant Boundary Line between Two (or More) Coastlines[J]. Applied Geomatics, 2009, 1(3): 49-58.

[12]SJ?BERG L E. The Three-point Problem of The Median Line Turning Point: on the Solutions for the Sphere and Ellipsoid[J]. International Hydrographic Review, 2002, 3(1): 81-87.

[13]李慶揚, 王能超, 易大義. 數(shù)值分析[M]. 第5版. 北京: 清華大學出版社, 2008.

LI Qingyang, WANG Nengchao, YI Dayi. Numerical Analysis[M]. 5th ed. Beijing: Tsinghua University Press, 2008.

[14]胡鵬, 游漣, 胡海. 地圖代數(shù)概論[M]. 北京: 測繪出版社, 2008.

HU Peng, YOU Lian, HU Hai. Introduction to Map Algebra[M]. Beijing: Publishing House of Surveying and Mapping, 2008.

[15]胡鵬, 范青松, 胡海. 橢球上的測地變換和Voronoi圖的生成——地理空間度量[J]. 武漢大學學報(信息科學版), 2007, 32(9): 825-828.

HU Peng, FAN Qingsong, HU Hai. Distance Transformation and Voronoi Generation on Earth Ellipsoid——Metrics of Geographic Space[J]. Geomatics and Information Science of Wuhan University, 2007, 32(9): 825-828.

[16]GOLD C, MOSTAFAVI M A. Towards the Global GIS[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2000, 55(3): 150-163.

[17]胡鵬, 劉沛蘭, 胡海, 等. 地球信息的度量空間和Global GIS[J]. 武漢大學學報(信息科學版), 2005, 30(4): 317-321.

HU Peng, LIU Peilan, HU Hai, et al. Metric Space of Earth Information and Global GIS[J]. Geomatics and Information Science of Wuhan University, 2005, 30(4): 317-321.

[18]KJENSTAD K. Construction and Computation of Geometries on the Ellipsoid[J]. International Journal of Geographical Information Science, 2011, 25(9): 1413-1437.

[19]HU Hai, LIU Xiaohang, HU Peng. Voronoi Diagram Generation on the Ellipsoidal Earth[J]. Computers & Geosciences, 2014, 73: 81-87.

[20]胡鵬, 楊傳勇, 胡海, 等. GIS的基本理論問題——地圖代數(shù)的空間觀[J]. 武漢大學學報(信息科學版), 2002, 27(6): 616-621.

HU Peng, YANG Chuanyong, HU Hai,et al. Space View of Map Algebra[J]. Geomatics and Information Science of Wuhan University, 2002, 27(6): 616-621.

[21]熊介. 橢球大地測量學[M]. 北京: 解放軍出版社, 1988.

XIONG Jie. Ellipsoidal Geodesy[M]. Beijing: Publishing House of Liberation Army, 1988.

[22]張楚賓. 球面三角學[M]. 北京: 高等教育出版社, 1959.

ZHANG Chubin. Spherical Trigonometry[M]. Beijing: Higher Education Press, 1959.

[23]陳健, 晁定波. 橢球大地測量學[M]. 北京: 測繪出版社, 1992.

CHEN Jian, CHAO Dingbo. Ellipsoidal Geodesy[M]. Beijing: Publishing House of Surveying and Mapping, 1992.

[24]劉相濱, 向堅持, 陽波. 基于八鄰域邊界跟蹤的標號算法[J]. 計算機工程與應(yīng)用, 2001, 37(23): 125-126, 132.

LIU Xiangbin, XIANG Jianchi, YANG Bo. A Labeling Algorithm Based on 8-Connected Boundary Tracking[J]. Computer Engineering and Applications, 2001, 37(23): 125-126, 132.

[25]孔祥元, 郭際明, 劉宗泉. 大地測量學基礎(chǔ)[M]. 武漢: 武漢大學出版社, 2005.

KONG Xiangyuan, GUO Jiming, LIU Zongquan. Foundation of Geodesy[M]. Wuhan: Wuhan University Press, 2005.

[26]宋樹華, 程承旗, 關(guān)麗, 等. 全球空間數(shù)據(jù)剖分模型分析[J]. 地理與地理信息科學, 2008, 24(4): 11-15.

SONG Shuhua, CHENG Chengqi, GUAN Li, et al. Analysis on Global Geodata Partitioning Models[J]. Geography and Geo-Information Science, 2008, 24(4): 11-15.

[27]胡海, 吳艷蘭. 橢球面上混合基線圖形的緩沖區(qū)和等比例線問題[J]. 測繪工程, 2011, 20(3): 1-4, 8.

HU Hai, WU Yanlan. Problems of Graphics Buffers and the Proportion Line of Ellipsoid Mixed Baseline (or Normal Curve)[J]. Engineering of Surveying and Mapping, 2011, 20(3): 1-4, 8.

(責任編輯：張艷玲)

First author： JIANG Huiping(1990—), male, postgraduate, majors in map algebra, spatial analysis on ellipsoid and global GIS.

E-mail： jianghuiping@whu.edu.cn; jianghp@lreis.ac.cn

Corresponding author： HU Hai

E-mail： huhai@whu.edu.cn

修回日期： 2015-08-10

Determination of Circumcenter of Triangle on Ellipsoidal Surface Based on Map Algebra

JIANG Huiping1,3,4, TAN Shudong2, HU Hai1

1. School of Resources and Environment Sciences, Wuhan University, Wuhan 430079, China; 2. National Marine Data and Information Service, Tianjin 300171, China; 3. State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101, China; 4. University of Chinese Academy of Sciences, Beijing 100049, China

Abstract：The geodesic distances from the circumcenter to 3 vertexes of the triangle on ellipsoidal surface are equal. The ellipsoid-oriented determination of circumcenter of triangle on ellipsoidal surface is applicable when it comes to generation of the Voronoi diagram and construction of the Delaunay triangulation net on the ellipsoidal earth, which can be considered as a solution of significance in computation of geometries and spatial analysis on the ellipsoid. Based on the idea of combining the raster and vector methods and the theory of map algebra, the working process can be described as below: firstly, initiate the geographical distance transformation and create the distance field with a high degree of accuracy; secondly, conduct boundary tracking and matching and then determinate the range of grids where the circumcenter of triangle locates; thirdly, interpolate the initial equidistant point; finally, approximate the circumcenter of triangle on earth ellipsoidal surface by means of numeric calculation. The positioning error of this algorithm is controlled less than 0.001 m within several thousand kilometers range of span. As regards the method proposed in the present paper, its computational efficiency is O(m) where m is the number of pixels in the image, i.e., grid resolution. In conclusion, this algorithm can be considered as both ellipsoid-oriented and not content-related, which is especially appropriate for complex geocomputation globally.

Key words：map algebra; geographical distance transformation; ellipsoidal surface; circumcenter of triangle; raster&vector-based

通信作者：胡海

作者簡介：第一 蔣會平(1990—)，男，碩士生，主要研究方向為地圖代數(shù)、橢球空間分析以及全球性地理信息系統(tǒng)。

收稿日期：2014-10-14

基金項目：國家自然科學基金 (41271443；41471328)；國家863計劃 (2009AA12Z224)

中圖分類號：P208

文獻標識碼：A

文章編號：1001-1595(2016)02-0241-09

引文格式：蔣會平，譚樹東，胡海.橢球面三角形外心的地圖代數(shù)解法[J].測繪學報，2016,45(2)：241-249.DOI:10.11947/j.AGCS.2016.20140503.

JIANG Huiping, TAN Shudong, HU Hai.Determination of Circumcenter of Triangle on Ellipsoidal Surface Based on Map Algebra[J]. Acta Geodaetica et Cartographica Sinica,2016,45(2):241-249.DOI:10.11947/j.AGCS.2016.20140503.