蘭勃特投影向高斯投影轉(zhuǎn)換程序編制及難點分析

胡圣武 鞠全泰 張其華

(河南理工大學(xué) 測繪與國土信息工程學(xué)院 河南 焦作 454000)

0 引言

蘭勃特(Lambert) 投影包括等角圓錐投影和等積方位投影,本文研究的蘭勃特投影是指等角圓錐投影。蘭勃特投影與高斯投影在我國應(yīng)用非常廣泛,我國規(guī)定大于等于1∶50萬地圖均采用高斯投影,據(jù)對高斯地圖投影變形分析和地圖精度的要求,我國規(guī)定比例尺大于和等于1∶100 000地圖采用3°帶投影,比例尺小于和等于1∶25 000地圖采用6°帶投影;而小于等于1∶100地圖,大多采用等角圓錐投影即蘭勃特投影[1-2]。由于二者均是我國地圖投影的主要方式,因此,在實際應(yīng)用中兩者之間坐標(biāo)的相互轉(zhuǎn)換是常見問題[3-5]。

目前國內(nèi)關(guān)于蘭勃特投影和高斯投影的轉(zhuǎn)換問題,主要分為數(shù)值法和解析法兩種方法,其中解析法是地圖投影轉(zhuǎn)換的主流方法,它又分為正解析法和反解析法。

反解析法就是通過中間過渡的方法,反解大地經(jīng)緯度坐標(biāo)(B,L),在帶入到新地圖投影公式中求解其平面坐標(biāo)的方法,即(x,y)→(B,L)→(X,Y)。對于以極坐標(biāo)形式的投影,如圓錐投影,還需要以極坐標(biāo)作為中間過渡。即(x,y)→(ρ,δ)→(B,L)→(X,Y)的方式進(jìn)行轉(zhuǎn)換[6]。這種轉(zhuǎn)化方式不受區(qū)域范圍的限制,在轉(zhuǎn)換過程中,理論上不會存在精度損失。缺點就是計算過程復(fù)雜,計算量大。依賴于使用計算機進(jìn)行數(shù)據(jù)計算[7]。

正解析法就是在復(fù)變函數(shù)理論的基礎(chǔ)上,將某一點的兩種投影的平面坐標(biāo)作級數(shù)展開。通過計算機求解多項式的系數(shù),從而得到兩種等角投影平面坐標(biāo)間的直接轉(zhuǎn)換公式。正解析法的優(yōu)點是不需要任何中間變量,可實現(xiàn)蘭勃特投影向高斯投影平面坐標(biāo)間的轉(zhuǎn)換。不足之處則是轉(zhuǎn)換的精度受轉(zhuǎn)換的方式和展開的階數(shù)有關(guān),也受到轉(zhuǎn)換區(qū)域大小的限制,轉(zhuǎn)換區(qū)域越大精度越低[8-9]。

本文主要利用反解析法,即利用大地經(jīng)緯度作為過渡,先將蘭勃特投影平面坐標(biāo)以極坐標(biāo)的形式表示,再將其轉(zhuǎn)換為大地經(jīng)緯度坐標(biāo),最后再將大地經(jīng)緯度坐標(biāo)代入高斯正算公式轉(zhuǎn)換為高斯平面坐標(biāo)進(jìn)而轉(zhuǎn)換為通用坐標(biāo),并通過編程實現(xiàn)二者的轉(zhuǎn)換,并分析了程序編制中的難點問題。

1 基本理論

1.1 蘭勃特投影

蘭勃特等角圓錐投影是J.H.蘭勃特于1772年所創(chuàng),根據(jù)其與旋轉(zhuǎn)橢球面的交線個數(shù)不同,將其分為蘭勃特切圓錐投影和蘭勃特割圓錐投影[9-10]。圓錐面與旋轉(zhuǎn)橢球的交線成為標(biāo)準(zhǔn)緯線。

蘭勃特投影的坐標(biāo)及變形計算一般公式為

(1)

式中,ρ為緯線投影半徑;ρs為最低緯度的投影半徑;B為緯度;ΔL是地球橢球面上兩條經(jīng)線的夾角;δ是兩條經(jīng)線夾角在平面上的投影;β是小于1的常數(shù),在同一蘭勃特投影中β是不變的,但在不同的蘭勃特投影中β是變化的;m,n為經(jīng)緯線長度比;a,b為極值長度比;ω為最大角度變形值;M是子午圈曲率半徑;(x,y)是平面坐標(biāo)。

1.2 高斯投影

在高斯投影的投影平面上,中央子午線和赤道的投影都是直線。在確定投影帶之后,以任意帶的中央經(jīng)線的投影作為X軸,Y軸則是赤道在高斯投影面內(nèi)的投影,交點O為坐標(biāo)原點,如圖1所示。高斯-克呂格投影的直角坐標(biāo)公式為

圖1 高斯投影示意圖

(2)

式中,s為由赤道至緯度B的經(jīng)線弧長;η=e2cosB;N為卯酉圈半徑。特別要注意的是式中的L是經(jīng)過中心化的經(jīng)度值,不是原經(jīng)度值。

1.3 兩種投影轉(zhuǎn)換實現(xiàn)過程

蘭勃特投影向高斯投影轉(zhuǎn)換流程如圖2所示。

圖2 蘭勃特投影向高斯投影轉(zhuǎn)換圖

2 蘭勃特投影向高斯投影轉(zhuǎn)換的難點問題

2.1 如何從蘭勃特投影的平面坐標(biāo)(x,y)求其大地經(jīng)緯度。

已知(x,y)求解大地坐標(biāo)經(jīng)緯度(B,L),其公式為

(3)

式中,ρ0為中心緯度的投影半徑;l為經(jīng)差;q為等量緯度;Δq為等量緯度差;ΔB為緯度差;b1,b2,b3,b4,b5與緯度有關(guān)的系數(shù)。從式(3)中不能求出緯度B。采用等量緯度反解大地緯度的方法來求解大地緯度B。

要想求解大地緯度B,首先要引入等量緯度q,即

(4)

把子午圈曲率半徑M和卯酉圈曲率半徑N帶入式(4)可得

(5)

經(jīng)過簡單的數(shù)學(xué)變換可得等量緯度與大地緯度的關(guān)系式為

(6)

顯然,當(dāng)B=0°時,q=0。

如果令e1=0,則地球橢球體變?yōu)榍蝮w,橢球面大地緯度B變?yōu)榍蛎婢暥圈?則有

(7)

將上式積分上下限分別為B和B0,經(jīng)過二次計算就可求得B和B0的等量為緯度q和q0。就可以得到由大地緯度計算等量緯度差Δq的公式。為方便實際計算,并不是采用式(7),而是直接采用大地緯度差ΔB來計算等量緯度差Δq的級數(shù)展開式[11-12]為

(8)

經(jīng)過泰勒級數(shù)展開式(8)并設(shè)

(9)

則可得

(10)

逐次求導(dǎo),可得

(11)

則可得

(12)

式中,

(13)

該方法可以簡化ρ的計算,使之更容易被計算機實現(xiàn)?？傻?/p>

(14)

2.2 在高斯正算公式中L的確定

知道了圓錐投影地圖上每點的(B,L),則可按式(2)求出平面坐標(biāo)。不過在代入公式之前要解決L的值的問題。在高斯正算中,每個點有兩種坐標(biāo),即3°度帶和6°度帶。因此,要先求出每個帶的中央經(jīng)度。

3°度帶帶號n與相應(yīng)的中央子午線經(jīng)度L0的關(guān)系為

(15)

6°度帶每帶的中央經(jīng)線度L0和代號n為

(16)

式中,[ ]表示商取整;L為某地點的經(jīng)度。

根據(jù)式(15)和式(16)求出相應(yīng)的中央經(jīng)度L0,則式(2)中L為

(17)

根據(jù)式(2)求出的平面坐標(biāo),x坐標(biāo)一定是正值,y坐標(biāo)有正有負(fù),這樣的坐標(biāo)稱為自然坐標(biāo)值。為了避免坐標(biāo)出現(xiàn)負(fù)值,把y坐標(biāo)進(jìn)行了處理,即加上500 km再加上帶號。這樣的坐標(biāo)稱之為通用坐標(biāo)。

3 利用C#實現(xiàn)蘭勃特投影向高斯投影轉(zhuǎn)換

3.1 蘭勃特投影反算

其流程圖如圖3所示。

圖3 蘭勃特投影反算流程圖

3.2 高斯投影正算

其流程圖如圖4所示。

圖4 高斯投影正算流程圖

3.3 程序界面

程序編寫過程中,考慮到蘭勃特切圓錐還是割圓錐的情況,通過單選框來讓用戶確認(rèn)。程序默認(rèn)了一些投影橢球參數(shù),也給定了自定義橢球參數(shù)的選項,以便用戶在不同參考橢球中都可以使用該程序進(jìn)行坐標(biāo)轉(zhuǎn)換工作,如圖5所示。

圖5 測量程序界面

4 工程實例

在河南省區(qū)圖中,取原點經(jīng)度L0=113°30′,原點緯度B0=35°。由我國各省所選的兩條標(biāo)準(zhǔn)緯線表可知,B1=32°30′,B2=35°30′。其蘭勃特平面坐標(biāo)為(3 021 513.241 9,978 993.339 4)參考橢球為1975國際橢球。通過本程序算得該點在3°帶和6°帶的坐標(biāo),如圖6和圖7所示。

圖6 三度帶成果圖

圖7 六度帶成果圖

5 結(jié)束語

本文對蘭勃特投影向高斯投影坐標(biāo)轉(zhuǎn)換進(jìn)行了系統(tǒng)的理論分析和公式推導(dǎo)。本文應(yīng)用C#語言編寫出視窗體(WinForm)程序,可通過計算機進(jìn)行兩種投影的單點計算。在確定Lambert投影的投影方式(切圓錐還是割圓錐),及投影參數(shù)(原點經(jīng)度和標(biāo)準(zhǔn)緯線緯度),通過計算即可求解得到該點在高斯投影當(dāng)中的坐標(biāo)和通用坐標(biāo)。程序給出了3°帶和6°帶的選項,通過自主選擇,可直接求得不同投影帶之間的坐標(biāo)。

程序編制的過程中,只給出了單點間的兩種投影之間的轉(zhuǎn)換,沒有給出批量點的轉(zhuǎn)換程序。