地图投影是建立地球椭球面(或球面)与地图平面之间坐标映射关系的数学方法，作为地图中不可或缺的要素之一，地图投影在描述和研究地理以及地理信息相关问题与现象的过程中至关重要^[1-4]。

按照可展曲面的不同，地图投影可以分为圆柱投影、圆锥投影、方位投影、伪圆柱投影、伪圆锥投影等不同类型^{[2, 5]}。由于地球曲面的不可展性，在地图投影的过程中，会不可避免地发生变形。按照投影变形特征，地图投影可以分为等角投影、等积投影、等距投影等类型^{[2, 6]}。

等角投影和等积投影是两类重要的地图投影方式。等角投影保持地图上各点在各方向的比例相同，即能够保留局部的形状；等积投影则保持地图要素的面积不变。不存在既等角又等积的地图投影^[1]。

地图投影有其数学结构和计算上的复杂性^[7]，研究地图投影的变形是地图制图领域的关键问题之一^[8]。考虑制图区域内的地图变形，将不可避免的变形控制在可以接受的范围之内，可以避免由于使用不恰当的地图投影导致变形量过大的问题，同时也能为地图的选择和设计提供依据^[9-14]。与此同时，地图投影变形也是地球格网化剖分中的重要研究问题^[15-20]。

底索指标(Tissot’s indicatrices)是一种直观的地图投影变形衡量方法。该方法将微分圆^[21]变换到投影平面上，投影之后的微分(正)圆通常会变成长短轴不等的椭圆。底索椭圆能够反映地图投影的距离拉伸、面积变化、形状变化等多个特征^[1]，但底索椭圆只反映椭圆所在位置的变形，而不是地图投影全部位置的变形。此外，底索椭圆只能表示椭圆圆心位置无限小区域的地图投影变形，而不是整个椭圆区域的变形，即底索椭圆可能导致误读等问题的出现^[22]。这些不足导致底索指标在应用上存在一定的局限性。

利用地图投影的微分性质能够研究地图投影的变形^[23]。基于微分的方法依据地图投影正解公式的偏导数，计算得到长度、面积或形状变形系数^[24]。微分指标是一种传统的方法，但该指标计算过程相对复杂(需要进行偏导数的计算)，且要求被评估的地图投影具有连续可微的正解公式。对于直接球面剖分(direct spherical subdivision)方法，其构造是通过递归的方式进行的，不存在任意等级格网的通用构造公式，即不存在连续的地图投影公式。文献[25]给出的结果也表明直接球面剖分的互补比率均值指标并不连续，因此微分指标在使用中也存在应用受限的问题。面积变形系数也可以根据球面微圆的边界参数方程，对球面微圆所在区域进行数值积分求得^[25]，但计算过程也较复杂，且积分运算一般不存在统一的形式。

本文拟提出一种不依赖于微分或积分运算的大圆面积指标，用于计算地图投影的面积变形系数，同时在文献[25]中提出的小圆指标(即互补比率均值指标)基础上简化得到大圆形状指标，用于计算地图投影的形状变形系数。基于大圆面积指标和大圆形状指标，结合随机采样方法，本文将针对不同类型的地图投影，开展评估、比较和分析，验证指标的有效性，并分析地图投影的质量。分析结果表明，大圆指标能够有效降低小圆指标和微分指标的计算量，并能够有效地评估地图投影质量。

1 不依赖微分计算的地图变形指标 1.1 互补比率均值

互补比率均值(averaged ratio between complementary profiles)^[25]ρ_c^s(文献[25]中以ρ_CP表示)是一种评价球面格网形状变形的方法，该指标由两项构成。为了区分小圆指标^[25](带有上标s)与本文中的大圆意义下指标(带有上标g)，本文中以ρ_c^s表示文献[25]中的ρ_CP，其计算公式为

式中，表示过P_i、P和P_j这3点的球面小圆弧长；P_iPP_j表示在小圆弧上过P点的单位切向量；点P为采样区域中心点对应到球面上的点(采样区域中心点的定义见第2节)；点P₁和P₃为采样区域某一条对角线的两端点对应到球面上的点；点P₂和P₄为采样区域另一条对角线的两端点对应到球面上的点。如图 1、图 2所示，图 1(a)中的球面四边形与图 2(c)中的正方形区域相对应，图 1中红色星号表示点P，点P与图 2(c)中心的黑色星号相对应，图 1中洋红色圆形表示大圆弧和的交点。图 1(a)中小圆弧通过点P；图 1(c)为图 1(b)的局部放大图。图 1(b)中的P₁P₃、P₂P₄、P₁PP₃和P₂PP₄并非长度为1的单位向量，只是为了示意各个向量的方向，计算互补比率均值时需要使用单位向量。

小圆指标^[25]需要计算两个球面小圆弧长与和两个切向量(P₁PP₃与P₂PP₄)，其中和P₁PP₃计算过程由式(2)给出^[25]

式中，；p、p₁和p₃分别为点P、P₁和P₃的笛卡儿坐标；·表示向量的点积；×表示向量的叉积；|*|表示向量的模长。计算P₂PP₄和时只需将上述公式中的p₁和p₃替换为p₂和p₄(为了计算上的便利，式(1)和式(2)中采用半径为1的单位球，下文中的式(3)—式(6)、式(9)—式(12)同样使用单位球进行计算)。

文献[25]利用互补比率均值，着重分析和比较了球面离散格网(而非地图投影)的质量。本文基于互补比率均值和球面大圆弧，适当地简化互补比率均值指标，同时提出一种不依赖于微分计算的面积变形指标(见1.2节)，并检验形状变形指标和面积变形指标对评价地图质量的有效性。

1.2 基于大圆弧的形状变形指标与面积变形指标

等角特征和等积特征是评估地图投影变形的两个独立的方面。本节将分别给出用于评估地图投影形状变形和面积变形的指标，关于两个指标的算例分析将在第2节中给出。

首先，本文针对文献[25]提出的小圆形状变形指标(式(1))进行简化，由式(3)给出球面大圆意义下的互补比率均值定义

式中，表示过P_i和P_j两点的球面大圆弧长；P_iP_j表示过P_i和P_j两点的单位切向量。点P₁到P₄的定义与小圆指标中的定义相同，大圆指标无须计算采样区域中心点的坐标，简化了计算过程。

大圆指标需要计算两个球面大圆弧长和和两个向量(P₁P₃和P₂P₄)，其中和P₁P₃计算过程由式(4)给出

对比式(2)与式(4)，两指标中计算两次弧长和两次向量的开销列于表 1中。

小圆形状变形指标的计算过程涉及较多的向量运算(点积、叉积、单位化等)，而大圆形状变形指标运算过程相对简洁。根据表 1的结果，相较于小圆指标，大圆指标能够大大降低计算的复杂度。

本文以球面大圆弧长(或小圆弧长)和向量(或小圆上的切向量)为基础，进一步提出一种衡量面积变形的小圆面积变形指标ρ^s_A和大圆面积变形指标ρ^g_A。ρ^s_A和ρ^g_A的计算依赖于采样区域在球面上的面积A^s和A^g。关于采样区域的定义见第2节，A^s和A^g分别由式(5)和式(6)给出

小圆面积变形指标ρ^s_A和大圆面积变形指标ρ^g_A的计算公式为

式中，S₁表示单个采样区域在投影平面上的面积。由于采样区域为正方形，故S₁=b²，b为采样区域边长。

本节给出的这些指标相比微分指标具有简捷的形式(见第2节算例分析中的比较)。此外，式(3)、式(4)、式(6)和式(8)采用大圆指标的形式，大圆指标相比小圆指标更为简捷。

2 地图投影变形的计算与分析算例 2.1 彭纳投影的投影公式

本节以具体的地图投影(彭纳投影)为例，给出评估地图投影形状变形和面积变形的过程。彭纳投影是一种等积伪圆锥投影，其正解公式^[24]为

式中, σ和ε为中间变量, 其公式为

式中，φ表示纬度；λ表示经度；φ₁为标准纬线；λ₀为中央子午线。

彭纳投影的反解公式^[24]为

式中，中间变量；ρ与φ₁的符号相同。

给定彭纳投影的正反解公式，本文取标准纬线为30°N(φ₁=30°)，中央子午线为本初子午线(λ₀=0°)作为后续计算的条件。

对于等积投影，其底索指标为面积相同的(椭)圆，如图 2(a)所示。在标准纬线和中央子午线上，底索指标显示为大小相同的(正)圆，表明彭纳投影只在这些区域没有(形状和面积)变形。在其他位置，底索椭圆能够反映彭纳投影的形状变形。

2.2 大圆指标的计算

为了使用小圆或大圆指标(见本节)以及微分指标(见2.3节)对彭纳投影的变形进行分析。首先应在地图投影平面上进行等间隔的采样，得到一系列的采样点，如图 2(b)所示，共有624个采样点。为了清晰地展示采样点，图 2中仅使用了较少数量的采样点。

给定采样点(图 2(c)中黑色圆点)后，将相邻的4个采样点所组成的正方形区域记为1个采样区域，(图 2(c)中正方形框)。对于每个采样区域，需要计算中心点(图 2(c)中心位置的黑色星号)的投影坐标。

为了后续计算过程在表述上的方便，记4个采样点或采样区域中心点的投影坐标为(x_p, y_p)(中心点坐标只在计算小圆指标时被使用)，将(x_p, y_p)代入式(11)，得到采样点或中心点的经纬度坐标(φ_p, λ_p)，再根据式(12)将(φ_p, λ_p)转换到单位球面上，得到笛卡尔坐标(X_p, Y_p, Z_p)

将采样点和中心点的笛卡尔坐标代入式(1)(或式(3))可得到形状变形指标，代入式(7)(或式(8))可得到面积变形指标，使用大圆指标(式(3)和式(8))无须计算中心点坐标。

2.3 微分指标的计算

已知彭纳投影的正解公式(式(9))，分别对纬度φ和经度λ求偏导，得到式(13)

将每个投影区域中心点坐标的经纬度值(φ_p, λ_p)代入式(13)，得到偏导数值，再将偏导数值代入式(14)，得到a′和b′

式中

随后，由a′和b′的值可以计算出最大缩放系数a和到最小缩放系数b的数值(式(16))

a与b的比值ρ_d=a/b表征点(φ_p, λ_p)位置的形状变形，a与b的乘积ρ_A=ab表征点(φ_p, λ_p)位置的面积变形。

根据本节的分析，微分方法计算相对复杂，且要求存在连续可微的投影变换公式。

2.4 不依赖微分计算的指标与微分指标的比较

按照不依赖微分计算的指标(ρ_c^g和ρ_c^s)和基于微分方法(ρ_d)得到的形状指标分布如3图所示，图中共计算了10⁶个采样区域。图 3中的结果表明ρ_c^g和ρ_c^s有相近的分布，且在中央子午线和标准纬线附近，彭纳投影的变形程度较低，远离这些区域时，变形程度增加。该结果与按照底索指标给出的彭纳投影特征相符合。

为衡量ρ_c^g与ρ_c^s之间的一致性，本文采用皮尔森积矩相关系数PPMCC(x, y)计算两个变量x与y之间的相关性。PPMCC(x, y)=cov(x, y)/σ_x/σ_y，cov(x, y)表示x和y的协方差，σ_x和σ_y分别为x和y的标准差。PPMCC(x, y)是介于-1与1之间的值，PPMCC(x, y)=1表示x和y完全正相关，PPMCC(x, y)=-1表示x和y完全负相关，PPMCC(x, y)=0表示x和y不相关。

表 2中给出了在不同采样区域数目的条件下ρ_c^s与ρ_c^g的一致性结果。由表 2可知，随着采样点数目的增加，ρ_c^s与ρ_c^g的一致性随之提高，并向1收敛。因此，在采样点足够多的情况下，完全能够用ρ_c^g替代ρ_c^s替代。在下文中，不再计算小圆指标，所有分析和评估的过程均使用大圆指标。

接下来，需要定量衡量指标的ρ_c^g与ρ_d之间的一致性。在10⁶个采样区域的情况下，ρ_c^g和ρ_d的范围分别为[1, 6.424 2]和[1, 12.217 1]，给定ρ_c^g的范围，分析对应的采样区域范围，得到ρ_c^g与ρ_d之间的PPMCC值列于表 3中。

由表 3可知，ρ_c^g与ρ_d的一致性较高，尤其是在ρ_c^g较小的采样区域，两者的一致性更好。因此，小圆或大圆指标能够较好地替代微分指标，衡量地图投影的形状变形。

接下来考虑彭纳投影的面积变形特性。由于彭纳投影是等积投影(没有面积变形)，经过计算，其微分指标ρ_A≡1，这与等积投影的特性相符。进一步计算得到ρ^g_A≈1。由于ρ_A≡1，若使用PPMCC指标分析ρ^g_A与ρ_A的一致性，则会出现cov(ρ^g_A, ρ_A)=0，σ_ρA=0，此时PPMCC中分母为0，无法计算PPMCC的值。

为了能够定量衡量ρ^g_A与1的接近程度，本文使用ρ^g_A的均值μ_A^g，用于衡量ρ^g_A与1的接近程度，使用ρ^g_A的变异系数χ_A^g，用于衡量ρ^g_A的分散(或集中)程度，式(17)给出了χ_A^g的计算公式

式中，N为采样区域数目。

经过计算，对于彭纳投影，本文提出ρ^g_A的均值μ_A^g=1.000 018 6，变异系数χ_A^g=0.000 86，即对于彭纳投影，大圆面积指标与1之间的误差不超过0.02‰，其变异系数值不超过0.9‰。

本节使用大圆指标对彭纳投影进行了定量分析，通过比较，大圆变形指标与微分变形指标结果具有较好的一致性。为了更加全面地验证大圆指标的有效性，下文将利用该指标全面评估不同类型的地图投影。

3 不同类型地图投影质量评估 3.1 等角地图投影的变形

本节将评估3种等角地图投影——墨卡托投影(属于圆柱投影)、兰伯特等角圆锥投影(属于圆锥投影)和方位立体投影(属于方位投影)。由于在接近两极的区域，墨卡托投影存在较为严重的面积变形，且在通常情况下，墨卡托投影并不用来表示两极附近的区域，因此对于墨卡托投影，本文只评估纬度小于85°的区域。兰伯特等角圆锥投影取标准纬线为20°N和60°N，中央子午线为0°，评估范围为纬度介于0°与80°N之间，经度小于80°的区域。方位立体投影取南极点为中心点，评估范围是南半球区域。

对于上述3种等角投影，经过计算ρ_c^s≡1，ρ_c^g≡1，且ρ_d≡1。进一步分析，对于等角地图投影，任一点的ρ_c^s、ρ_c^g和ρ_d均为常数1。因此，使用ρ_c^g能够很好地衡量等角地图投影的等角特性。

图 4中给出了3种等角投影的ρ^g_A分布特性，若ρ^g_A偏离1，表示存在面积变形。由于ρ^g_A和ρ_A具有相近的分布特征，因此图 4中只给出了ρ^g_A的分布。由图 4可知在高纬度区域，墨卡托投影存在较为严重的面积变形，兰伯特等角圆锥投影在靠近80°N的区域面积变形较大，方位立体投影在靠近南极的区域面积变形较大。

进一步计算ρ^g_A与ρ_A的一致性，得到两者的PPMCC值列于表 4中。表中计算结果考虑了约10⁶个采样区域(即使调整采样点间距，也只能使采样区域数目近似为10⁶，但无法保证严格等于10⁶)。

根据表 4中的结果，对于3种等角地图投影，ρ^g_A能够替代ρ_A衡量地图投影的面积变形特征。

3.2 等积地图投影的变形

本节将对10种等积地图投影进行评价，包括3种圆柱投影(兰伯特圆柱投影、贝尔曼圆柱投影和高尔-彼得斯投影)、一种圆锥投影(阿尔伯斯投影)、一种方位投影(兰伯特等积方位投影)、一种改良的方位投影(汉麦尔-埃托夫投影)、两种伪圆柱投影(摩尔维特投影和正弦曲线投影)和两种伪圆锥投影(彭纳投影和Bottomley投影)，其中，彭纳投影已在第2节进行了详细分析。3种圆柱投影的评估范围是纬度小于85°的区域。阿尔伯斯投影取标准纬线为20°N和60°N，中央子午线为0°，评估范围为纬度介于65°S与85°N之间，经度小于80°的区域。兰伯特等积方位投影取南极点为中心点，评估范围是南半球区域。其余5种投影的评估范围为全球。

图 5给出了被评估等积投影的ρ_c^g指标与ρ_d指标的分布情况，进一步计算得到形状变形指标之间的PPMCC值列于表 5中。图 5只对3种圆柱投影中的高尔-彼得斯投影进行了经纬度值的文字标注(经度线之间的间隔为60°，纬度线之间的间隔为40°)，另外两种圆柱投影的经纬度间隔与之相同。

根据图 5和表 5中的结果，ρ_c^g与ρ_d具有相似的分布和较好的一致性(彭纳投影和Bottomley投影的PPMCC值大于0.988 4，阿尔伯斯投影和兰伯特等积方位投影的PPMCC值大于0.991 1，其余地图投影的PPMCC值均大于0.996 4)。

等积投影能够保持地图要素的面积不变，为了验证这一特征，接下来使用不同的面积变形指标对被评估的地图投影进行计算。经过计算，对于上述等积投影的ρ_A≡1，ρ^g_A≈1。为了衡量本文提出的ρ^g_A指标评估等积投影的有效性(与1的接近程度)，分别计算ρ^g_A的均值和变异系数，计算结果列于表 6中。

表 6中的结果表明，ρ^g_A的均值μ_A^g与1之间的偏差小于0.06‰，ρ^g_A的变异系数χ_A^g小于3‰。由此可见，对于被评估的等积投影，ρ^g_A能够较好地替代ρ_A衡量地图投影的面积变形特征。

3.3 任意投影的变形

本节将针对4种任意投影，分析其形状变形和面积变形特征，其中透视方位投影取南极点为中心点，评估范围是南半球区域。米勒投影(属于圆柱投影)和另外两种伪圆柱投影(卡夫拉依斯基Ⅶ投影和Wagner Ⅵ投影)虽然能够表示全球的范围，但这3种投影在高纬度区域，会存在较大的变形，因此对于这3种投影，只评估纬度小于85°的区域。

图 6和图 7给出了被评估任意地图投影的ρ_c^g分布和ρ^g_A分布，图中的卡夫拉依斯基Ⅶ投影和Wagner Ⅵ投影的经纬度网格线具有相同的间隔(经度线之间的间隔为60°，纬度线之间的间隔为40°)。

为验证任意地图投影的形状变形指标和面积变形指标的有效性，表 7给出了大圆指标与微分指标的PPMCC值。根据表 7中的结果，对于被评估的任意地图投影，ρ^g_A能够替代ρ_A衡量地图投影的面积变形(PPMCC值大于1-5×10^-7)，ρ_c^g也能够在一定程度上较好地表征地图投影的形状变形(PPMCC值大于0.988)。

3.4 地图投影变形指标质量评估小结

根据本节对不同地图投影变形指标的计算结果，大圆形状变形指标ρ_c^g以及大圆面积变形指标ρ^g_A与微分指标之间(ρ_d或ρ_A)的一致性结果如下：对于等角投影(见3.1节)，ρ_c^g≡1；对于等积投影(见3.2节)，ρ^g_A的均值和ρ_A的均值之间的偏差小于0.06‰；对于非等角投影(见3.2节和3.3节)，PPMCC(ρ_c^g, ρ_d)>0.988；对于非等积投影(见3.1节和3.3节)，PPMCC(ρ^g_A, ρ_A)>1-5×10^-7。

4 地图投影变形指标的回归分析

第3.1节至3.3节的计算过程是针对每种地图投影单独进行的，本节将综合考虑多个不同的地图投影，进行回归分析，并计算回归误差。对于ρ_c^g和ρ^g_A，线性回归应具有式(18)给出的形式。

针对第3.2节给出的等积投影和第3.3节给出的任意地图投影，采用最小二乘法计算得到的回归系数为k_c^g=0.423 932 79，b_c^g=0.757 870 25。针对第3.1节给出的等角投影和第3.3节给出的任意地图投影，采用最小二乘法计算得到的回归系数为k_A^g=0.995 586 25，b_A^g=0.185 500 11。

按照上述给定的回归系数及其近似值(k_c=0.42，b_c=0.76，k_A=1，b_A=0.18)，由ρ_d和ρ_A值计算得到ρ_c^g和ρ^g_A值与实际的ρ_c^g和ρ^g_A之间的均方根误差列于表 8中。根据表 8中的结果，ρ_c^g平均误差小于1.10‰，ρ^g_A的平均误差小于1.08‰。由此可知，本文提出的指标与现有微分指标具有较好的线性关系，见式(19)

5 地图投影变形指标的随机采样方法

上文给出的变形指标计算过程使用完全采样方法，该过程在投影平面上等间隔地划分采样区域并计算变形指标。本节将在两个坐标维度上，使用随机采样方法，评估地图投影的变形。

蒙特卡罗^[26-27]是一类随机方法的统称，在模拟运算、近似求解、误差分析等问题中有着广泛的应用。本文给出的随机采样方法首先在单位球面上随机生成近似均匀分布的采样点(记采样点数目为t_s)，并使用投影正解公式将球面采样点投影到投影平面，然后在投影平面上以采样点为中心建立正方形的投影区域。记正方形的边长为d_s，则式(7)和式(8)中的S₁=d_s²，并将投影区域的4个顶点按照投影反解公式转换到单位球面上，最后根据式(3)、式(4)、式(6)和式(8)计算变形指标。

由于不同地图投影的投影区域不同，使用完全采样方法需要在投影平面针对每种地图投影进行采样，而随机采样方法只需在球面进行采样，能够避免完全采样方法中采样点不统一和采样过程烦琐等问题。

图 8给出了3.3节中3种任意地图投影ρ_c^g指标和ρ^g_A指标的二维分布情况(取d_s=2×10^-6)。图中椭圆表示对应地图投影大圆变形指标的90%置信区域，3种地图投影的评估范围是纬度小于85°的区域，这与上文内容保持一致，纬度小于85°的区域覆盖了球面上的绝大部分(超过99.6%)，只评估小于85°的区域不会对本节评估的结果和结论带来显著的影响。

对于每种地图投影和给定的采样随机采样点数目t_s(从10³个到10⁵个)，图 8均显示了10次计算过程的90%置信椭圆。图 8(a)由于采样次数(10³)较少，10次计算得到置信椭圆的位置存在一定的偏差，随着采样次数的增加，10次计算得到置信椭圆趋于一致，图 8(c)中的10条置信椭圆重叠度最高。

根据上述分析结果，使用随机采样方法能够使用较少的采样点较好地表征地图投影的失真特征，且结果表明使用约10⁵个随机采样点即能够有效分析地图投影的变形情况。

6 结论

本文提出了一种基于球面大圆弧的地图投影(大圆)面积变形指标，并在互补比率均值的基础上改进得到了地图投影(大圆)形状变形指标。提出和改进的指标与现有的微分指标一致，也与互补比率均值(即小圆指标)一致，但所提出的指标不依赖于微分计算，也不要求地图投影具有连续可微的投影正解公式，且具有较低的计算复杂度。通过回归分析，本文提出的指标与现有微分指标具有较好的线性关系。最后，本文给出了采用随机采样方法进行地图投影质量评估的计算过程，该过程能够有效地降低采样点的数量，使用约10⁵个随机采样点即能够有效分析地图投影的变形情况。