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

李 家，肖益松

(1. 遼寧省自然地理與空間信息科學(xué)實(shí)驗(yàn)室，遼寧 大連 116027； 2. 遼寧師范大學(xué) 城市與環(huán)境學(xué)院，遼寧 大連 116029)

一、投影的基本公式和設(shè)計(jì)參數(shù)

等差分緯線多圓錐投影已在我國(guó)編制的各種比例尺的世界地圖及其他類(lèi)型世界地圖中得到廣泛應(yīng)用，并獲得了較好的效果[1]。由于在ArcGIS、MapInfo、SuperMap等GIS軟件中沒(méi)有該投影的正反解變換功能，在實(shí)際工作中只能通過(guò)二次開(kāi)發(fā)實(shí)現(xiàn)該功能。

根據(jù)文獻(xiàn)[1]，計(jì)算該投影中央經(jīng)線上x(chóng)0、y0的公式為

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

極坐標(biāo)公式為

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

計(jì)算某一點(diǎn)投影后直角坐標(biāo)的公式為

文獻(xiàn)[1]同時(shí)給出了該投影的設(shè)計(jì)參數(shù)，包括直角坐標(biāo)值(見(jiàn)表1)、緯圈半徑ρ和邊緣經(jīng)線上的極角δφi的值(見(jiàn)表2)。

圖1 等差分緯線多圓錐投影直角坐標(biāo)計(jì)算示意圖

表1 等差分緯線多圓錐投影的直角坐標(biāo)值

續(xù)表

表2等差分緯線多圓錐投影的緯圈半徑ρ和邊緣經(jīng)線上的極角δφ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中的數(shù)據(jù)是基于1∶10 000 000的世界地圖的，(x，y)坐標(biāo)和ρ的單位為厘米。

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

二、邊緣經(jīng)線的擬合函數(shù)模型

國(guó)家測(cè)繪地理信息局在網(wǎng)上正式發(fā)布的世界地圖沒(méi)有給出具體的投影參數(shù)和正反解變換公式[2]，因此有些文獻(xiàn)提出了各自采用多項(xiàng)式擬合確定邊緣經(jīng)線函數(shù)式的方法。

文獻(xiàn)[2]從國(guó)家測(cè)繪地理信息局發(fā)布的等差分緯線多圓錐投影的世界地圖上選擇參考點(diǎn)，采用三次擬合確定邊緣經(jīng)線的函數(shù)式

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

文獻(xiàn)[4—5]介紹了蘇聯(lián)中央測(cè)繪科學(xué)研究所廣義多圓錐投影方案(1950年)中的邊緣經(jīng)線的設(shè)計(jì)函數(shù)，其縱坐標(biāo)x的方程式為緯度的奇次多項(xiàng)式，橫坐標(biāo)y的方程式為緯度的偶次多項(xiàng)式

雖然文獻(xiàn)[4—5]中沒(méi)有指出式(6)可以用于等差分緯線多圓錐投影，但因?yàn)榈炔罘志暰€多圓錐投影的邊緣經(jīng)線與該方案的邊緣經(jīng)線的對(duì)稱(chēng)性質(zhì)相同，因此可將此擬合函數(shù)模型用于等差分緯線多圓錐投影。

分析上述三種邊緣經(jīng)線擬合函數(shù)模型可以發(fā)現(xiàn)，在式(4)、式(5)基礎(chǔ)上通過(guò)移軸并考慮到邊緣經(jīng)線的對(duì)稱(chēng)性后，可推導(dǎo)得出如下結(jié)果

將投影后的直角坐標(biāo)原點(diǎn)水平移至邊緣經(jīng)線的中點(diǎn)，則去除了常數(shù)項(xiàng)，式(7)變?yōu)?/p>

設(shè)x=fx(φ)，根據(jù)邊緣經(jīng)線的對(duì)稱(chēng)性，有

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等偶次項(xiàng)的系數(shù)只能為0。

同理，設(shè)y=fy(φ)，根據(jù)邊緣經(jīng)線的對(duì)稱(chēng)性，有

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等奇次項(xiàng)的系數(shù)只能為0。這樣式(7)就推導(dǎo)成為式(6)。

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

三、擬合誤差對(duì)比

在邊緣經(jīng)線的多項(xiàng)式擬合中，一般采用最小二乘擬合，如文獻(xiàn)[3]中利用1stOpt軟件進(jìn)行邊緣經(jīng)線多項(xiàng)式擬合。也可以采用分段擬合方法，如文獻(xiàn)[6]中選用與投影點(diǎn)的緯度值較近的3～4個(gè)已知點(diǎn)進(jìn)行多項(xiàng)式擬合，這樣做雖然在已知點(diǎn)上的擬合誤差較小，但因?yàn)槭欠侄螖M合，邊緣經(jīng)線的連續(xù)性稍差。

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

已知點(diǎn)數(shù)據(jù)：對(duì)表2中的14個(gè)點(diǎn)，根據(jù)邊緣經(jīng)線對(duì)稱(chēng)性，計(jì)算出其緯度為負(fù)值的對(duì)稱(chēng)點(diǎn)，再加上表1中的Δλ為180°、φ為0°的點(diǎn)，共29個(gè)點(diǎn)作為已知點(diǎn)。

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

擬合后式(6)的多項(xiàng)式系數(shù)為：

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項(xiàng)時(shí)擬合后的多項(xiàng)式系數(shù)為：

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項(xiàng)時(shí)擬合后的多項(xiàng)式系數(shù)為：

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)的兩種擬合結(jié)果中，縱坐標(biāo)x的a0、a2等偶次項(xiàng)系數(shù)和橫坐標(biāo)y的b1、b3等奇次項(xiàng)系數(shù)均近于0，可見(jiàn)關(guān)于橫軸對(duì)稱(chēng)的已知點(diǎn)數(shù)據(jù)迫使擬合后的多項(xiàng)式對(duì)稱(chēng)于橫軸。

已知點(diǎn)與對(duì)應(yīng)的擬合后的點(diǎn)之間的距離誤差(本文中稱(chēng)擬合誤差)見(jiàn)表3。由于式(6)與擬合數(shù)據(jù)特性相同，也關(guān)于橫軸對(duì)稱(chēng)，而式(7)無(wú)此特性，因此表3中式(7)的擬合誤差遠(yuǎn)大于式(6)，其平均值在多項(xiàng)式項(xiàng)數(shù)相同條件下是式(6)的9倍，在最高次項(xiàng)方次相近的條件下是式(6)的4倍。

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

四、結(jié)束語(yǔ)

在等差分緯線多圓錐投影正解變換中，邊緣經(jīng)線擬合函數(shù)模型的選擇對(duì)擬合結(jié)果影響很大。文獻(xiàn)[2—3]中用于另一種廣義多圓錐投影的擬合函數(shù)模型，考慮到了邊緣經(jīng)線的對(duì)稱(chēng)性，因?yàn)榈炔罘志暰€多圓錐投影的邊緣經(jīng)線具有同樣性質(zhì)，因此應(yīng)該選用該擬合函數(shù)據(jù)模型進(jìn)行正解變換。對(duì)比結(jié)果表明，這種擬合函數(shù)模型與擬合數(shù)據(jù)的契合度高，因此擬合誤差較小。

參考文獻(xiàn)：

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

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

[3] 葉遠(yuǎn)智. 柵格數(shù)據(jù)的等差分緯線多圓錐投影轉(zhuǎn)換[J]. 測(cè)繪通報(bào)，2012(9)：68-70.

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

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

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

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

[8] 鐘業(yè)勛. 不等分經(jīng)緯線多圓錐投影的設(shè)計(jì)與解析計(jì)算方法[J]. 測(cè)繪學(xué)報(bào)，1965，8(3)：219.

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

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