如何利用三点生成平滑的卫星轨迹大圆弧？

我正尝试编写程序，利用卫星过地平线时的起始方位角与仰角、最大方位角与仰角、终止方位角与仰角这三个点，生成卫星轨迹的俯视极坐标图。目标效果类似目视过境示意图，但为俯视视角。

目前我使用球面线性插值（SLERP）分别生成起点到最高点、最高点到终点的两段大圆弧，再拼接成单段弧线，但生成的弧线在接近最高点处存在凸起，并非平滑曲线。我猜测是两段SLERP弧线拼接导致了该问题，请问是否有更合适的公式，可利用这三个点生成单段平滑的大圆弧？

附测试代码：

import numpy as np
import matplotlib.pyplot as plt

def spherical_to_cartesian(az_deg, el_deg, radius=1.0):
    """将球面坐标（方位角、仰角、半径）转换为笛卡尔坐标。"""
    az = np.radians(az_deg)
    el = np.radians(el_deg)
    
    x = radius * np.cos(el) * np.sin(az)
    y = radius * np.cos(el) * np.cos(az)
    z = radius * np.sin(el)
    return np.array([x, y, z])

def slerp(v0, v1, t):
    """在两个向量间执行球面线性插值（SLERP）。"""
    v0 = v0 / np.linalg.norm(v0)
    v1 = v1 / np.linalg.norm(v1)
    dot = np.clip(np.dot(v0, v1), -1.0, 1.0)
    omega = np.arccos(dot)
    
    if np.isclose(omega, 0):
        return v0
    return (np.sin((1 - t) * omega) * v0 + np.sin(t * omega) * v1) / np.sin(omega)

def interpolate_hemisphere_arc(start_az, start_el, mid_az, mid_el, end_az, end_el, samples=50, radius=1.0):
    """在北半球上生成三点（方位角/仰角）间的插值弧线。"""
    start = spherical_to_cartesian(start_az, start_el, radius)
    mid   = spherical_to_cartesian(mid_az, mid_el, radius)
    end   = spherical_to_cartesian(end_az, end_el, radius)

    path_xy = []

    # 插值起点→最高点
    for t in np.linspace(0, 1, samples // 2 + 1):
        p = slerp(start, mid, t)
        path_xy.append((p[0], p[1]))

    # 插值最高点→终点（跳过重复的最高点）
    for t in np.linspace(0, 1, samples // 2 + 1)[1:]:
        p = slerp(mid, end, t)
        path_xy.append((p[0], p[1]))

    return path_xy

if __name__ == "__main__":
    # 测试用例：绘制单段插值卫星轨迹
    path = interpolate_hemisphere_arc(
        start_az=169.67,
        start_el=0,
        mid_az=256.75,
        mid_el=75.92,
        end_az=345.34,
        end_el=0
    )

    xs, ys = zip(*path)
    xs = np.array(xs)
    ys = np.array(ys)

    r = 90 - 90 * np.sqrt(xs**2 + ys**2)
    theta = np.arctan2(xs, ys)

    plt.figure(figsize=(6, 6))
    ax = plt.subplot(111, polar=True)
    ax.plot(theta, r)

    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)
    ax.set_title("卫星插值轨迹的俯视极坐标图")
    ax.set_rlim(90.0, 0.0)

    plt.tight_layout()
    plt.show()

问题核心是两段SLERP拼接的曲线分属两个不同大圆，卫星实际过境轨迹是单一大圆的一部分，因此拼接处会出现曲率突变，导致凸起。正确做法是先通过三个点确定唯一大圆，再在该大圆上生成平滑轨迹。

实现思路 #

1. 计算大圆法向量：利用三个点的笛卡尔坐标，通过两次叉乘得到大圆的单位法向量。
2. 罗德里格斯旋转生成轨迹：以起点为基准，绕法向量旋转不同角度，生成大圆上的连续点，确保轨迹经过起点、最高点和终点。
3. 坐标转换适配绘图：将笛卡尔坐标投影到xy平面，转换为俯视极坐标所需的参数。

修改后的代码 #

import numpy as np
import matplotlib.pyplot as plt

def spherical_to_cartesian(az_deg, el_deg, radius=1.0):
    az = np.radians(az_deg)
    el = np.radians(el_deg)
    x = radius * np.cos(el) * np.sin(az)
    y = radius * np.cos(el) * np.cos(az)
    z = radius * np.sin(el)
    return np.array([x, y, z])

def great_circle_angle(a, b):
    """计算两个单位向量间的大圆角度"""
    dot = np.dot(a, b)
    return np.arccos(np.clip(dot, -1.0, 1.0))

def great_circle_interpolation(start, mid, end, samples=100):
    # 计算大圆法向量
    vec1 = mid - start
    vec2 = end - start
    normal = np.cross(vec1, vec2)
    normal = normal / np.linalg.norm(normal)

    # 计算起点到终点的总大圆角度
    total_angle = great_circle_angle(start, end)
    # 生成均匀分布的角度参数
    t_values = np.linspace(0, total_angle, samples)

    path = []
    for t in t_values:
        # 罗德里格斯旋转公式：绕法向量旋转t角度
        cos_t = np.cos(t)
        sin_t = np.sin(t)
        rotated = start * cos_t + np.cross(normal, start) * sin_t + normal * np.dot(normal, start) * (1 - cos_t)
        rotated = rotated / np.linalg.norm(rotated)
        path.append(rotated)
    
    return np.array(path)

def interpolate_satellite_track(start_az, start_el, mid_az, mid_el, end_az, end_el, samples=100):
    start = spherical_to_cartesian(start_az, start_el)
    mid = spherical_to_cartesian(mid_az, mid_el)
    end = spherical_to_cartesian(end_az, end_el)
    
    # 生成单段大圆轨迹
    great_circle_points = great_circle_interpolation(start, mid, end, samples)
    # 提取xy平面投影坐标
    path_xy = [(p[0], p[1]) for p in great_circle_points]
    return path_xy

if __name__ == "__main__":
    path = interpolate_satellite_track(
        start_az=169.67,
        start_el=0,
        mid_az=256.75,
        mid_el=75.92,
        end_az=345.34,
        end_el=0
    )

    xs, ys = zip(*path)
    xs = np.array(xs)
    ys = np.array(ys)

    r = 90 - 90 * np.sqrt(xs**2 + ys**2)
    theta = np.arctan2(xs, ys)

    plt.figure(figsize=(6, 6))
    ax = plt.subplot(111, polar=True)
    ax.plot(theta, r)

    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)
    ax.set_title("卫星插值轨迹的俯视极坐标图")
    ax.set_rlim(90.0, 0.0)

    plt.tight_layout()
    plt.show()

代码说明 #

• great_circle_interpolation函数：通过罗德里格斯旋转公式生成单一大圆上的连续点，确保轨迹平滑且经过三个指定点。
• 避免拼接问题：直接生成起点到终点的单段弧线，自然过渡经过最高点，无曲率突变。
• 保持绘图逻辑：沿用原有的极坐标转换和绘图代码，确保输出符合需求。

内容的提问来源于stack exchange，提问作者Blake