【2025年认证杯数学中国数学建模网络挑战赛】A题解题思路与模型代码

news2025/4/19 2:03:58

【2025年认证杯数学建模挑战赛】A题

该题为典型的空间几何建模+轨道动力学建模+预测问题。

⚙ 问题一：利用多个天文台的同步观测，确定小行星与地球的相对距离

问题分析

已知若干地面天文台的观测数据：方位角 (Azimuth) 和 高度角 (Altitude)；
所有天文台在同一时间观测同一颗小行星；
目标是根据这些观测结果三角测量出小行星的空间坐标，进而得出相对于地球的位置（特别是距离）。

解题思路

获取每个天文台的地理坐标（经纬度、海拔）；
将每个天文台的观测方位角、高度角转换为单位观测方向向量（以地心为参考）；
将地面站的位置用地心坐标表示；
利用两台或多台天文台的观测方向和坐标，建立**多射线最近点法（最小距离点）**来确定目标的小行星的位置。

可能用到的数学模型与方法

地理坐标系 → 地心坐标系转换；
球面坐标系 → 笛卡尔单位向量；
多射线几何定位（三维空间中线与线之间的最短距离问题）；
最小二乘法求解最优交点。

初步代码框架（Python）

import numpy as np

def geo_to_ecef(lat, lon, alt):
    # WGS84参数
    a = 6378137.0  # 地球长半轴
    e = 8.1819190842622e-2  # 偏心率
    lat, lon = np.radians(lat), np.radians(lon)
    N = a / np.sqrt(1 - e**2 * np.sin(lat)**2)
    x = (N + alt) * np.cos(lat) * np.cos(lon)
    y = (N + alt) * np.cos(lat) * np.sin(lon)
    z = (N * (1 - e**2) + alt) * np.sin(lat)
    return np.array([x, y, z])

def azel_to_direction(az, el, lat, lon):
    az, el = np.radians(az), np.radians(el)
    x = np.cos(el) * np.sin(az)
    y = np.cos(el) * np.cos(az)
    z = np.sin(el)
    # 本地坐标转为地心坐标（简略版）
    return np.array([x, y, z])  # 需进一步旋转变换

def triangulate(positions, directions):
    # 最小二乘法解多条射线的最近点
    A, b = [], []
    for p, d in zip(positions, directions):
        d = d / np.linalg.norm(d)
        I = np.eye(3)
        A.append(I - np.outer(d, d))
        b.append((I - np.outer(d, d)) @ p)
    A, b = np.sum(A, axis=0), np.sum(b, axis=0)
    return np.linalg.lstsq(A, b, rcond=None)[0]

# 示例：三台天文台
stations = [
    {'lat': 34.0, 'lon': -118.0, 'alt': 1000, 'az': 130, 'el': 45},
    {'lat': 40.0, 'lon': -75.0,  'alt': 500,  'az': 120, 'el': 40},
    {'lat': 35.0, 'lon': 139.0,  'alt': 20,   'az': 150, 'el': 50},
]

positions, directions = [], []
for s in stations:
    pos = geo_to_ecef(s['lat'], s['lon'], s['alt'])
    dir = azel_to_direction(s['az'], s['el'], s['lat'], s['lon'])  # 简化
    positions.append(pos)
    directions.append(dir)

asteroid_pos = triangulate(positions, directions)
distance = np.linalg.norm(asteroid_pos) / 1.496e+11  # 转为AU单位
print("小行星距离地球中心：", distance, "AU")

⚙ 问题二：利用多次观测数据，预测小行星未来15或30天的轨道及天文台观测角

问题分析

输入：多台天文台对某颗小行星的多时刻观测数据（时刻、方位角、高度角）；
输出：未来15或30天内，每天某固定时刻的方位角和高度角预测值。

解题思路

将观测数据反演为轨道初始状态（即轨道六要素或位置-速度）；
根据轨道力学模型（两体问题/摄动模型），进行轨道积分，得到未来轨道；
对于预测日的每个时间点，将轨道位置转换成各天文台的视角坐标，计算方位角和高度角。

可能用到的数学模型与方法

Gauss法或Laplace法进行轨道确定（初轨道估计）；
开普勒轨道传播（Kepler Propagation）；
JPL HORIZONS、NASA SPICE工具可模拟（也可调用简化库如 Skyfield）；
天球坐标计算公式（向量到方位角、高度角转换）；

初步代码框架（Skyfield库）

from skyfield.api import load, Topos

ts = load.timescale()
ephemeris = load('de421.bsp')
earth = ephemeris['earth']

# 假设已知小行星轨道初值（可从MPC或观测数据拟合）
from skyfield.keplerlib import KeplerOrbit
# 伪代码: 构造轨道
# orbit = KeplerOrbit(elements, epoch=ts.utc(2024, 1, 1))

# 假设某台站
station = earth + Topos(latitude_degrees=34.0, longitude_degrees=-118.0, elevation_m=1000)

# 预测未来30天
t = ts.utc(2024, 1, range(1, 31))
# astrometric = station.at(t).observe(orbit)  # 实际需从轨道构造天体
# alt, az, d = astrometric.apparent().altaz()

# print("每天12:00的高度角：", alt.degrees)
# print("每天12:00的方位角：", az.degrees)

✅ 总结建议

问题	本质	方法建议	数据需求
问题一	三维空间定位	三角测量 + 多射线最短距离	多台站的同时观测数据
问题二	轨道传播与预测	轨道反演 + 两体运动传播 + 坐标转换	多次观测 + 天文台坐标