UTM投影和Gauss-Krüger投影及其变换实现

UTM投影和Gauss-Krüger投影及其变换实现UTM投影和Gauss-Krüger投影及其变换实现

摘要：本文主要介绍了UTM投影和Gauss-Krüger投影的基本原理和应用，以及它们之间的转换方法和实现过程。其中UTM投影是一种针对局部区域的投影系统，适用于纬度范围在84°S至84°N之间的地区；Gauss-Krüger投影则是一种基于椭球体模型的投影系统，适用于整个地球表面。两种投影系统都可以通过一定的数学公式进行坐标变换，实现它们之间的转换。文章还介绍了常见的UTM和Gauss-Krüger坐标转换算法，以及它们在实际应用中的具体实现方法。

关键词：UTM投影；Gauss-Krüger投影；坐标转换；实现过程

UTM投影和Gauss-Krüger投影简介

UTM（UniversalTransverseMercator）投影最早由美国陆军和英国王家地图局于1947年联合研制而成，是一种针对局部区域的投影系统。该系统将地球表面划分成60个跨度为6度的投影带，每个带都有一个中央经线，将带内的地图平面投影到以中央经线为纵轴、赤道为横轴的平面上。UTM投影适用于纬度范围在84°S至84°N之间的地区，因此它在中低纬度地区的应用比较广泛。

Gauss-Krüger投影是一种基于椭球体模型的投影系统，由德国天文学家高斯和克吕格尔于19世纪中叶研制而成。该系统将地球表面划分成一系列带状区域，每个带在横向上跨越3度带，并沿纵向采用等角度投影。它适用于整个地球表面，并且在高纬度地区的精度更高，因此在极地和高纬度地区的应用比较广泛。

UTM投影和Gauss-Krüger投影的转换

由于UTM投影和Gauss-Krüger投影系数不同，因此它们之间的坐标需要进行转换。下面介绍两种常见的转换方法：

方法一：由UTM转换为Gauss-Krüger

UTM协调制下的坐标(x,y)可以通过以下公式转换为Gauss-Krüger下的坐标(N,E)：

其中，Δλ为中央经线λ0与自定义经度λ之差；e2为椭球体第一偏心率的平方，f为椭球体扁率。

方法二：由Gauss-Krüger转换为UTM

Gauss-Krüger协调制下的坐标(N,E)可以通过以下公式转换为UTM下的坐标(x,y)：

其中，k0为中央子午线比例因子；α为中央子午线上一点的方位角；μ为坐标点与赤道的夹角余切值。

实现过程

UTM投影和Gauss-Krüger投影的转换可以通过程序实现。下面介绍一种基于Matlab的实现方法。

方法一：由UTM转换为Gauss-Krüger

（1）定义输入变量：UTM投影下的坐标值x、y；中央经线λ0；椭球体参数a、b。

（2）计算各项参数：

k0=0.9996;

f=(a-b)/a;

e2=2*f-f^2;

n=(a-b)/(a+b);

A=a/(1+n)*(1+n^2/4+n^4/64);

α=[(x-500000)/(A*k0)+λ0]*pi/180;

t=tan(α);

η2=e2*cos(α)^2;

M=A*[((1-e2/4-3*e2^2/64-5*e2^3/256)*α-(3*e2/8+3*e2^2/32+45*e2^3/1024)*sin(2*α)+(15*e2^2/256+45*e2^3/1024)*sin(4*α)-35*e2^3/3072*sin(6*α))];

ρ=A*(1-e2)/(1-e2*sin(α)^2)^(3/2);

ν=A/sqrt(1-e2*sin(α)^2);

（3）计算高斯投影坐标N、E：

N=M+t/2*ρ*ν^2*(1-ρ/12*(5+3*t^2+η2-9*t^2*η2));

E=k0*ν*[t*ρ*(1+η2)/6+ρ^3*(1-t^2+η2+15*t^2*η2)/120];

方法二：由Gauss-Krüger转换为UTM

（1）定义输入变量：Gauss-Krüger投影下的坐标值N、E；中央子午线λ0；椭球体参数a、b；带号zone。

（2）计算各项参数：

f=(a-b)/a;

e2=2*f-f^2;

n=(a-b)/(a+b);

A=a/(1+n)*(1+n^2/4+n^4/64);

ω=(E-500000)/(k0*A);

sinω=sin(ω);

cosω=cos(ω);

cosω2=cosω^2;

cosω4=cosω2^2;

cosω6=cosω4^2;

sinω2=sinω^2;

sinω4=sinω2^2;

sinω6=sinω4^2;

τ=tan(ω^2);

τ2=τ^2;

τ4=τ2^2;

τ6=τ4^2;

η2=e2*cosω2;

η4=η2^2;

η6=η4*η2;

Nf=(a-b)/(a+b);

α=[1+Nf^2/4+Nf^4/64]*A;

β=[3*Nf/2-27*Nf^3/32]*sin(2*ω);

γ=[21*Nf^2/16-55*Nf^4/32]*sin(4*ω);

δ=[151*Nf^3/96]*sin(6*ω);

ε=[1097*Nf^4/512]*sin(8*ω);

N=α*(ω+β*sinω2+γ*sinω4+δ*sinω6+ε*sinω4);

（3）计算UTM坐标x、y：

y=N-[k0*A*(1-η2+e2*cosω2)/(1-e2*sinω2)^(3/2)]*τ2/2+k0*A*(1-η2+e2*cosω2)/(1-e2*sinω2)^(5/2)*τ4/24*[5+3*τ2+10*η2-4*η4-9*τ2*η2];

x=[k0*A/(1-e2*sinω2)^(1/2)]*[1+τ2/(2*(1-e2*sinω2))*[-1-η2+3*cosω2+η4-9*τ2*η2]+τ4/(24*(1-e2*sinω2))^2*[1+15*τ2-15*η2-15*τ2*η2+5*η4+1*η6-90*τ2*η4]+τ6/(720*(1-e2*sinω2))^3*[-1-28*τ2+24*τ4+3*τ6+16*η2+192*τ2*η2-142*η4-131*τ2*η4]];

x=x+500000+(zone-1)*1000000;

结论

本文介绍了UTM投影和Gauss-Krüger投影的原理和应用，以及它们之间坐标变换的基本方法和实现过程。我们发现，两种投影系统之间的转换方法虽然有些复杂，但通过计算机程序可以实现自动化，大大提高了工作效率。建立合理的坐标变换模型是地理信息系统应用的基础，对于跨越不同坐标系的数据整合和地图制作具有重要意义。UTM投影和Gauss-Krüger投影是地图学和测量学领域中，常用的两种投影方法。两种投影方法都在地图制作、地理测量、导航、军事等领域得到广泛应用。但是两种投影方法是不同的，而且分别适用于不同的区域，因此需要进行坐标转换。

UTM投影

UTM投影是一种局部投影系统，主要用于小范围地图制作。它将地球表面分成60个以6度为跨度、由中央子午线为分界线的狭长条带。每个带内部采用横向等角投影，即沿着子午线方向采用墨卡托投影，沿着赤道方向采用极圆投影。UTM投影系统划分的60个带覆盖的纬度范围是-80度至+84度。因此，UTM投影适用于纬度范围在84度南至84度北之间的地区。

UTM的坐标系统采用以120度为中央子午线的经度坐标和相对于赤道的平面坐标。坐标单位为米。UTM的坐标系统又分为北半球（N）和南半球(S)两个区域。

UTM的坐标转换公式：

UTM坐标可以通过以下公式转换为经纬度坐标：

其中，k0为中央子午线比例因子；E和N为UTM投影下的坐标值；λ0为中央子午线经度。

Gauss-Krüger投影

Gauss-Krüger投影是一种全球性的投影方法。它是基于高斯椭球体模型构建的，可同时适用于整个地球表面的任何区域。Gauss-Krüger投影的特别之处在于，它消除了地球椭球体的影响，使得地图上的距离尽可能地接近实际的地面距离。因此，在高纬度地区的精度要比UTM投影高，特别适合用于极地和高纬度地区。

Gauss-Krüger投影的坐标系统是基于高斯-克吕格尔坐标系统，采用横向等角投影，即沿着中央子午线方向采用高斯圆柱投影。Gauss-Krüger投影按照3度分带，从东经0度开始，以3度间隔往西分带，以每个带的中央经线与本初子午线之间的纬度带的宽度为3度。每个带内使用相同的投影系数和椭球体参数。

Gauss-Krüger的坐标转换公式：

Gauss-Krüger坐标可以通过以下公式转换为经纬度坐标：

其中，N和E为高斯-Krüger坐标，λ0为中央子午线经度，a、f为WGS84参考椭球体参数。k0为中央子午线的投影尺度因子。

坐标转换的实现过程

坐标转换的实现可以通过编写程序实现。常用的编程语言有Matlab、Python和C++等。由于Matlab编程简单，读者可以借助Matlab工具箱和函数库来实现。

Matlab中实现UTM和Gauss-Krüger的坐标转换需要掌握Matlab的基础语法和公式计算。下面以Matlab为例，介绍UTM和Gauss-Krüger的坐标转换。

一、由UTM转换为经纬度

1.输入变量：

x,y：UTM投影下的坐标值

λ0：中央子午线经度

a：参考椭球体的长半轴

e：参考椭球体的离心率

2.计算各项参数：

f=1/298.257223563;

b=a*(1-f);

e_=sqrt(e^2/(1-e^2));

n=(a^2-b^2)/(a^2);

p0=a*(1-n)*((1+n^2/4+n^4/64)/(1+n))^(-1);

k0=0.9996;

x_=x-k0*500000;

y_=y;

M=y_./k0;

mu=M./(p0*(1+(1+n+5/4*n^2+5/4*n^3)*(pi/180)*mu.^2+...

3*(1+3/2*n+45/16*n^2+175/64*n^3)*(pi/180)*mu.^4/4+...

(15/8*n^2+15/8*n^3)*(pi/180)*mu.^6/6+(35/24*n^3)*(pi/180)*mu.^8/8));

e1=(1-sqrt(1-e^2))./(1+sqrt(1-e^2));

phi1=mu+(3*e1/2-27/32*e1^3)*sin(2*mu)+(21/16*e1^2-55/32*e1^4)*sin(4*mu)+...

(151/96*e1^3)*sin(6*mu)+(1097/512*e1^4)*sin(8*mu);

C1=e^2*cos(phi1).^2./(1-e^2);

T1=tan(phi1).^2;

N1=a./sqrt(1-e^2.*sin(phi1).^2);

R1=N1.*(1-e^2)./(1-e^2.*sin(phi1).^2);

D=x_./(N1.*k0);

phi=phi1-(N1.*tan(phi1)./R1).*(D.^2/2-(5+3*T1+10*C1-4*C1.^2-9*e^2).*(D.^4)/24+(61+90*T1+298*C1+45*T1.^2-252*e^2-3*C1.^2).*(D.^6)/720);

lambda=λ