近年来，随着全球气候变暖，能源形势日益紧张，北极地区的航道价值、科学价值和军事价值日益凸显^[1-2]，研究极区航海导航技术，保证载体在极区的安全航行至关重要，而选择合适的投影方式绘制极区航海图具有重要的意义^[3]。在中低纬地区，航海图通常采用墨卡托投影方法，但随着纬度升高，墨卡托投影长度变形逐渐增大，在极区精度较低^[4]。早期极区航海通常采用极球面投影和横墨卡托投影^[5-7]，在极点附近区域，横墨卡托投影与极区横向坐标导航方法均基于地球球体模型，可实现无缝配合^[8]。由于地球是一个近似的旋转椭球体，随着导航精度要求越来越高，基于椭球体的横向坐标导航方法研究也逐渐成熟^[9]。这就导致球体横墨卡托投影不利于与现代基于地球椭球模型的导航设备（如卫星导航、惯性导航）交互使用，用于极区航行时存在固有原理性误差。因此，需要研究基于地球椭球模型的极区横墨卡托投影。

通用横墨卡托投影或者高斯-克吕格投影是基于椭球体的等角椭圆柱投影^[10]，为了保持精度，通常对该投影进行分带处理，其应用范围也限制在84°N~80°S^[11]，研究已经较为成熟^[12-13]，但其应用并没有拓展到极地地区。并且，极地地区，尤其是近极点区域，面积较小，经线快速收敛，再以经差进行分带不利于完整表达整个极地地形，因此需要研究非分带情形下的极区椭球横向墨卡托投影。

类比于双重投影极球面投影方法^[14]，可尝试用双重投影解决极区椭球横墨卡托投影的问题，但该方法在近极点地区会存在计算奇异和计算溢出问题。通过研究和解决这些数学问题，本文给出完整的横墨卡托投影极区应用情况下的改进方法。

1 双重投影用于极区横墨卡托投影的问题分析

在地图投影中，椭球面在球面上的等角投影可分为完整等角投影和局部等角投影，主要区别反映在经度坐标变换的数学关系上，即：

$ l = \alpha \lambda $

(1)

式中，l为球面经度；λ为椭球面经度；α为比例常数, α=1对应于完整等角投影，α > 1对应于局部等角投影。

由于北极经线快速收敛于极点，北极投影图中经度跨度较大甚至涵盖所有经线，采用局部等角投影可能会出现经度超过180°的情形，不符合经度的定义原则。因此，本文采用完整等角投影方法。

将地球椭球面完整投影到球面上一般要求地球椭球面赤道所在的平面与球面赤道所在的平面共面，且两个赤道同心^[15]。地球椭球面上的子午线投影在球面上为球面子午线，等纬线圈投影到球面上也为等纬线圈。按照这种规定，设地球椭球面投影到球面上的投影变化公式为：

$ \left\{ \begin{array}{l} l = \lambda \\ \varphi = f\left( B \right) \end{array} \right. $

(2)

式中，φ、l为球面纬度、经度；B、λ为椭球面纬度、经度；f(·)为纬度投影坐标变换函数。

由式（2）可求得经线长度比m和纬线长度比n：

$ \left\{ \begin{array}{l} m = \frac{{R{\rm{d}}\varphi }}{{M{\rm{d}}B}}\\ n = \frac{{R{\rm{cos}}\varphi }}{{N{\rm{cos}}B}} \end{array} \right. $

(3)

式中，R为球体半径；M和N分别为子午圈曲率半径和卯酉圈曲率半径。

根据等角投影的条件m=n，利用数学积分的方法可得投影坐标变换公式为：

$ \left\{ \begin{array}{l} l = \lambda \\ \tan \left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{\varphi }{2}} \right) = \tan \left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){\left( {\frac{{1 - e\sin B}}{{1 + e\sin B}}} \right)^{e/2}} \end{array} \right. $

(4)

式中，e为地球椭球体第一偏心率。

设置投影的基准纬度为B₀，对应球面上的纬度为φ₀，利用在该纬线上的长度比为1的条件，可得球体半径R为：

$ R = \frac{{{N_0}\cos {B_0}}}{{\cos {\varphi _0}}} $

(5)

可进一步求得地球椭球面投影到球面上的长度比μ₁为：

$ {\mu _1} = \frac{{R\cos \varphi }}{{N\cos B}} $

(6)

将椭球面等角投影到球面后，再采用球面横墨卡托投影方法把球面投影到平面上，此处直接给出经典的横墨卡托投影公式：

$ \left\{ \begin{array}{l} x = R\arctan \left( {\tan \varphi \sec l} \right)\\ y = \frac{1}{2}R\ln \left( {\frac{{1 + \cos \varphi \sin l}}{{1 - \cos \varphi \sin l}}} \right)\\ {\mu _2} = \frac{1}{{\sqrt {1 - {{\cos }^2}\varphi {{\sin }^2}l} }}\\ {P_2} = \frac{1}{{1 - {{\cos }^2}\varphi {{\sin }^2}l}} \end{array} \right. $

(7)

式中，（x, y）为球面横墨卡托投影后的平面直角坐标；μ₂、P₂分别为球面横墨卡托投影的长度比和面积比。

中低纬度地区通常需要根据经差进行分带。但是极区面积较小，并且经线快速收敛于极点，尤其是近极点地区，如果继续采用分带的方式，会导致极区难以形成完整的表达，因此极区投影图设计将不再进行分带处理。

此外，该投影方法还存在以下两个问题：一是近极点地区的极区投影图通常选北极点为投影中心，即基准纬度为90°，用上述投影方法求取球体半径及接近极点处长度比和面积比时会存在计算奇异的问题，具体表现在式（5）和式（6）中作为分母的余弦函数值为零；二是极点处纬度为90°，接近极点处计算投影坐标时均会计算溢出，并且极点处经度具有多值性，具体表现在式（4）和式（7）中正切函数的计算溢出。

2 横墨卡托投影极区应用改进方法

上述分析表明, 双重投影用于极区横墨卡托投影设计会存在计算溢出和计算奇异问题。为实现完整的极区椭球横墨卡托投影方法，需要有针对性地解决以下4种问题。

1）椭球面到球面等角投影计算溢出问题

由于式（4）纬度坐标变换公式采用正切函数，当纬度B=90°时，出现计算溢出，因此将式（4）中第二式等号左侧变换为正弦函数的形式，即：

$ \begin{array}{l} \sin \left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{\varphi }{2}} \right) = \\ \;\;\;\;\;\frac{1}{{\sqrt {1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{\varphi }{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} }} \end{array} $

(8)

进一步可推导出：

$ \varphi = 2\arcsin \left( {\frac{1}{{\sqrt {1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} }}} \right) - \frac{{\rm{ \mathsf{ π} }}}{2} $

(9)

依据地球椭球面上的子午线投影在球面上为球面子午线，地球椭球面上的纬线圈投影到球面上仍为纬线圈这一规定可知，椭球面上的极点对应于球面上的极点，即B=90°对应φ=90°，经验证满足式（9）。由于式（9）由式（4）变换得到，因此B≠90°时，也必然满足式（9），故式（9）能够解决北极地区投影变换的计算溢出的问题。

2）球体半径、长度比和面积比公式问题

北极极区尤其是近极点地区投影图通常以北极点作为投影中心，即极点处长度变形比为1。由于极点处没有纬线圈，采用纬度长度比n=1为条件求取球体半径不再适用。考虑到极点处任意方向均沿经线方向，因此可以采用经线长度比m=1为条件求取球体半径，具体过程如下：要使

$ m = \frac{{R{\rm{d}}\varphi }}{{M{\rm{d}}B}} = 1 $

(10)

因为：

$ M = \frac{{a\left( {1 - {e^2}} \right)}}{{{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{3/2}}}} $

(11)

式中，a为地球椭球体长半径。可得：

$ \frac{{{\rm{d}}\varphi }}{{{\rm{d}}B}} = \frac{{a\left( {1 - {e^2}} \right)}}{{R{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{3/2}}}} $

(12)

将式（9）看作φ关于B的函数，对该函数求微分可得：

$ \begin{array}{l} \frac{{{\rm{d}}\varphi }}{{{\rm{d}}B}} = \frac{{{{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^{e/2}}}}{{1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}}} \cdot \\ \;\;\;\;\;\;\;\;\;\left[ {{{\csc }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right) - 2{e^2}\frac{{\cos B\cot \left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right)}}{{1 + {e^2}{{\sin }^2}B}}} \right] \end{array} $

(13)

联立式（12）和式（13），并令B=90°，可得：

$ R = \frac{a}{{\sqrt {1 - {e^2}} }}{\left( {\frac{{1 - e}}{{1 + e}}} \right)^{e/2}} $

(14)

式（14）即为椭球面投影到球面上的球体半径。

由于长度比和面积比公式在极点处存在计算奇异，因此推导北极地区长度比公式时，将式（11）、式（13）和式（14）代入式（3）中的经线长度比公式，并利用等角投影经纬线长度比相等这一条件可得：

$ {\mu _1} = m = n = \frac{{R{\rm{d}}\varphi }}{{M{\rm{d}}B}} = \frac{{{{\left( {\frac{{1 - e}}{{1 + e}}} \right)}^{e/2}}{{\left( {\frac{{1 - e\sin B}}{{1 + e\sin B}}} \right)}^{e/2}}{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{3/2}}}}{{{{\left( {1 - {e^2}} \right)}^{3/2}}\left[ {1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} \right]}}\left[ {{{\csc }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right) - 2{e^2}\frac{{\cos B\cot \left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right)}}{{1 - {e^2}{{\sin }^2}B}}} \right] $

(15)

可以看出，式（15）在北极点处没有计算奇异，经验证当B≠90°，式（15）与式（6）等价。

3）球面横墨卡托投影计算溢出问题

由于式（7）中x坐标变换公式采用正切函数、余割函数和反正切函数，当纬度B=90°或l=±90°时，会出现计算溢出，因此，将式（7）右侧变换为反余切函数的形式，即：

$ x = \left\{ {\begin{array}{*{20}{c}} {R\arcsin \frac{1}{{\sqrt {1 + {{\cot }^2}\varphi {{\cos }^2}l} }},}&{\left| l \right| \le {{90}^ \circ }}\\ { - R\arcsin \frac{1}{{\sqrt {1 + {{\cot }^2}\varphi {{\cos }^2}l} }},}&{\left| l \right| > {{90}^ \circ }} \end{array}} \right. $

(16)

上述变换关系的投影原点为赤道，为了使用方便，北极地区投影通常以北极点为投影原点，因此，在式（16）的基础上进行坐标平移，可得：

$ x = \left\{ {\begin{array}{*{20}{c}} {R\arcsin \frac{1}{{\sqrt {1 + {{\cot }^2}\varphi {{\cos }^2}l} }} - \frac{{\rm{ \mathsf{ π} }}}{2}R,}&{\left| l \right| \le {{90}^ \circ }}\\ { - R\arcsin \frac{1}{{\sqrt {1 + {{\cot }^2}\varphi {{\cos }^2}l} }} + \frac{{\rm{ \mathsf{ π} }}}{2}R,}&{\left| l \right| > {{90}^ \circ }} \end{array}} \right. $

(17)

由式（17）可知，B=90°或者l=±90°时，不会出现计算溢出，因此在北极地区不会出现计算溢出。

4）极点处经度多值性问题

由于极点处经线收敛于一点，会出现经度多值性的问题，但由前文公式可知, 在极点处经度取任意值均不会影响计算结果，因此，可以假设极点处经度为零，椭球面投影到球面上的经度变换公式为：

$ \left\{ \begin{array}{l} l = \lambda = 0,\;\;\;\;极点处\\ l = \lambda ,\;\;\;\;\;\;\;\;\;非极点处 \end{array} \right. $

(18)

综上所述，式（9）、式（14）、式（15）和式（18）构成了极区椭球面到球面的投影公式，式（17）和式（7）中后3个公式构成了极区球面横墨卡托投影公式，两者共同构成了基于双重投影的横墨卡托投影极区应用改进方法。

3 子午线收敛角计算

在极区航海导航中，通常需要进行地理真航向和极区航向之间的转换，这势必需要计算子午线收敛角。

由文献[16]可知球面横墨卡托投影的子午线收敛角γ_S为:

$ {\gamma _S} = \arctan \left( {\tan /\sin \varphi } \right) $

(19)

为避免计算奇异，采用等效数学表达式：

$ {\gamma _S} = \arccos \left( {\frac{{\cos l}}{{\sqrt {{{\cos }^2}l + {{\sin }^2}l{{\sin }^2}\varphi } }}} \right) $

(20)

将式（9）代入式（20），可得极区椭球横墨卡托投影子午线收敛角γ_E的计算公式：

$ {\gamma _E} = \arccos \left( {\frac{{\cos l\left[ {1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} \right]}}{{\sqrt {{{\cos }^2}l{{\left[ {1 + {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} \right]}^2} + {{\sin }^2}l{{\left[ {1 - {{\cot }^2}\left( {\frac{{\rm{ \mathsf{ π} }}}{4} + \frac{B}{2}} \right){{\left( {\frac{{1 + e\sin B}}{{1 - e\sin B}}} \right)}^e}} \right]}^2}} }}} \right) $

(21)

4 算例分析

为验证本文改进方法的可行性，采用WGS84椭球模型，长半径a=6 378 137 m，扁率f=1/298.257 223 563，并设置纬度范围为85°N~90°N，经度范围为180°W~180°E，使用mathematica软件可画出该投影的经纬线网图，如图 1所示。从坐标变换公式和仿真过程中可以得出，改进方法计算中不存在计算奇异和计算溢出问题，可以完整表达整个极地地区。

为了比较与球面横墨卡托投影的差异，取球体半径为极点处的平均曲率半径，在上述仿真条件下可得两种方法的平面投影x坐标和y坐标差值dx和dy随经纬度变化的三维分布图，如图 2和图 3所示。

从图 2和图 3可以看出，本文改进方法与球面横墨卡托投影x坐标和y坐标具有较大差别，最大可达5 000 m。由于改进方法采用更加精确的地球椭球模型，因此与球面横墨卡托投影相比精度更高。

在同样仿真条件下，对两种投影方法下子午线收敛角进行计算，计算结果如表 1所示。从表 1可以看出，本文改进方法同一经线上子午线收敛角随纬度升高而增大，近似于经度值。与球面墨卡托投影方法相比，两者子午线收敛角差别较小，误差为秒级，近极点区域甚至小于1″。

本文改进方法与航行阶段导航设备采用的椭球体模型保持一致，投影图上坐标参数与导航设备输出的参数保持一致，投影图和导航设备可以较好地交互使用。而球面横墨卡托投影与导航设备模型不一样，因此必然会引起误差。通过以上分析，极区使用球体横墨卡托投影引起的误差主要表现为平面投影坐标误差较大、子午线收敛角较小的差异。因此，改进方法和导航设备采用一致的地球模型，可以解决由于地球模型不同引起的误差，提高航海绘算精度。

在上述仿真区域内对改进方法的长度变形进行仿真，结果如图 4所示。

从图 4和图 5可以看出，极区椭球面到球面的等角投影长度变形较小，改进方法的长度变形主要是由球面到平面投影的长度变形决定的。由于假设极点为基准纬度，因此长度变形随纬度升高而减小，在纬度88°N到极点的范围内变形较小，如果设置其他纬度为基准纬度，则在基准纬度附近变形较小。因此该投影方法可以配合极区导航技术用于极区航行。

5 结语

本文研究了基于双重投影的横墨卡托投影极区应用改进方法，推导了投影坐标变换公式、球体半径公式、长度比变形公式和子午线收敛角公式。理论分析和对比算例表明，横墨卡托投影极区应用改进方法与球体横墨卡托投影存在明显差异。而改进方法与导航设备采用的椭球体模型相同，不会引起由于地球模型差异导致的误差，可以提高航海绘算精度。改进算法同样适用于极区航海图的程序化设计，可为近极点地区横墨卡托投影和横向坐标导航技术结合使用提供理论基础。