【Virtual Globe 渲染技术笔记】2 椭球体基础

椭球体基础

球体在三维空间中由球心( c )和半径( r )定义，所有与球心距离为( r )的点构成了球体表面。为方便起见，球体通常以原点为中心，其隐式方程为：

$$x_s^2+y_s^2+z_s^2=r^2\text{\hspace{1em}(1.1)}$$

满足方程（1.1）的点 $(x_s,y_s,z_s)$ 位于球体表面。这里使用下标 $s$ 表示该点在表面上，以区别于可能在表面内外的任意点$(x,y,z)$。

在某些情况下，将地球建模为球体是合理的，但下一节会提到，椭球体能提供更高的精度和灵活性。以原点 $(0,0,0)$ 为中心的椭球体由沿x、y、z轴的三个半径 $(a,b,c)$ 定义。若点 $(x_s,y_s,z_s)$ 位于椭球体表面，则满足方程（1.2）：

$$\frac{x_s^2}{a^2}+\frac{y_s^2}{b^2}+\frac{z_s^2}{c^2}=1\text{\hspace{1em}(1.2)}$$

当 $a=b=c$ 时，方程（1.2）简化为方程（1.1），因此椭球体即为球体。扁球体是一种椭球体，特别适用于地球建模。扁球体有两个等长的半径（例如 $a=b$ ）和一个较短的第三半径（例如 $c<a$、$c<b$）。等长的较长半径称为半长轴，较短的半径称为半短轴。图1.4展示了具有不同半短轴的扁球体：半短轴与半长轴的差值越大，扁球体就越“扁”。

2.1 WGS84椭球体

在许多应用中，尤其是游戏领域，将地球或其他行星表示为球体是可以接受的。事实上，有些天体几乎呈球形，例如月球——其赤道处的半长轴为1738.1千米，两极处的半短轴为1736千米[180]。而另一些天体则与球体相去甚远，比如火星的卫星火卫一（Phobos），其半径为27×22×18千米。

尽管地球的形状不像火卫一那样奇特，但它也并非完美的球体。地球最适合用扁球体来表示：其赤道半径为6378137米（即半长轴），极地半径为6356752.3142米（即半短轴），这意味着地球赤道处比两极处大约长21384米。

这种地球的椭球体表示被称为WGS84椭球体。截至本书撰写时，它是美国国家地理空间情报局（NGA）最新的地球模型（该模型始于1984年，最后一次更新在2004年）。

WGS84椭球体应用广泛：我们在STK和Insight3D中使用它，许多虚拟地球也采用这一模型。甚至一些游戏也会用到，例如微软的《飞行模拟器》。

在代码中处理地球形状时，最灵活的方法是使用一个可由用户自定义半径构建的通用椭球体类。这使得代码既能支持WGS84椭球体，也能支持其他天体的椭球体模型，如月球、火星等的椭球体。在OpenGlobe中，Ellipsoid类就具备这样的功能（见代码清单1.3）。

/*----------------------------代码清单1.3-----------------------------------*/
class Ellipsoid {
public:
  static const Ellipsoid Wgs84;
  static const Ellipsoid ScaledWgs84;
  static const Ellipsoid UnitSphere;

  Ellipsoid(double x, double y, double z);
  Ellipsoid(const glm::dvec3& radii);

  static glm::dvec3 centricSurfaceNormal(const glm::dvec3& positionOnEllipsoid);
  glm::dvec3 geodeticSurfaceNormal(const glm::dvec3& positionOnEllipsoid) const;
  glm::dvec3 geodeticSurfaceNormal(const Geodetic3D& geodetic) const;

  glm::dvec3 radii() const;
  glm::dvec3 radiiSquared() const;
  glm::dvec3 oneOverRadiiSquared() const;
  double minimumRadius() const;
  double maximumRadius() const;

  std::vector<double> intersections(const glm::dvec3& origin, const glm::dvec3& direction) const;

  glm::dvec3 toDVec3(const Geodetic2D& geodetic) const;
  glm::dvec3 toDVec3(const Geodetic3D& geodetic) const;

  std::vector<Geodetic3D> toGeodetic3D(const std::vector<glm::dvec3>& positions) const;
  std::vector<Geodetic2D> toGeodetic2D(const std::vector<glm::dvec3>& positions) const;

  Geodetic2D toGeodetic2D(const glm::dvec3& positionOnEllipsoid) const;
  Geodetic3D toGeodetic3D(const glm::dvec3& position) const;

  glm::dvec3 scaleToGeodeticSurface(const glm::dvec3& position) const;
  glm::dvec3 scaleToGeocentricSurface(const glm::dvec3& position) const;

  std::vector<glm::dvec3> computeCurve(const glm::dvec3& start, const glm::dvec3& stop, double granularity) const;

private:
  glm::dvec3 _radii;
  glm::dvec3 _radiiSquared;
  glm::dvec3 _radiiToTheFourth;
  glm::dvec3 _oneOverRadiiSquared;
};

2.2 椭球体表面法向量

计算椭球体表面某点的外向表面法向量有诸多用途，包括着色计算以及精确定义地理坐标中的高度。对于球体表面上的点，只需将该点视为向量并进行归一化，即可得到表面法向量。对椭球体表面上的点执行相同操作，得到的是地心表面法向量。称之为“地心”是因为它是从椭球体中心穿过该点的归一化向量。若椭球体不是完美球体，对于大多数点而言，这个向量实际上并非表面的法向量。

另一方面，大地表面法向量才是椭球体表面某点的实际表面法向量。可以想象一个与椭球体在该点相切的平面，大地表面法向量便是垂直于这个平面的向量，如图2.5所示。对于球体，地心表面法向量与大地表面法向量是一致的。而对于更扁的椭球体（如图1.5(b)和图1.5©所示），大多数表面点的地心法向量与大地法向量会出现显著偏差。

大地表面法向量的计算仅比地心法向量稍复杂一些：

椭球方程为隐函数形式：
$$F(x,y,z)=\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}-1=0$$

对于隐函数定义的曲面$F(x,y,z)=0$，其在任一点的法向量方向与函数$F$在该点的梯度（$\nabla F$）方向一致。梯度的定义为：
$$\nabla F=\left(\frac{\partial F}{\partial x},\frac{\partial F}{\partial y},\frac{\partial F}{\partial z}\right)$$

计算椭球的梯度：

• 对$x$的偏导数：$\frac{\partial F}{\partial x}=\frac{2x}{a^2}$
• 对$y$的偏导数：$\frac{\partial F}{\partial y}=\frac{2y}{b^2}$
• 对$z$的偏导数：$\frac{\partial F}{\partial z}=\frac{2z}{c^2}$

因此，椭球面在点$s=(x_s,y_s,z_s)$处的梯度为：
$$\nabla F(s)=\left(\frac{2x_s}{a^2},\frac{2y_s}{b^2},\frac{2z_s}{c^2}\right)$$

梯度方向垂直于椭球面，即梯度是椭球面在该点的法向量（方向向量）。由于法向量的“方向”是核心，常数因子（如这里的2）不影响方向，因此法向量可简化为：

$$\mathbf{m}=\left(\frac{x_s}{a^2},\frac{y_s}{b^2},\frac{z_s}{c^2}\right)$$
$${\hat{\mathbf{n}}}_s=\frac{\mathbf{m}}{\Vert \mathbf{m}\Vert }$$

其中，$(a,b,c)$ 是椭球体的半径，$(x_s,y_s,z_s)$ 是表面上的点，${\hat{\mathbf{n}}}_s$ 是所得的表面法向量。

在实际应用中，$\left(\frac{1}{a^2},\frac{1}{b^2},\frac{1}{c^2}\right)$ 会预先计算并随椭球体一同存储。计算大地表面法向量只需将这个预计算值与表面点进行分量式乘法，然后归一化即可，如代码清单1.4中的 Ellipsoid.GeodeticSurfaceNormal 所示。

代码清单1.5展示了一个非常类似的GLSL函数。传递给 oneOverEllipsoidRadiiSquared 的值通过统一变量（uniform）提供给着色器，因此它在CPU上预先计算一次，即可在GPU上用于多次计算。通常来说，我们会寻找预计算值的方法来提升性能，尤其是在像这样内存开销极小的情况下。

/*------------代码清单1.4：在CPU上计算椭球体大地表面法向量---------------*/
private:
    glm::dvec3 oneOverRadiiSquared_;

// ...

glm::dvec3 Ellipsoid::geodeticSurfaceNormal(const glm::dvec3& positionOnEllipsoid) const {
   return glm::normalize(positionOnEllipsoid * oneOverRadiiSquared_);
}

/*------------代码清单1.5: 在着色器中计算椭球体大地表面法向量---------------*/

vec3 geodeticSurfaceNormal(vec3 p, vec3 oneOverEllipsoidRadiiSquared) {
    return normalize(p * oneOverEllipsoidRadiiSquared);
}

2.3 大地纬度与高度

基于对大地表面法向量的理解，我们可以精确定义地理坐标中的纬度和高度。大地纬度是赤道平面（例如WGS84坐标系中的xy平面）与某点的大地表面法向量之间的夹角。而地心纬度是赤道平面与从原点到该点的向量之间的夹角。在地球的大多数点上，大地纬度与地心纬度并不相同，如图1.6所示。除非另有说明，否则“纬度”一词均指大地纬度。

高度的度量应沿该点的大地表面法向量进行。沿地心法向量度量高度会引入误差，尤其是在较高高度（如太空资产所处的高度）下。地心法向量与大地法向量之间的夹角越大，误差就越大。这种角度差异取决于纬度：在WGS84椭球体上，大地纬度与地心纬度的最大角度差异出现在约45°纬度处。

参考：

• Cozi, Patrick; Ring, Kevin. 3D Engine Design for Virtual Globes. CRC Press, 2011.
• 梯度（Gradient）是微积分中描述多元函数变化率的重要概念，本质上是一个向量，它包含了函数在某一点处沿各个坐标轴方向的偏导数信息。