等差分緯線多圓錐投影正解變換的擬合函數模型及誤差分析

李 家，肖益松

(1. 遼寧省自然地理與空間信息科學實驗室，遼寧 大連 116027； 2. 遼寧師范大學 城市與環境學院，遼寧 大連 116029)

一、投影的基本公式和設計參數

等差分緯線多圓錐投影已在我國編制的各種比例尺的世界地圖及其他類型世界地圖中得到廣泛應用，并獲得了較好的效果[1]。由于在ArcGIS、MapInfo、SuperMap等GIS軟件中沒有該投影的正反解變換功能，在實際工作中只能通過二次開發實現該功能。

根據文獻[1]，計算該投影中央經線上x0、y0的公式為

式中，μ0=10 000 000；R=6 371 116 m；φ為要變換的點的緯度值(以弧度計)。

極坐標公式為

式中，(xn,yn)是邊緣經線上(Δλ=180°)根據所設計的經、緯線草圖所量取的直角坐標值；b=1.1；C=0.000 505 050 5；δφi是某一緯線上各經線的極角；λn是邊緣經線與中央經線的經差；λi為各點經線與中央經線的經差(如圖1所示)。

計算某一點投影后直角坐標的公式為

文獻[1]同時給出了該投影的設計參數，包括直角坐標值(見表1)、緯圈半徑ρ和邊緣經線上的極角δφi的值(見表2)。

圖1 等差分緯線多圓錐投影直角坐標計算示意圖

表1 等差分緯線多圓錐投影的直角坐標值

續表

表2等差分緯線多圓錐投影的緯圈半徑ρ和邊緣經線上的極角δφi的值

φρ/cmδφi10°2399.9433°54′10.5″ 15°1643.3485°38′9.5″ 20°1201.1747°35′46.0″ 回歸線1024.1598°47′52.0″ 30°807.63410°51′26.0″ 40°618.06013°29′8.5″ 45°556.62914°31′58.5″ 50°507.56115°24′15.5″ 60°430.76316°39′4.0″ 極圈388.63317°12′31.5″ 70°375.47817°5′27.0″ 75°372.03016°11′0.5″ 80°388.52514°28′57.5″ 90°508.0629°20′42.5″

表1、表2中的數據是基于1∶10 000 000的世界地圖的，(x，y)坐標和ρ的單位為厘米。

對地圖上某一點從經緯度坐標變換為直角坐標(正解變換)的過程是：根據該點的緯度值和邊緣經線的函數式求出該點所在緯線與邊緣經線的交點(xn,yn)，這樣就可以確定相應的緯圈半徑ρ；根據ρ和yn計算交點(xn，yn)的極角δφn(如圖1所示)；利用式(2)計算該點的極角δφi；最后利用式(3)計算其直角坐標值。

二、邊緣經線的擬合函數模型

國家測繪地理信息局在網上正式發布的世界地圖沒有給出具體的投影參數和正反解變換公式[2]，因此有些文獻提出了各自采用多項式擬合確定邊緣經線函數式的方法。

文獻[2]從國家測繪地理信息局發布的等差分緯線多圓錐投影的世界地圖上選擇參考點，采用三次擬合確定邊緣經線的函數式

文獻[3]根據文獻[1]中給出的該投影的設計參數，分析了投影后(x，y)的分布情況，綜合考慮精度和計算量，確定x的擬合曲線應為奇次方程，采用一元五次多項式擬合；y的擬合曲線應為偶次方程，采用一元六次多項式擬合

文獻[4—5]介紹了蘇聯中央測繪科學研究所廣義多圓錐投影方案(1950年)中的邊緣經線的設計函數，其縱坐標x的方程式為緯度的奇次多項式，橫坐標y的方程式為緯度的偶次多項式

雖然文獻[4—5]中沒有指出式(6)可以用于等差分緯線多圓錐投影，但因為等差分緯線多圓錐投影的邊緣經線與該方案的邊緣經線的對稱性質相同，因此可將此擬合函數模型用于等差分緯線多圓錐投影。

分析上述三種邊緣經線擬合函數模型可以發現，在式(4)、式(5)基礎上通過移軸并考慮到邊緣經線的對稱性后，可推導得出如下結果

將投影后的直角坐標原點水平移至邊緣經線的中點，則去除了常數項，式(7)變為

設x=fx(φ)，根據邊緣經線的對稱性，有

fx(-φ)=-fx(φ)-a1φ+a2φ2-a3φ3+a4φ4-…=

-a1φ-a2φ2-a3φ3-a4φ4-…a2φ2+a4φ4+…=

-a2φ2-a4φ4+…

因此有

a2φ2=-a2φ2

a4φ4=-a4φ4

……

為使上述等式成立，a2、a4、a6等偶次項的系數只能為0。

同理，設y=fy(φ)，根據邊緣經線的對稱性，有

fy(-φ)=fy(φ)-b1φ+b2φ2-b3φ3+b4φ4-…=

b1φ+b2φ2+b3φ3+b4φ4+…-b1φ-b3φ3+…=

b1φ+b3φ3+…

因此有

-b1φ=b1φ

-b3φ3=b3φ3

……

為使上述等式成立，b1、b3、b5等奇次項的系數只能為0。這樣式(7)就推導成為式(6)。

實際上，即使在邊緣經線擬合中不考慮其對稱性，采用式(7)這種一般的曲線擬合模型，由于邊緣經線上已知點數據本身的對稱性，使得解出的待定系數中，a2,a4,a6，…和b1,b3,b5，…的系數也接近于0，表現出了其對稱性。但這相當于減少了參與擬合的多項式的項數，因此與使用相同已知點數據和多項式項數但采用式(6)進行邊緣經線擬合的結果相比，式(7)的擬合誤差明顯增大。即使為保證兩式的最高次項的方次相近，將式(7)的項數增至式(6)的兩倍，式(7)的擬合誤差仍大于式(6)。

三、擬合誤差對比

在邊緣經線的多項式擬合中，一般采用最小二乘擬合，如文獻[3]中利用1stOpt軟件進行邊緣經線多項式擬合。也可以采用分段擬合方法，如文獻[6]中選用與投影點的緯度值較近的3～4個已知點進行多項式擬合，這樣做雖然在已知點上的擬合誤差較小，但因為是分段擬合，邊緣經線的連續性稍差。

下面對式(6)和式(7)這兩種不同的擬合函數模型進行最小二乘擬合誤差對比。程序中采用正規方程組法解出最小二乘擬合后的多項式系數a0,a1,…,aM-1,b0,b1,…,bM-1，根據這些系數確定了式(6)、式(7)的邊緣經線多項式，再傳入已知點的緯度值求出擬合后的點，計算它與已知點的距離作為擬合誤差進行對比。

已知點數據：對表2中的14個點，根據邊緣經線對稱性，計算出其緯度為負值的對稱點，再加上表1中的Δλ為180°、φ為0°的點，共29個點作為已知點。

多項式項數和方次：經過測試，式(6)取5項，項數過多時高次項會導致多項式擬合中的震蕩現象，造成投影結果變形。式(7)在兩種條件下與式(6)進行對比。第一種條件下兩式的項數相同，同為5項。第二種條件下兩式的x和y的最高次項方次盡可能相近，測試中x的最高次項的方次可以與式(6)一樣取9，而y的最高次項方次與式(6)相同取10時，擬合結果出現震蕩現象，因此y的最高次項方次取9次，這樣x和y的多項式項數均為10項，為式(6)的二倍，計算時間也相應地增加。

擬合后式(6)的多項式系數為：

a0=1.676 261 862 750 15

a1=-0.000 110 734 568 663 173

a2=1.954 477 061 124 52E-8

a3=-2.668 551 264 636 22E-12

a4=1.344 174 925 330 45E-16

b0=-0.017 050 852 472 941 4

b1=3.941 869 059 706 5E-6

b2=-1.099 928 618 801 22E-9

b3=1.439 009 388 859 4E-13

b4=6.823 147 812 808 01E-18

式(7)取5項時擬合后的多項式系數為：

a0=2.271 571 248 610 32E-14

a1=1.616 105 121 194 76

a2=2.224 643 055 051 45E-17

a3=-4.983 172 137 650 04E-5

a4=2.82518 709 304 736E-21

b0=163.845 946 891 866

b1=-1.469 684 263 135 46E-17

b2=-0.012 379 068 753 132 9

b3=2.030 902 730 958 03E-20

b4=2.952 177 662 935 74E-7

式(7)取10項時擬合后的多項式系數為：

a0=4.745 577 324 584 05E-13

a1=1.673 862 577 481 37

a2=-1.689 756 693 362 57E-15

a3=-0.000 111 052 837 763 137

a4=1.029 503 409 078 15E-18

a5=1.989 551 012 425 34E-8

a6=-2.033 715 263 461 02E-22

a7=-2.650 154 786 566 88E-12

a8=1.251 941 337 157 33E-26

a9=1.269 199 712 877 04E-16

b0=164.593 495 751 958

b1=-9.532 554 139 680 61E-14

b2=-0.014 617 851 141 999

b3=1.618 156 840 527 21E-16

b4=1.418 656 989 887 4E-6

b5=-7.333 870 717 170 98E-20

b6=-1.870 472 805 776 91E-10

b7=1.229 392 976 687 05E-23

b8=9.911 129 417 695 13E-15

b9=-6.829 140 038 805 23E-28

式(7)的兩種擬合結果中，縱坐標x的a0、a2等偶次項系數和橫坐標y的b1、b3等奇次項系數均近于0，可見關于橫軸對稱的已知點數據迫使擬合后的多項式對稱于橫軸。

已知點與對應的擬合后的點之間的距離誤差(本文中稱擬合誤差)見表3。由于式(6)與擬合數據特性相同，也關于橫軸對稱，而式(7)無此特性，因此表3中式(7)的擬合誤差遠大于式(6)，其平均值在多項式項數相同條件下是式(6)的9倍，在最高次項方次相近的條件下是式(6)的4倍。

表3 式(6)與式(7)的擬合誤差對比 cm

四、結束語

在等差分緯線多圓錐投影正解變換中，邊緣經線擬合函數模型的選擇對擬合結果影響很大。文獻[2—3]中用于另一種廣義多圓錐投影的擬合函數模型，考慮到了邊緣經線的對稱性，因為等差分緯線多圓錐投影的邊緣經線具有同樣性質，因此應該選用該擬合函數據模型進行正解變換。對比結果表明，這種擬合函數模型與擬合數據的契合度高，因此擬合誤差較小。

參考文獻：

[1] 胡毓鉅，龔劍文，黃偉. 地圖投影[M]. 北京：測繪出版社，1981：159-163.

[2] 董曼，李勝樂. 世界地圖等差分緯線多圓錐投影的正反解變換[J]. 大地測量與地球動力學，2008，28(2)：95-99.

[3] 葉遠智. 柵格數據的等差分緯線多圓錐投影轉換[J]. 測繪通報，2012(9)：68-70.

[4] 黃國壽. 地圖投影[M]. 北京：測繪出版社，1983：199-201.

[5] 孫達，蒲英霞. 地圖投影[M]. 南京：南京大學出版社，2005：119-121.

[6] 李家，趙晴晴. 等差分緯線多圓錐投影正解變換的參數和方法[J]. 地理空間信息，2012，10(2)：38-40.

[7] 國家測繪地理信息局.《世界地圖》[EB/OL]. [2013-05-20].https:∥map.sbsm.gov.cn/mcp.

[8] 鐘業勛. 不等分經緯線多圓錐投影的設計與解析計算方法[J]. 測繪學報，1965，8(3)：219.

[9] 鐘業勛，童新華，李占元. 橢圓邊經線多圓錐投影[J]. 測繪工程，2012，21(2)：1-5，8.

[10] 張曉盼，齊歡. 地圖投影的最小二乘二元多項式擬合的誤差估計[J]. 測繪科學，2003，28(3)：49-51.