天和空间站凌月凌日时间计算的尝试

本文主要是介绍天和空间站凌月凌日时间计算的尝试，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

最近对空间站过境、凌日月时间计算产生了兴趣，之前写了空间站过境的代码，接着该空间站凌月时间计算了。在网上搜索没有现成的，暂时也没找到相关算法的介绍。网上找到只有这个可以查询国际空间站凌日月时间的网站。找不到就自己尝试办法吧。凌月无非就是在地面观察点看到两者的位置重合了，两条思路：1.空间站和月球在天球的赤经和赤纬差值极小即可，2.两者从观察点看的夹角为0，根据这两条思路尝试写代码。首先可以把查询时间段内能过境本地的时间筛选出来，这个已经在空间站过境时间计算的代码里实现，并且可已简化计算，只要把过境最高点大于10°（太低了看不到）的升降到地平线时的时间算出来即可。然后分别计算各时间段内两者赤经赤纬差值最小值或夹角最小值。值越小说明两者距离越近。
下面是自己尝试写的代码，代码比较乱，没怎么整理，计算后在stellarium上模拟。（编程和天文都是纯业余爱好，轻喷）

"""
使用skyfield计算空间站凌日月时间
delcomp 2021-06-25
"""
# Topos 已弃用 建议使用wgs84
from skyfield.api import load, wgs84
from skyfield.searchlib import find_maxima, find_minima, _find_discrete, find_discrete
from skyfield.trigonometry import position_angle_of
from skyfield.constants import tau, DAY_S
from pytz import timezone
import numpy as np
import os
import time
from datetime import datetime, timedelta
from my_function import BD09_to_GCJ02, LatLng2DMS, GCJ02_to_GPS84, BD09_to_WGS84start_time = time.perf_counter()# 百度地图经纬度
LAT = xx.xx1
LNG = xxx.xxSTATION = 'TIANHE'
# STATION = 'ISS (ZARYA)'
DAYS = 300
ELEVATION = 0  # 观察点高度
ALT_DEF = 10  # 指定的最低高度
DT_FSTR = '%Y-%m-%d %H:%M:%S'  # 时间格式化串beijing = timezone('Asia/Shanghai')  # 中国时区
stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
stations_file = 'stations.txt'
ts = load.timescale()
# satellites = load.tle_file(stations_url, reload=True)
if not os.path.exists(stations_file):print('本地文件不存在，下载stations文件...')satellites = load.tle_file(stations_url, reload=True)
else:satellites = load.tle_file(stations_file)print('加载本地stations文件！')# 根据卫星名称构建卫星字典by_name = {sat.name: sat for sat in satellites}tianhe = by_name[STATION]# 文件旧了，更新t = ts.now()days = t - tianhe.epochprint('{:.3f} days away from epoch'.format(days))if abs(days) > 0.5:satellites = load.tle_file(stations_url, reload=True)print('stations文件已更新！')# 根据卫星编号构建卫星字典
# by_number = {sat.model.satnum: sat for sat in satellites}
# # 25544：ISS
# # 48274：TIANHE
# satellite = by_number[25544]by_name = {sat.name: sat for sat in satellites}
tianhe = by_name[STATION]print(tianhe)
print('BD09坐标: %.5f %.5f' % (LAT, LNG))
# 转为GCJ02坐标
lat_gcj, lng_gcj = BD09_to_GCJ02(LAT, LNG)
lat_d, lat_m, lat_s = LatLng2DMS(lat_gcj)
lon_d, lon_m, lon_s = LatLng2DMS(lng_gcj)
print('GCJ02坐标: %.5f %.5f' % (lat_gcj, lng_gcj),'(%d°%d′%.1f″  %d°%d′%.1f″)' % (lat_d, lat_m, lat_s, lon_d, lon_m, lon_s))
# 转为WGS84
lat_wgs, lng_wgs = BD09_to_WGS84(LAT, LNG)
lat_d, lat_m, lat_s = LatLng2DMS(lat_wgs)
lon_d, lon_m, lon_s = LatLng2DMS(lng_wgs)
print('WGS84坐标: %.5f %.5f' % (lat_wgs, lng_wgs),'(%d°%d′%.1f″  %d°%d′%.1f″)' % (lat_d, lat_m, lat_s, lon_d, lon_m, lon_s))mypos = wgs84.latlon(lat_wgs, lng_wgs, elevation_m=ELEVATION)
eph = load('de/de421.bsp')
sun, earth, moon = eph['sun'], eph['earth'], eph['moon']# 使用本地时区
now = beijing.localize(datetime.now())
# t0 = now.replace(year= 2022, month=1, day=14, hour=0, minute=0, second=0, microsecond=0)
t0 = now.replace(hour=0, minute=0, second=0, microsecond=0)
t1 = t0 + timedelta(days=DAYS)
print('查询时间间隔(北京)：%s 至 %s' % (t0, t1))
# find_events参数为UTC，将北京时间转为UTC
t0 = ts.from_datetime(t0)
t1 = ts.from_datetime(t1)
print('查询时间间隔(UTC)：%s 至 %s' % (t0.utc_strftime(DT_FSTR), t1.utc_strftime(DT_FSTR)))
difference= tianhe - myposdef cheat(t):"""Avoid computing expensive values that cancel out anyway."""t.gast = t.tt * 0.0t.M = t.MT = np.identity(3)def below_horizon_at(t):cheat(t)return difference.at(t).altaz()[0].degrees < ALT_DEFdef altitude_at(t):cheat(t)return difference.at(t).altaz()[0].degrees# 找到符合条件（中天高度至少10度）的过境时间
half_second = 0.5 / DAY_S
orbits_per_minute = tianhe.model.no_kozai / tau
orbits_per_day = 24 * 60 * orbits_per_minute
step_days = 0.05 / max(orbits_per_day, 1.0)
if step_days > 0.25:step_days = 0.25
# step_days = 1 / DAY_S
altitude_at.step_days = step_days
# 找到最高点的时间及高度
tmax, altitude = find_maxima(t0, t1, altitude_at, half_second, 12)# 筛选高度>10的tmax 对应_find_discrete(高度<10)，
tmax = tmax.tt[altitude >= ALT_DEF]doublets_0 = np.repeat(np.concatenate(((t0.tt,), tmax, (t1.tt,))), 2)
jdo_0 = (doublets_0[:-1] + doublets_0[1:]) / 2.0# 在地平线的时间，一次升起(event=0)，一次落下(event=1)
# 注意：可能碰到一种特殊情况，在指定的时间段内，有可能第一个时间为降落至10°，即event=1 例如：[1,0,1,0,1,0,1,0,1]
# 或者最后一个时间为升高至10°，即event=0，例如：[0,1,0,1,0,1,0,1,0]
# 当上述两种情况，不同时出现时，在后面计算pass_pairs时，无法配对过境时间
# 如果两种情况同时存在时，[1,0,1,0,1,0,1,0,1,0]，可以配对过境时间(升-降，0-1)，但配对将是错误的, [1-0, 1-0, 1-0, ... ]
# 正确配对应该是[0-1, 0-1, 0-1, ... ]
trs_0, event_0 = _find_discrete(t0.ts, jdo_0, below_horizon_at, half_second, 8)
# 判读是否出现上面注释提到的特殊情况，删除
if event_0[0] == 1:  # 第一种情况trs_0 = trs_0[1:]
if event_0[-1] == 0:  # 第二种情况trs_0 = trs_0[:-1]
# 将符合条件的过境时间配对
pass_pairs = trs_0.reshape(len(trs_0) // 2, 2)
# exit(-1)
# 计算在这些时间区间内是否凌月
pos = earth + mypos# def calc_transit(t):
#     m = pos.at(t).observe(moon).apparent()
#     sat = pos.at(t).observe(earth + tianhe).apparent()
#     return position_angle_of(m.altaz(), sat.altaz()).degrees# 角距
def calc_transit_moon(t):m = pos.at(t).observe(moon)sat = pos.at(t).observe(earth + tianhe)return m.separation_from(sat).degreesdef calc_transit_sun(t):m = pos.at(t).observe(sun)sat = pos.at(t).observe(earth + tianhe)return m.separation_from(sat).degrees# 赤经赤纬差
# def calc_transit_1(t):
#     moon_ra, moon_dec, _ = pos.at(t).observe(moon).apparent().radec()
#     station_ra, stations_dec, _ = pos.at(t).observe(earth + tianhe).apparent().radec()
#     return abs(moon_ra.hours - station_ra.hours) + abs(moon_dec.degrees - stations_dec.degrees)print('总过境次数', len(pass_pairs))
calc_transit_moon.step_days = 1 / DAY_S
for p in pass_pairs:tm0 = ts.tt_jd(p[0])tm1 = ts.tt_jd(p[1])times, angle = find_minima(tm0, tm1, calc_transit_moon)for t, ang in zip(times, angle):# 挑选夹角小于5的显示出来，角度太大的肯定不会凌月if ang < 1:t_str = str(t.astimezone(beijing))[:19]alz, az, dist = difference.at(t).altaz()unit = '°'print('{} \t角距:{:>5.3f}{} \t高度:{:>5.1f}° \t方位角:{:>6.1f}° \t距离:{:>7.1f}km'.format(t_str, ang, unit, alz.degrees, az.degrees, dist.km))print('elapsed time:', time.perf_counter() - start_time)

简单说明一下，代码里查询范围是当日至以后300天，这个明显是不可能的，空间站的轨道信息每天都会更新，所以关于空间站的计算，即使没有大的变轨操作，计算结果在1周后误差也会越来越大，这里只当是做个数学游戏。计算后在stellarium里看看是不是符合。另外在百度地图里查到的坐标，需要转换为WGS84坐标。在stellarium里设定自己坐标时，填这个转换后的坐标。
下面是运行结果：
在这里插入图片描述