GNSS 误差椭圆：测量平差原理与 PPP 定位精度评估

坐标约定：水平面采用 East–North（E–N），E 为横轴，N 为纵轴，单位为 m；公式中的角度如无说明均为弧度。
术语约定：严格区分经验散布椭圆、总体均值置信椭圆、单历元形式误差椭圆、新观测预测椭圆。四者虽都可以画成椭圆，却不是同一类概率问题。

首先区分：描述误差、估计均值，还是评价某次定位？

同一批 PPP 历元误差可以产生两种中心相同但统计含义不同的椭圆：

• 历元误差的经验散布椭圆：描述所选时间段中，各个历元的水平误差围绕总体/样本平均偏差如何分散。它关注“单个历元误差会怎样变化”，不是某一个特定历元滤波器输出的协方差。
• 平均偏差的均值置信椭圆：描述根据这批历元误差估计未知总体平均偏差，估计有多大不确定性。它关注“我们对平均偏差有多大把握”，不表示 95% 的历元误差点会落在其中。

一句话记忆：散布椭圆看“误差点群有多散”；均值置信椭圆看“平均偏差估得有多确定”；形式误差椭圆看“解算器对某次坐标解的不确定度估计”；预测椭圆看“新观测的误差可能落在哪里”。 这四个问题不能用同一个 95% 概念混答；如果误差均值明显偏离零点，椭圆再小也不代表定位准确。

容易混淆的问题	正确认识
“单历元散布椭圆”是否就是某一历元的形式椭圆？	不是。散布椭圆由多个历元误差的样本协方差得到，是跨历元分布的统计摘要；形式椭圆来自某一历元解算器估计协方差。
“95% 散布椭圆”是什么意思？	在指定的二维正态边缘分布和参数已知的理想模型下，随机抽取一个历元的误差向量落入椭圆的概率为 95%；从同一批样本估计参数后，只能称名义 95%。
“95% 均值置信椭圆”是什么意思？	在统计假设成立时，重复独立抽样并重新构造置信椭圆，约 95% 的随机椭圆会包含固定但未知的总体平均偏差。它不是单历元位置误差的覆盖域。
为什么均值椭圆通常更小？	独立同分布时，均值协方差是总体协方差的 $1/n$；散布协方差则描述单历元的波动。PPP 有时间相关时不能照搬 $1/n$。
散布椭圆与均值椭圆中心相同，是否等价？	不等价。中心相同不代表半轴、覆盖对象、统计假设相同。
95% 经验散布椭圆就是未来单历元的 95% 预测区域吗？	不是严格等价。估计均值、协方差带来的有限样本不确定性，需要预测区域专门处理。

1. 先理解：为什么需要误差椭圆？

误差椭圆（error ellipse）是二维位置误差或估计不确定性的几何表达。 它不是新的定位算法，而是把二维协方差矩阵转成三个直观量：长半轴、短半轴和长轴方向。[1]

例如某测站东向标准差为 2 cm、北向标准差为 1 cm。两个标准差可以说明东西方向更分散，但无法完整说明 E、N 两方向的误差是否相关，也就不能确定真正的最大离散方向。误差椭圆正是为了表达这种方向性和相关性而引入的。

问题	误差椭圆能否回答	原因
水平不确定性主要沿哪个方向伸展？	能	长轴方向表示协方差的最大主方向
哪个正交方向相对更稳定？	能	短轴表示另一主方向
东向与北向误差是否相关？	能	椭圆倾斜反映交叉协方差
结果是否整体偏离真值？	仅靠椭圆尺寸不能	必须同时报告椭圆中心偏差
高程精度怎么样？	不能	二维水平椭圆不含 U 向
故障、粗差或长尾发生的概率？	不能单独判断	二阶统计量无法描述完整分布及时间过程

典型用途包括控制网平差、测站几何精度比较、GNSS 解算形式精度表达、PPP/RTK 不同算法的误差散布对比、形变位移显著性分析。用于形变监测时，应计算位移差的协方差，不能直接以单点椭圆判断位移是否显著。

2. 为什么会得到椭圆，而不是圆？

2.1 二维误差与协方差

定义水平误差向量及其均值：

$$\mathbf{e}=\left[\begin{array}{c}{e}_E\\ e_N\end{array}\right],\mathbf{\mu }=E[\mathbf{e}].$$

对应二维协方差矩阵：

$${\mathbf{\Sigma }}_{EN}=E\left[(\mathbf{e}-\mathbf{\mu })(\mathbf{e}-\mathbf{\mu }{)}^T\right]=\left[\begin{array}{cc}{\sigma }_E^2& {\sigma }_{EN}\\ {\sigma }_{EN}& {\sigma }_N^2\end{array}\right].$$

符号	含义	单位
$e_E,e_N$	东、北方向的定位误差	m
$\mathbf{\mu }$	二维误差期望（平均偏差）	m
${\sigma }_E,{\sigma }_N$	东、北误差标准差	m
${\sigma }_{EN}$	E、N 误差协方差	m²
${\mathbf{\Sigma }}_{EN}$	二维误差协方差矩阵	m²
$E[\cdot ]$、$T$	数学期望、矩阵转置	—

相关系数为：

$${\rho }_{EN}=\frac{{\sigma }_{EN}}{{\sigma }_E{\sigma }_N},-1\le {\rho }_{EN}\le 1,$$

其中相关系数仅在 ${\sigma }_E>0$、${\sigma }_N>0$ 时有定义。

• ${\rho }_{EN}=0$：协方差主轴与 E/N 轴对齐（若两个标准差相等，则是圆）。不相关不必然独立；若联合正态则可推出独立。
• ${\rho }_{EN}\ne 0$：主轴一般倾斜；倾斜方向还取决于 ${\sigma }_E^2,{\sigma }_N^2$。
• ${\sigma }_E={\sigma }_N$ 且 ${\sigma }_{EN}=0$：协方差各向同性，椭圆退化为圆形，长轴方向没有唯一意义。

2.2 特征分解：找出真正的长短轴

对实对称、半正定的协方差矩阵进行特征分解：

$${\mathbf{\Sigma }}_{EN}=\mathbf{V}\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]{\mathbf{V}}^T,{\lambda }_1\ge {\lambda }_2\ge 0.$$

其中 ${\lambda }_1,{\lambda }_2$ 为协方差特征值（m²）；$\mathbf{V}=[{\mathbf{v}}_1,{\mathbf{v}}_2]$ 的两列是单位正交特征向量（无量纲）。

$${\lambda }_{1,2}=\frac{{\sigma }_E^2+{\sigma }_N^2}{2}\pm \frac{\sqrt{({\sigma }_E^2-{\sigma }_N^2{)}^2+4{\sigma }_{EN}^2}}{2}.$$

大特征值对应最大离散方向，小特征值对应正交的最小离散方向。长短轴的几何形状由协方差决定，而不是人为挑选东或北方向。

3. 1σ 椭圆、95% 椭圆分别表示什么？

3.1 椭圆数学公式

假设 $\mathbf{e}$ 服从非退化二维正态分布 $N_2(\mathbf{\mu },{\mathbf{\Sigma }}_{EN})$，则满足：

$$(\mathbf{e}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}_{EN}^{-1}(\mathbf{e}-\mathbf{\mu })\le {\chi }_{2,p}^2$$

的区域具有理论覆盖概率 $p$。这里 ${\chi }_{2,p}^2$ 是自由度为 2 的卡方分布的 $p$ 分位数；左侧为平方马氏距离。由此得到椭圆半轴：

$$\boxed{a_p=\sqrt{{\chi }_{2,p}^2{\lambda }_1}},\boxed{b_p=\sqrt{{\chi }_{2,p}^2{\lambda }_2}}.$$

其中 $a_p,b_p$ 分别为椭圆长、短半轴（m），$p\in (0,1)$ 为理论覆盖概率；计算逆矩阵要求协方差正定（${\lambda }_2>0$）。

二维卡方分布有一个方便的闭式：

$${\chi }_{2,p}^2=-2\ln (1-p),c_p=\sqrt{-2\ln (1-p)}.$$

$c_p$ 是相对于“1σ 主轴椭圆”的半轴放大系数。

名称	半轴放大系数 $c_p$	二维正态理论覆盖率
1σ 主轴椭圆	1.000	39.35%
68.27% 椭圆	约 1.515	68.27%
2σ 主轴椭圆	2.000	86.47%
95% 椭圆	约 2.448	95.00%
99% 椭圆	约 3.035	99.00%

容易混淆的地方：“一维 $\pm 1\sigma$ 覆盖约 68.27%”不能照搬到二维；二维 1σ 主轴椭圆只覆盖约 39.35%。二维 95% 也不是半轴乘 1.96。[1][4]

3.2 椭圆方向角怎样定义？

设长轴单位特征向量 ${\mathbf{v}}_1=[v_E,v_N{]}^T$，在 E 横轴、N 纵轴的坐标系下：

$$\theta =atan2(v_N,v_E),A=(90^{\circ }-{\theta }_{\mathrm{deg}})\mathrm{mod}180^{\circ }.$$

其中 $\theta$ 是从正东逆时针的数学方向角（rad），${\theta }_{\mathrm{deg}}$ 是它换算成度的数值；$A$ 是从正北顺时针的长轴方位角（度）。因为长轴没有箭头，$A$ 在 $[0^{\circ },180^{\circ })$ 内表示即可。

当 ${\lambda }_1\approx {\lambda }_2$ 时椭圆近似圆形，长轴方向容易因微小扰动而大幅改变，此时不宜解读方位角差异。

3.3 “95%”并不自动成立

上述概率是基于误差服从二维正态、均值与协方差已知且正确，并围绕真实误差均值居中给出的理论覆盖率。真实误差均值可以不为零，不能把“中心在均值处”误写成“算法无偏”。实际 PPP 结果可能具有漂移、偏差、时间相关、错固定、粗差或非正态尾部；此时乘以 2.448 只能得到名义 95% 椭圆，不能未经验证就宣称实测覆盖率为 95%。[4]

另外，“某次位置估计的形式置信区域”与“多个历元的误差散布域”外形相似，但协方差来源不同；统计含义见第 6 章。

4. 数值算例：如何计算长短轴？

设一组水平误差近似满足：${\sigma }_E=0.020$ m，${\sigma }_N=0.010$ m，${\rho }_{EN}=0.60$。

于是：

$${\mathbf{\Sigma }}_{EN}=\left[\begin{array}{cc}0.000400& 0.000120\\ 0.000120& 0.000100\end{array}\right]{\mathrm{m}}^2.$$

结果	数值
大特征值 ${\lambda }_1$	0.00044209 m²
小特征值 ${\lambda }_2$	0.00005791 m²
1σ 半长轴 / 半短轴	2.103 cm / 0.761 cm
95% 半长轴 / 半短轴	5.147 cm / 1.863 cm
长轴方位角 $A$	约 70.67°（北起顺时针，模 180°）

虽然东向原始标准差最大，真正的主方向因 E/N 正相关而倾斜。只给 E、N 标准差会漏掉这部分信息。

5. 测量平差中，误差椭圆从哪里来？

线性化最小二乘观测模型：

$$\mathbf{l}=\mathbf{H}\mathbf{x}+\mathbf{\epsilon },E[\mathbf{\epsilon }]=0,\mathrm{Cov}(\mathbf{\epsilon })=\mathbf{R}.$$

当 $\mathbf{H}$ 满列秩、$\mathbf{R}$ 正定且随机模型正确时，广义最小二乘得到：

$$\hat{\mathbf{x}}=({\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H}{)}^{-1}{\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{l},$$

$$\boxed{{\mathbf{P}}_{\hat{x}}=({\mathbf{H}}^T{\mathbf{R}}^{-1}\mathbf{H}{)}^{-1}}.$$

其中 $\mathbf{l}$ 为线性化观测量（m 或其他一致单位），$\mathbf{H}$ 为设计矩阵，$\mathbf{x}$ 为待估参数，$\hat{\mathbf{x}}$ 为估计值，$\mathbf{\epsilon }$ 为观测误差，$\mathbf{R}$ 为观测协方差，${\mathbf{P}}_{\hat{x}}$ 为估计值的形式协方差。设计矩阵可能含有不同单位的参数，参数协方差需按相应量纲理解。

若只有相对权阵 $\mathbf{W}$，并采用 $\mathbf{R}={\sigma }_0^2{\mathbf{W}}^{-1}$，则：

$${\mathbf{P}}_{\hat{x}}={\sigma }_0^2({\mathbf{H}}^T\mathbf{W}\mathbf{H}{)}^{-1}.$$

${\sigma }_0^2$ 是单位权方差；若其未知而由残差估计，必须区分已知方差下的卡方近似与有限自由度下估计方差的置信域。在独立高斯观测及适用的经典平差条件下，后者可涉及 $F$ 分布放大，而不能无条件照搬卡方系数。[1]

工程注意：PPP 常包含钟差、对流层、模糊度等未知参数及状态先验；应读取与最终解状态一致的完整后验协方差矩阵。不能把其他参数删掉后直接对缩减的法方程重新求逆，冒充坐标边缘协方差；固定模糊度后也必须注意正确固定条件下的协方差与错固定风险。

5.1 从 ECEF 坐标协方差变为水平椭圆

如果 PPP 输出完整参数协方差矩阵，先提取其中 ECEF 坐标 $(X,Y,Z)$ 对应的 $3\times 3$ 边缘子矩阵 ${\mathbf{P}}_{XYZ}$。提取这个子矩阵时，钟差、对流层和模糊度的相关性已经通过完整协方差体现，不是先删除这些状态再反求协方差。

设参考点的大地纬度为 $\phi$、大地经度为 $\lambda$（rad），ECEF 到 ENU 的旋转矩阵为：

$$\mathbf{T}=\left[\begin{array}{ccc}-\sin \lambda & \cos \lambda & 0\\ -\sin \phi \cos \lambda & -\sin \phi \sin \lambda & \cos \phi \\ \cos \phi \cos \lambda & \cos \phi \sin \lambda & \sin \phi \end{array}\right].$$

$${\mathbf{P}}_{ENU}=\mathbf{T}{\mathbf{P}}_{XYZ}{\mathbf{T}}^T,{\mathbf{P}}_{EN}={\left({\mathbf{P}}_{ENU}\right)}_{1:2,1:2}.$$

这里 $\mathbf{T}$ 为无量纲的正交旋转矩阵，${\mathbf{P}}_{ENU}$ 为 ENU 三维位置协方差（m²），${\mathbf{P}}_{EN}$ 为其左上角 E/N 二维子矩阵（数学记号采用从 1 开始的下标）。必须使用大地纬度，不能误用地心纬度。[2][3]

对 ${\mathbf{P}}_{EN}$ 进行特征分解，即可画出某一历元定位解的形式误差椭圆。它反映模型给出的不确定性，并不自动包含未建模系统误差、参考坐标误差、错误整数固定或完好性风险。

5.2 两测站的相对误差椭圆

控制网或相对定位还可能关心测站 B 相对于 A 的坐标差。其协方差为：

$${\mathbf{P}}_{B-A}={\mathbf{P}}_B+{\mathbf{P}}_A-{\mathbf{P}}_{BA}-{\mathbf{P}}_{AB},{\mathbf{P}}_{BA}={\mathbf{P}}_{AB}^T.$$

其中 ${\mathbf{P}}_A,{\mathbf{P}}_B$ 为两测站坐标协方差，${\mathbf{P}}_{AB}=\mathrm{Cov}({\hat{\mathbf{r}}}_A,{\hat{\mathbf{r}}}_B)$；各矩阵须在同一坐标系、相同分量顺序下使用。共同误差引起的相关性不能直接忽略。

6. PPP 误差序列可以使用误差椭圆吗？

可以，但首先必须确定自己评估的是哪一种不确定性。

椭圆类型	输入数据	中心	描述对象	是否需要真值
形式误差椭圆	单次解的 ${\mathbf{P}}_{EN,k}$	当前坐标估计值	给定模型下该次坐标估计误差的不确定性	不必
经验散布椭圆	多历元相对真值的 $(e_{E,k},e_{N,k})$	误差样本均值	一般历元误差的横向散布，而非某一指定历元的形式协方差	需要
均值置信椭圆	多历元误差与均值估计的统计分布	样本均值	总体平均偏差（或静态平均位置）的估计不确定性	评估平均偏差时需要
新历元预测椭圆	历史误差样本、预测假设	样本均值	一个尚未观测的新历元误差可能出现的区域	评价定位误差时通常需要

最容易犯的错误是：把经验散布椭圆称为“PPP 平均坐标的 95% 置信椭圆”。 散布、均值置信和新点预测分别回答不同问题，缩放公式也不同。对 PPP 定位误差序列做算法性能评估，通常应优先使用经验散布椭圆 + Bias + RMS + R95 + 收敛指标。

6.1 准备误差：统一真值和坐标框架

静态站在每个历元使用一致的参考坐标；动态站必须有严格时间同步的参考轨迹：

$${\mathbf{e}}_k^{ENU}={\mathbf{T}}_k\left({\hat{\mathbf{r}}}_k^{ECEF}-{\mathbf{r}}_k^{ref,ECEF}\right).$$

其中 $k$ 为历元编号，${\hat{\mathbf{r}}}_k^{ECEF}$ 是 PPP ECEF 解，${\mathbf{r}}_k^{ref,ECEF}$ 为对应历元的参考坐标；${\mathbf{T}}_k$ 为在参考位置计算的 ECEF→ENU 旋转矩阵。静态小范围可以使用固定参考 ENU；大范围动态轨迹宜逐历元定义局部 ENU。如果将不同地理位置的局部 E/N 误差直接汇总到一张椭圆，应明确它表示各时刻“局部东/北”误差统计，而不是同一个固定 ECEF 方向的空间散布。

必须统一坐标框架/参考历元、时间系统、天线参考点及所需的固体潮等位移约定。否则可能把基准不一致造成的偏差错误解释为 PPP 算法误差。

6.2 历元误差的均值与协方差

设水平误差 ${\mathbf{e}}_k=[e_{E,k},e_{N,k}{]}^T$，有效历元数为 $n$：

$$\bar{\mathbf{e}}=\frac{1}{n}\sum _{k=1}^n{\mathbf{e}}_k,\mathbf{S}=\frac{1}{n-1}\sum _{k=1}^n({\mathbf{e}}_k-\bar{\mathbf{e}})({\mathbf{e}}_k-\bar{\mathbf{e}}{)}^T.$$

$\bar{\mathbf{e}}$ 是样本均值/偏差估计（m）；$\mathbf{S}$ 为样本协方差（m²）。在独立同分布的常规假设下，$\mathbf{S}$ 是总体协方差的常用估计；PPP 时间相关或非平稳时，应将其理解为选定时间段内的描述性散布统计，不宜直接依赖独立样本推断性质。

将 $\mathbf{S}$ 代入第 3 章的特征值公式，就能画出名义 $p$ 散布椭圆。这是样本拟合后的描述性结果；同样一批样本中的点落入椭圆的比例不必恰好为 $p$，独立未来历元的真实覆盖概率也需要另行验证。绘图应同时标记：

• 误差原点 $(0,0)$：代表真值。
• 误差均值 $\bar{\mathbf{e}}$：代表平均偏差，也是散布椭圆中心。

若 PPP 长期稳定偏东 5 cm、随机波动仅约 1 cm，椭圆依然可能很小。“稳”不等于“准”。

6.3 散布椭圆与均值置信椭圆：统计对象不同

为避免概念偷换，先用同一批水平误差样本 ${\mathbf{e}}_1,\dots ,{\mathbf{e}}_n$ 定义两种随机对象：

• 单历元误差变量 $\mathbf{e}$：随机选取一个符合所分析工况的历元，其误差具有总体均值 $\mathbf{\mu }$ 和总体协方差 $\mathbf{\Sigma }$。
• 样本均值估计量 $\bar{\mathbf{e}}$：将 $n$ 个历元平均得到的统计量，其估计目标是总体平均偏差 $\mathbf{\mu }$，而非某一历元的误差。

其中 $\mathbf{\mu }$（m）是未知的平均偏差，$\mathbf{\Sigma }$（m²）是单历元误差的总体协方差，$n$ 为样本历元数。两种椭圆不应只用一个“95% 误差椭圆”的名称概括。

A. 经验散布椭圆（dispersion ellipse）

用样本均值 $\bar{\mathbf{e}}$ 和样本协方差 $\mathbf{S}$ 描述单历元误差的散布，其名义 95% 区域为：

$$(\mathbf{e}-\bar{\mathbf{e}}{)}^T{\mathbf{S}}^{-1}(\mathbf{e}-\bar{\mathbf{e}})\le {\chi }_{2,0.95}^2\approx \mathrm{5.991.}$$

式中 $\mathbf{e}$ 是任意一个待描述的二维历元误差（m），${\chi }_{2,0.95}^2$ 是自由度为 2 的卡方 95% 分位数。将估计出的 $\bar{\mathbf{e}},\mathbf{S}$ 代回理论卡方公式只产生近似/名义覆盖域；它不是严格的有限样本预测区域。若观测含明显漂移、非正态或跨工况混合，名义 95% 的概率解释还会失效。

B. 均值置信椭圆（confidence ellipse for the mean）

当 ${\mathbf{e}}_1,\dots ,{\mathbf{e}}_n$ 独立同分布、服从二维正态，且 $\mathbf{\Sigma }$ 未知时，利用 Hotelling $T^2$ 构造总体均值 $\mathbf{\mu }$ 的 95% 联合置信椭圆：

$$(\mathbf{\mu }-\bar{\mathbf{e}}{)}^T{\mathbf{S}}^{-1}(\mathbf{\mu }-\bar{\mathbf{e}})\le \frac{2(n-1)}{n(n-2)}{F}_{2,n-2;0.95},n>2.$$

其中 $F_{2,n-2;0.95}$ 为分子自由度 2、分母自由度 $n-2$ 的 $F$ 分布 95% 分位数。其中心是样本均值，要覆盖的对象是未知的总体均值 $\mathbf{\mu }$，而不是历元误差点。此式只在上述独立二维正态假设成立且 $\mathbf{S}$ 可逆时适用。[5][6]

对二维情况，可将 $F$ 分位数进一步化简，因此无需调用 SciPy 也能计算均值置信椭圆：

$$\boxed{c_{\mathrm{mean},p}^2=\frac{n-1}{n}\left[(1-p{)}^{-2/(n-2)}-1\right]},a_{\mathrm{mean},p}=c_{\mathrm{mean},p}\sqrt{{\lambda }_1},b_{\mathrm{mean},p}=c_{\mathrm{mean},p}\sqrt{{\lambda }_2}.$$

这里 $c_{\mathrm{mean},p}$ 是由二维 Hotelling 置信域推出的半轴系数；${\lambda }_1\ge {\lambda }_2>0$ 为样本协方差 $\mathbf{S}$ 的特征值（m²）；$p$ 为置信水平。对比之下，经验散布椭圆采用 $c_{\mathrm{scatter},p}=\sqrt{-2\ln (1-p)}$，不含 $n$ 引起的均值精度提升因子。

6.4 如果要预测“下一历元”，还需要预测椭圆

经验散布椭圆描述已经观测到的历史点群；均值置信椭圆描述总体平均偏差。如果真正想回答“下一次独立抽取的水平误差，可能落在哪个区域？”，应使用预测椭圆（prediction ellipse）。[6]

若 $n$ 个训练样本及未来新样本 ${\mathbf{e}}_{new}$ 均为独立同分布二维正态，且总体均值、协方差未知，有限样本的预测区域可写为：

$$({\mathbf{e}}_{new}-\bar{\mathbf{e}}{)}^T{\mathbf{S}}^{-1}({\mathbf{e}}_{new}-\bar{\mathbf{e}})\le \frac{2(n+1)(n-1)}{n(n-2)}{F}_{2,n-2;p},n>2.$$

对应预测半轴 $a_{pred,p}=c_{pred,p}\sqrt{{\lambda }_1}$、$b_{pred,p}=c_{pred,p}\sqrt{{\lambda }_2}$，其中

$$c_{pred,p}^2=(n+1)c_{mean,p}^2=\frac{(n+1)(n-1)}{n}\left[(1-p{)}^{-2/(n-2)}-1\right].$$

这里 ${\mathbf{e}}_{new}$ 是尚未观测的一个新历元误差（m），$F_{2,n-2;p}$ 是 $F$ 分布的 $p$ 分位数，$n$ 为训练样本数；${\lambda }_1,{\lambda }_2$ 为样本协方差特征值。这里 $c_{pred,p}$ 是无量纲的预测半轴放大系数。该区域额外考虑了从有限历史样本估计分布参数带来的不确定性；它的覆盖率是针对重复获取训练样本并独立预测下一样本而言，不能理解为给定某一批训练样本后必定有 95% 的条件覆盖率。**它不等于均值置信域，也不等于历史散点外包络。**实际 PPP 高频历元有时间相关、漂移和工况变化，不能未经建模就套用上述 iid 公式。

6.5 用一组数字理解：同一中心、不同半轴

假设有 $n=100$ 个独立同分布的教学样本，统计得到东向与北向样本标准差分别为 2 cm、1 cm，且相关系数为 0；误差均值为东 5 cm、北 0 cm。

指标	名义 95% 散布椭圆	95% 均值置信椭圆	95% 新历元预测椭圆
描述对象	历元误差的散布	总体平均偏差	下一次独立观测误差
中心	$(5,0)$ cm	$(5,0)$ cm	$(5,0)$ cm
半长轴	4.895 cm	0.500 cm	5.021 cm
半短轴	2.448 cm	0.250 cm	2.511 cm
是否覆盖 95% 的现有样本	不作严格保证	不是此含义	不是此含义

三种椭圆都以偏东 5 cm 的样本均值为中心。均值置信椭圆更小，只说明对“平均偏东多少”的估计更确定，并不表示单历元误差更小。 在这个教学例子中，均值置信椭圆与零误差原点不相交，说明在 iid 正态模型下可拒绝“总体平均水平误差为零”的假设；但不能因此忽略 PPP 的时间相关性和参考真值误差。新历元预测椭圆略大于名义散布椭圆，是因为还需考虑参数估计的不确定性。真实 PPP 高频连续历元常相关，不能直接把上述 $n=100$ 的独立样本推断用到实际数据。

6.6 PPP 历元相关时，为什么不能写 $\mathbf{S}/n$？

若不同历元之间存在相关性，则均值协方差的一般形式为：

$$\boxed{\mathrm{Cov}(\bar{\mathbf{e}})=\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n\mathrm{Cov}({\mathbf{e}}_i,{\mathbf{e}}_j)}.$$

其中 $\mathrm{Cov}({\mathbf{e}}_i,{\mathbf{e}}_j)$（m²）是第 $i,j$ 个历元水平误差的交叉协方差。当各历元独立且同分布时，只有 $i=j$ 项不为零，才简化为：

$$\mathrm{Cov}(\bar{\mathbf{e}})=\mathbf{\Sigma }/n.$$

若进一步假设误差序列是二阶平稳的，并定义滞后 $h$ 的矩阵自协方差 $\mathbf{\Gamma }(h)=\mathrm{Cov}({\mathbf{e}}_k,{\mathbf{e}}_{k+h})$（m²），则

$$\mathrm{Cov}(\bar{\mathbf{e}})=\frac{1}{n^2}\left[n\mathbf{\Gamma }(0)+\sum _{h=1}^{n-1}(n-h)(\mathbf{\Gamma }(h)+\mathbf{\Gamma }(h{)}^T)\right].$$

由此可见，时间相关影响的是整个二维协方差矩阵，不能简单给 E/N 统一套一个“有效样本数”就认为得到严格的联合置信椭圆。

PPP 误差往往包含慢变化的轨钟残差、对流层残差、随机过程滤波记忆和收敛过程。相邻误差相关时，不能直接用 $\mathbf{S}/n$ 或上述 Hotelling 系数来给均值椭圆标注严格的 95% 置信水平。对近似平稳、相关随滞后衰减且样本足够长的序列，可估计时间滞后协方差以评估均值不确定性；也可在满足条件时使用移动块自助法等重采样方法。这里即使估计到了均值协方差，也不意味着直接套用原来的 Hotelling 有限样本 $F$ 阈值就有严格的 95% 覆盖率；仍需选用与相关过程相符的统计推断方法。明显非平稳或存在持续系统偏差时，应先分阶段、诊断趋势及偏差来源。

工程上更常见的目的是评价 PPP 定位误差是否准确、稳定：先画散布椭圆，并同时报告 Bias、RMS、R95；只有真的需要对平均偏差或对多历元坐标取平均所得的静态站位置作统计推断时，再谨慎构造均值置信椭圆。多历元坐标的算术平均值与 PPP 滤波器的末历元状态估计并不等价，不能将前者的均值置信椭圆直接当作后者的形式协方差。

6.7 四类椭圆的选择流程（Mermaid）

7. 为什么不能仅用误差椭圆评价 PPP 精度？

误差椭圆主要呈现去均值后的水平散布及相关方向，不直接包含平均偏差，也不描述高程或时间先后。因此应搭配以下指标：

指标	回答的问题	是否包含平均偏差
E/N/U Bias	各方向是否长期偏移？	直接体现
E/N/U STD	各方向围绕均值波动多大？	不包含
E/N/U RMS	相对于真值的总体误差多大？	包含
水平 RMS（RMS_H）	综合水平位置误差多大？	包含
水平 R95	95% 样本的水平径向误差不超过多大？	包含
水平散布椭圆	水平误差分布朝哪里伸展？	需另看椭圆中心
收敛时间/收敛率	多久达到约定精度？	取决于判据
异常率/最大误差	尾部风险如何？	可体现，需定义阈值

定义 $e_{E,k},e_{N,k}$ 为第 $k$ 个历元的水平误差（m），则：

$${\mathrm{RMS}}_E^2=\frac{1}{n}\sum _{k=1}^n{e}_{E,k}^2,{\mathrm{RMS}}_N^2=\frac{1}{n}\sum _{k=1}^n{e}_{N,k}^2,$$

$$\boxed{{\mathrm{RMS}}_H=\sqrt{{\mathrm{RMS}}_E^2+{\mathrm{RMS}}_N^2}},R_{95}=Q_{0.95}\left(\sqrt{e_{E,k}^2+e_{N,k}^2}\right).$$

其中 $Q_{0.95}$ 是样本 95% 分位数算子（具体值受分位数插值规则影响）。R95 是以真值为中心的水平半径分位数，而 95% 椭圆给的是围绕均值的主轴半径；它们不是同一指标。

7.1 RMS、Bias 与样本 STD 的精确关系

若样本标准差 $s_E,s_N$ 采用分母 $n-1$（即 ddof=1），则精确恒等式是：

$${\mathrm{RMS}}_E^2={\bar{e}}_E^2+\frac{n-1}{n}{s}_E^2,{\mathrm{RMS}}_N^2={\bar{e}}_N^2+\frac{n-1}{n}{s}_N^2,$$

$$\boxed{{\mathrm{RMS}}_H^2=\Vert \bar{\mathbf{e}}{\Vert }^2+\frac{n-1}{n}\mathrm{tr}(\mathbf{S})}.$$

其中 ${\bar{e}}_E,{\bar{e}}_N$ 为样本均值，$\Vert \cdot \Vert$ 为欧氏范数，$\mathrm{tr}(\mathbf{S})$ 为样本协方差对角线之和（m²）。

**注意分母：**只有标准差使用分母 $n$ 定义时，才能直接写 ${\mathrm{RMS}}^2={\mathrm{Bias}}^2+{\mathrm{STD}}^2$ 而无需额外因子；当 numpy.std(..., ddof=1) 时应保留 $(n-1)/n$。

7.2 误差椭圆为什么不能替代水平径向准确度？

即使两组 PPP 解的 $\mathbf{S}$ 完全相同，只要均值 $\bar{\mathbf{e}}$ 不同，它们相对真值的 RMS 和 R95 就可能显著不同。误差椭圆的轴长反映的是去均值后的离散特征，中心位置反映平均偏差，两者必须一起查看。在算法论文或测试报告中，不宜把“95% 椭圆半长轴”直接写成“95% 水平定位精度”。

8. 形式误差椭圆怎样与真实 PPP 误差对照？

如果解算器逐历元输出位置形式协方差 ${\mathbf{P}}_{EN,k}$，并有可信的真实误差 ${\mathbf{e}}_k$，可构造归一化平方误差（NEES 在该二维位置子空间的形式）：

$$\boxed{d_k^2={\mathbf{e}}_k^T{\mathbf{P}}_{EN,k}^{-1}{\mathbf{e}}_k}.$$

其中 ${\mathbf{e}}_k$ 必须是估计坐标减真值的误差（不是去样本均值的偏差），${\mathbf{P}}_{EN,k}$ 为正定的该历元估计误差协方差（m²）。在误差针对真值无偏、协方差校准正确且每历元二维高斯模型成立等条件下，$d_k^2$ 的边缘分布为 ${\chi }_2^2$，所以理论 95% 阈值约为 5.991。

可统计 $d_k^2\le 5.991$ 的历元比例，作为形式协方差的经验覆盖率诊断。不要混淆两个问题：

• 历元间可以相关，不一定改变各历元正确的边缘覆盖概率。
• 但相关会减少用于推断覆盖率的有效独立样本数；不能把 1 万个高频历元当成 1 万次相互独立实验。

此外，若参考坐标本身存在不确定度，比较的误差协方差还应考虑参考坐标误差及其与 PPP 解的交叉协方差；否则可能误判形式协方差过于乐观。这里使用的 ${\mathbf{P}}_{EN,k}$ 必须是坐标估计误差的协方差，不是输入观测值协方差，也不是把 E/N/U 三个 STD 简单平方后忽略交叉项得到的矩阵。

对于 PPP-AR，错误整数固定可能造成“小形式椭圆、大实际偏差”；因此椭圆不能替代固定正确性检验或完好性分析。

9. PPP 工程评估的必要条件

1. 先分阶段，再统计：未收敛段、收敛中、稳定段、重收敛段应分别分析。若报告全程结果，需明确含有哪几段。
2. 真值和参考框架一致：包括坐标框架、参考历元、接收机天线 ARP/APC、时间同步及所采用的地球形变模型。
3. 静态与动态分开：静态站用统一真值；动态站须逐时刻对齐参考轨迹。移动轨迹的坐标散布绝不等于定位误差散布。
4. 保留时间信息：单个椭圆无法显示缓慢漂移、突变、周期误差、重新收敛；最好搭配 ENU 误差时间序列。
5. 不随意剔除离群点：粗差处理规则、比例及修正前后结果应透明说明。协方差和椭圆尤其容易受极端值影响。
6. 别误读近圆形椭圆：近圆时长轴方位角数值不稳定，不适合用于算法优劣排序。
7. 样本不能太少或退化：计算二维协方差需至少 3 个非共线样本点；数值上还需保证矩阵正定、条件数可接受。
8. 统一比较口径：不同算法比较时，采用同一真值、时间窗、历元匹配规则、置信水平与剔除策略。

10. 误差椭圆的优点和局限

维度	优点	局限及应对
协方差信息	同时表达大小、主方向、相关性	仅是二阶统计摘要，不代表完整概率分布
测量平差	能展示不同点的精度差异	形式精度依赖随机模型、约束与解算状态
PPP 误差分析	能比较散布形态和方向性	不含平均偏差；同时报告 Bias、RMS、R95
可视化	简洁直观，便于对比	不能表示时间演化、周跳、慢漂移
置信解释	二维正态下有清晰数学公式	非正态/错固定/模型检验选择可破坏名义覆盖率
数值实现	只需协方差和特征分解	需检查秩、病态及单位、轴向定义
测量设计	识别薄弱方向，帮助优化几何	最优设计还需考虑偏差、观测成本与可靠性

11. Python 3.6 示例：绘制 PPP 水平误差椭圆

以下代码既可用模拟数据演示，也可读取真实 PPP 的 ENU 误差（不是原始坐标），输出 Bias、STD、RMS、R95、名义 $p$ 散布椭圆半轴和样本内覆盖率。可选开关 --iid-mean-ci 和 --iid-prediction 分别展示均值置信椭圆与下一独立样本的预测椭圆；这两个开关都仅用于 iid 二维正态教学数据，不应直接用于未经相关性处理的真实 PPP 序列。

输入 CSV 列名：e_m,n_m,u_m，单位均为 m。代码适用于 Python 3.6 语法；请安装与该 Python 版本兼容的 numpy、matplotlib。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""GNSS PPP 水平误差散布椭圆（输入 ENU 误差，单位 m）。"""
import argparse
import csv

import matplotlib.pyplot as plt
import numpy as np

DEFAULT_CONFIDENCE = 0.95
DEFAULT_OUTPUT = "ppp_error_ellipse.png"


def load_errors(filename):
    if not filename:
        # 可重复的测试数据；非真实 PPP 性能数据。
        rng = np.random.RandomState(7)
        cov = [[0.020 ** 2, 0.00012, 0.0],
               [0.00012, 0.010 ** 2, 0.0],
               [0.0, 0.0, 0.030 ** 2]]
        return rng.multivariate_normal([0.050, -0.010, 0.008], cov, 2000)

    with open(filename, "r", encoding="utf-8-sig", newline="") as fp:
        reader = csv.DictReader(fp)
        columns = ("e_m", "n_m", "u_m")
        if not reader.fieldnames or any(c not in reader.fieldnames for c in columns):
            raise ValueError("CSV 需要列名：e_m,n_m,u_m")
        data = [[float(row[c]) for c in columns] for row in reader]
    return np.asarray(data, dtype=float).reshape(-1, 3)


def main():
    parser = argparse.ArgumentParser(description="PPP error ellipse")
    parser.add_argument("--input", default="", help="ENU 误差 CSV；为空时用模拟数据")
    parser.add_argument("--output", default=DEFAULT_OUTPUT)
    parser.add_argument("--confidence", type=float, default=DEFAULT_CONFIDENCE)
    parser.add_argument("--iid-mean-ci", action="store_true",
                        help="仅 iid 二维正态样本：绘制均值置信椭圆")
    parser.add_argument("--iid-prediction", action="store_true",
                        help="仅 iid 二维正态样本：绘制下一独立样本预测椭圆")
    args = parser.parse_args()
    if not 0.0 < args.confidence < 1.0:
        parser.error("--confidence 必须在 0 与 1 之间")

    enu = load_errors(args.input)
    enu = enu[np.isfinite(enu).all(axis=1)]
    n = len(enu)
    if n < 3:
        raise ValueError("至少需要 3 个有限数值的历元")

    en = enu[:, :2]
    bias = enu.mean(axis=0)
    mean_en = bias[:2]
    cov = np.cov(en, rowvar=False, ddof=1)
    values, vectors = np.linalg.eigh(cov)
    order = values.argsort()[::-1]
    values, vectors = values[order], vectors[:, order]
    if values[1] <= 1e-12 * max(values[0], 1e-20):
        raise ValueError("水平误差样本近共线，协方差奇异或严重病态")

    chi2_threshold = -2.0 * np.log(1.0 - args.confidence)
    semiaxes = np.sqrt(chi2_threshold * values)
    major = vectors[:, 0]
    azimuth = (90.0 - np.degrees(np.arctan2(major[1], major[0]))) % 180.0

    t = np.linspace(0.0, 2.0 * np.pi, 361)
    circle = np.vstack((np.cos(t), np.sin(t)))
    boundary = mean_en[:, None] + vectors.dot(semiaxes[:, None] * circle)

    centered = en - mean_en
    rotated = centered.dot(vectors)
    d2 = np.sum(rotated ** 2 / values, axis=1)
    sample_coverage = (d2 <= chi2_threshold).mean()

    std = np.std(enu, axis=0, ddof=1)
    rms = np.sqrt(np.mean(enu ** 2, axis=0))
    horizontal = np.linalg.norm(en, axis=1)
    rms_h = np.sqrt(np.mean(horizontal ** 2))
    r95 = np.percentile(horizontal, 95.0)

    print("Epochs:", n)
    print("Bias ENU (cm):", np.round(bias * 100, 3))
    print("STD  ENU (cm):", np.round(std * 100, 3))
    print("RMS  ENU (cm):", np.round(rms * 100, 3))
    print("Horizontal RMS / R95 (cm): {:.3f} / {:.3f}".format(
        rms_h * 100, r95 * 100))
    print("{:.1f}% ellipse semiaxes (cm):".format(args.confidence * 100),
          np.round(semiaxes * 100, 3))
    if (values[0] - values[1]) / values[0] < 0.01:
        print("Major-axis azimuth: undefined/unstable (nearly circular)")
    else:
        print("Major-axis azimuth from North (deg): {:.2f}".format(azimuth))
    print("In-sample ellipse coverage (not independent validation): {:.1f}%"
          .format(sample_coverage * 100))

    # Hotelling T²：仅 iid 二维正态样本；均值置信域 != 历元散布域。
    mean_boundary, prediction_boundary = None, None
    if args.iid_mean_ci or args.iid_prediction:
        mean_c2 = (n - 1.0) / n * np.expm1(
            -2.0 * np.log1p(-args.confidence) / (n - 2.0))
        if args.iid_mean_ci:
            ci_axes = np.sqrt(mean_c2 * values)
            mean_boundary = mean_en[:, None] + vectors.dot(
                ci_axes[:, None] * circle)
            print("IID mean CI semiaxes (cm):", np.round(ci_axes * 100, 3))
        if args.iid_prediction:
            predict_axes = np.sqrt((n + 1.0) * mean_c2 * values)
            prediction_boundary = mean_en[:, None] + vectors.dot(
                predict_axes[:, None] * circle)
            print("IID prediction semiaxes (cm):", np.round(predict_axes * 100, 3))

    fig, ax = plt.subplots(figsize=(8, 7))
    ax.scatter(en[:, 0] * 100, en[:, 1] * 100, s=7, alpha=0.22,
               label="Epoch errors")
    ax.plot(boundary[0] * 100, boundary[1] * 100, lw=2,
            label="Nominal {:.1f}% dispersion ellipse".format(args.confidence * 100))
    if mean_boundary is not None:
        ax.plot(mean_boundary[0] * 100, mean_boundary[1] * 100,
                lw=2, linestyle="--", label="IID mean confidence ellipse")
    if prediction_boundary is not None:
        ax.plot(prediction_boundary[0] * 100, prediction_boundary[1] * 100,
                lw=2, linestyle=":", label="IID new-epoch prediction ellipse")
    ax.scatter([0], [0], s=130, marker="+", label="Truth (zero error)")
    ax.scatter([mean_en[0] * 100], [mean_en[1] * 100], s=100,
               marker="x", label="Sample mean (bias)")
    ax.set(xlabel="East error (cm)", ylabel="North error (cm)",
           title="PPP horizontal error ellipse")
    ax.set_aspect("equal", adjustable="box")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="best")
    fig.tight_layout()
    fig.savefig(args.output, dpi=200)
    plt.close(fig)
    print("Saved:", args.output)


if __name__ == "__main__":
    main()

用法：

# 演示：模拟数据
python plot_ppp_error_ellipse.py

# 实际数据：输入必须是相对真值的 ENU 误差
python plot_ppp_error_ellipse.py --input ppp_enu_errors.csv --output ppp_95.png

# 调整名义椭圆覆盖概率
python plot_ppp_error_ellipse.py --input ppp_enu_errors.csv --confidence 0.99

# 仅作 iid 二维正态教学演示：均值置信椭圆与新历元预测椭圆
python plot_ppp_error_ellipse.py --iid-mean-ci --iid-prediction

CSV 示例：

e_m,n_m,u_m
0.012,-0.007,0.021
0.018,-0.004,0.024
-0.003,0.006,-0.013

解读提醒：--iid-mean-ci 虚线描述平均偏差的联合置信区域，而 --iid-prediction 点线描述新独立历元的预测区域；两者都不是已有样本的“95% 实测覆盖保证”。未经时间相关性处理的真实 PPP 数据，不应对这些教学曲线赋予严格的 95% 概率解释。In-sample ellipse coverage 是同一组数据拟合椭圆、又拿同一组数据计数的覆盖率，不能当作独立验证的概率，也不代表单次 PPP 解算器形式协方差的覆盖率。

12. PPP 评价流程

13. 与二维误差椭圆相关的其他概念

13.1 水平 R95、CEP、2DRMS 为什么不能直接当成椭圆？

• R95：围绕真值原点、按水平距离排序后的 95% 经验分位数（m）。
• CEP50：水平径向误差达到 50% 覆盖的圆半径，具体定义需说明分布中心及真值约定。
• 2DRMS：术语经常指 $2\sqrt{{\sigma }_E^2+{\sigma }_N^2}$，但应注明标准差和分布假设；它不是一般情况下的严格 95% 圆半径。
• 误差椭圆：由协方差决定的方向性二维区域，通常中心位于误差均值或估计坐标，而不是简单的圆半径。

对于非各向同性、有偏或非正态数据，不存在一个普适系数可以在 R95、CEP、2DRMS 和椭圆半轴之间自由换算。

13.2 三维误差椭球

把 U 向也加入，得到三维协方差 ${\mathbf{\Sigma }}_{ENU}$ 和误差向量 ${\mathbf{e}}_{ENU}$。对于非退化三维正态模型：

$$({\mathbf{e}}_{ENU}-{\mathbf{\mu }}_{ENU}{)}^T{\mathbf{\Sigma }}_{ENU}^{-1}({\mathbf{e}}_{ENU}-{\mathbf{\mu }}_{ENU})\le {\chi }_{3,p}^2.$$

其中 ${\mathbf{\mu }}_{ENU}$ 是三维误差均值，${\chi }_{3,p}^2$ 是自由度 3 的卡方分位数；三条主半轴为 $\sqrt{{\chi }_{3,p}^2{\lambda }_j}$，${\lambda }_j$ 为三维协方差特征值（m²）。三维 95% 的主轴放大系数约为 2.795，不是二维的 2.448。 直接将三维 95% 椭球投影到水平面，其投影区域也不等同于水平边缘分布的二维 95% 椭圆。

14. 记住六个要点

1. 误差椭圆是二维协方差的几何表示，不是新的定位解算方法。
2. 长短轴说明离散大小，方向角说明相关主方向；近圆形时方向角不稳定。
3. 二维 1σ 主轴椭圆正态覆盖率约 39.35%，二维 95% 半轴放大系数约 2.448。
4. PPP 误差序列可以画经验散布椭圆，但必须同时报告偏差、RMS、R95 和收敛表现。
5. 散布椭圆看单历元误差的统计散布；均值置信椭圆看总体平均偏差的不确定性；预测椭圆关注一个新历元。 三者的“95%”不是同一含义。
6. 形式椭圆、经验散布椭圆、均值置信椭圆不能混用；PPP 时间相关时不能直接套用独立样本的均值置信公式。

参考文献

[1] Ghilani, C. D. (2017). Error Ellipse. In Adjustment Computations: Spatial Data Analysis, 6th ed., Chapter 19. John Wiley & Sons. https://doi.org/10.1002/9781119390664.ch19
引用范围：控制网平差的标准误差椭圆、主轴方向、精度可视化。出版社在线章节首次发表为 2017 年，图书版权年份为 2018 年。

[2] Sanz Subirana, J., Juan Zornoza, J. M., & Hernández-Pajares, M. (2011). Positioning Error. ESA Navipedia. https://gssc.esa.int/navipedia/index.php/Positioning_Error
引用范围：GNSS 形式精度、测量精度及水平位置协方差概念；并非 PPP 散布椭圆的专门实验研究。

[3] Sanz Subirana, J., Juan Zornoza, J. M., & Hernández-Pajares, M. (2011). Transformations between ECEF and ENU coordinates. ESA Navipedia. https://gssc.esa.int/navipedia/index.php/Transformations_between_ECEF_and_ENU_coordinates
引用范围：ECEF/ENU 旋转矩阵、使用大地纬度和经度的要求。

[4] Zaminpardaz, S., & Teunissen, P. J. G. (2022). On the computation of confidence regions and error ellipses: a critical appraisal. Journal of Geodesy, 96(2), Article 10. https://doi.org/10.1007/s00190-022-01596-y
引用范围：测量估计与假设检验结合时置信域覆盖率可能偏离名义水平；不是本文 Hotelling 公式的直接来源。

[5] The Pennsylvania State University, Department of Statistics. STAT 505, Lesson 7: Inferences Regarding Multivariate Population Mean. https://online.stat.psu.edu/stat505/Lesson07
引用范围：多元正态样本、Hotelling $T^2$、总体均值联合置信域及 $F$ 分布。

[6] SAS Institute (2007). Parametric Curves: Confidence Ellipses. SAS/INSIGHT User’s Guide. https://support.sas.com/documentation/cdl/en/sgug/59902/HTML/default/fit_sect40.htm
引用范围：二维正态下均值置信椭圆与新独立观测预测椭圆的有限样本公式及区别。

本教程中的 PPP 统计流程、模拟数据、Python 示例由上述理论组织而成，不是对任何文献实验结果的复现。名义概率解释均依赖相应模型假设，实际工程评估应报告验证样本、时间窗和质量控制规则。