868 lines
37 KiB
Python
868 lines
37 KiB
Python
#!/usr/bin/env python3
|
||
"""Tianwen-1 frame analysis notebook converted to a Python script.
|
||
|
||
本脚本由 `Tianwen-1 frame analysis.ipynb` 转换而来，保留原 notebook
|
||
的执行顺序，并把原 Markdown 标题和交互式观察步骤整理为注释。
|
||
|
||
脚本主要功能：
|
||
1. 读取 2020-07-30 的天问一号 AOS 帧二进制数据。
|
||
2. 使用本目录下的 `ccsds.py` 解析 AOS 主头、插入区和 Space Packet。
|
||
3. 按 spacecraft ID、virtual channel 和 APID 对数据分组。
|
||
4. 绘制帧计数、时间戳、载荷字节和若干遥测通道的变化曲线。
|
||
|
||
注意：
|
||
- 原 notebook 假设当前工作目录就是 `Tianwen/`。本脚本会自动切换到
|
||
文件所在目录，以便相对路径 `tianwen1_frames_20200730.u8` 能正常工作。
|
||
- notebook 中单独显示表达式结果的单元，在脚本中改为 `print(...)`。
|
||
- 绘图默认使用 matplotlib 当前后端；如果在无图形界面环境运行，可在
|
||
执行前设置 `MPLBACKEND=Agg`。
|
||
"""
|
||
|
||
|
||
# sys 用来修改 Python 的模块搜索路径。脚本在 test/ 目录下，
|
||
# 但要导入 Tianwen/ 目录里的模块，所以需要把项目根目录加进去。
|
||
import sys
|
||
|
||
# pathlib.Path 是处理文件路径的标准库工具。这里取别名 PathlibPath，
|
||
# 是为了避免后面 `from construct import *` 导入的 Path 把它覆盖掉。
|
||
from pathlib import Path as PathlibPath
|
||
|
||
# os 用来切换当前工作目录、设置环境变量。
|
||
import os
|
||
|
||
# re 是正则表达式库，用来把图标题清洗成安全的文件名。
|
||
import re
|
||
|
||
# 获取项目根目录。当前文件在 test/ 下，所以 parent.parent 就是 CCSDS_study/。
|
||
PROJECT_ROOT = PathlibPath(__file__).resolve().parent.parent
|
||
|
||
# 把项目根目录加入 sys.path 后，Python 才能通过 `import Tianwen.xxx`
|
||
# 找到 /home/zjz/CCSDS_study/Tianwen 这个包目录。
|
||
sys.path.append(str(PROJECT_ROOT))
|
||
|
||
# collections.Counter 可以方便地统计某个字段不同取值出现了多少次。
|
||
import collections
|
||
|
||
# struct 用来按指定字节序把 bytes 解析成整数；这里用于解析 packet 内部时间戳。
|
||
import struct
|
||
|
||
# Matplotlib 会写缓存文件。WSL 或某些 sandbox 环境下用户配置目录可能不可写，
|
||
# 因此把缓存目录放到 /tmp/matplotlib，避免启动时出现权限警告。
|
||
os.environ.setdefault("MPLCONFIGDIR", "/tmp/matplotlib")
|
||
|
||
# 先导入 matplotlib 主模块，设置后端后再导入 pyplot。
|
||
import matplotlib
|
||
|
||
# Agg 是非交互式绘图后端，只负责把图片画到文件中，不需要 GUI 窗口。
|
||
# 这对 WSL 很重要，因为 WSL 下直接 plt.show() 经常会因为没有显示服务而报错。
|
||
matplotlib.use("Agg")
|
||
|
||
# pyplot 是 matplotlib 最常用的画图接口，习惯上命名为 plt。
|
||
import matplotlib.pyplot as plt
|
||
|
||
# 脚本会创建很多张图。这个参数关闭“打开超过 20 个 figure”的提示，
|
||
# 因为这里的大量 figure 是预期行为，最后会统一保存并关闭。
|
||
plt.rcParams["figure.max_open_warning"] = 0
|
||
|
||
# numpy 用于高效读取二进制文件、重塑数组、计算差分和处理时间序列。
|
||
import numpy as np
|
||
|
||
# scipy.signal 是 notebook 原始导入；当前脚本没有直接使用，保留它是为了尽量
|
||
# 与原 notebook 环境一致，后续如果扩展信号处理步骤可直接使用。
|
||
import scipy.signal # 保留 notebook 原始导入；当前脚本未显式调用。
|
||
|
||
# construct 是解析二进制协议结构的库。ccsds.py 中的帧结构定义依赖它。
|
||
# 这里使用星号导入是 notebook 原始写法；新项目中通常更建议显式导入需要的名字。
|
||
from construct import * # noqa: F401,F403 - ccsds 解析结构依赖 construct 风格对象。
|
||
|
||
# ccsds.py 定义了 AOSFrame、SpacePacketPrimaryHeader 和 Space Packet 重组逻辑。
|
||
import Tianwen.ccsds as ccsds
|
||
|
||
# tianwen1_tm.py 定义了天问一号时间戳到 numpy.datetime64 的转换函数。
|
||
import Tianwen.tianwen1_tm as tianwen1_tm
|
||
|
||
|
||
# `construct import *` 会导入名为 Path 的对象，因此 pathlib.Path 必须使用别名。
|
||
# notebook 使用相对路径读取 Tianwen 目录下的数据文件；脚本位于 test/ 目录时，
|
||
# 运行前需要切换到真正的数据目录。
|
||
DATA_DIR = PROJECT_ROOT / "Tianwen"
|
||
|
||
# 所有 PNG 图片统一保存到这个目录，避免和源码文件混在一起。
|
||
OUTPUT_DIR = DATA_DIR / "frame_analysis_figures"
|
||
|
||
|
||
def safe_filename(text, max_length=120):
|
||
"""把图标题转换成适合作为文件名的短字符串。"""
|
||
|
||
# 去掉标题首尾空白，并统一转成小写，让文件名风格一致。
|
||
text = text.strip().lower()
|
||
|
||
# 把不适合放进文件名的字符替换为下划线。
|
||
text = re.sub(r"[^0-9a-zA-Z._() -]+", "_", text)
|
||
|
||
# 把连续空白压缩成一个下划线。
|
||
text = re.sub(r"\s+", "_", text)
|
||
|
||
# 去掉文件名开头/结尾的点、下划线、短横线和空格。
|
||
text = text.strip("._- ")
|
||
|
||
# 如果标题为空，就用 untitled；再截断到最大长度，避免文件名过长。
|
||
return (text or "untitled")[:max_length]
|
||
|
||
|
||
def save_all_figures(output_dir):
|
||
"""保存当前 matplotlib 会话中的所有 figure。
|
||
|
||
文件名格式为 `figure_序号_标题.png`。标题来自当前坐标轴的 title；
|
||
如果没有标题，则使用 `untitled`。保存后关闭 figure，降低长脚本的内存占用。
|
||
"""
|
||
|
||
output_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
# saved_paths 用来记录保存出的每一个 PNG 路径，方便最后返回或调试。
|
||
saved_paths = []
|
||
|
||
# plt.get_fignums() 返回当前还打开着的所有 figure 编号。
|
||
for figure_number in plt.get_fignums():
|
||
# 根据编号取回 figure 对象。
|
||
fig = plt.figure(figure_number)
|
||
|
||
# 一个 figure 里可以有一个或多个坐标轴 axes。
|
||
axes = fig.get_axes()
|
||
|
||
# 默认标题为空；后面会尝试从坐标轴里找第一个非空标题。
|
||
title = ""
|
||
|
||
# 遍历坐标轴，拿到用于生成文件名的图标题。
|
||
for ax in axes:
|
||
title = ax.get_title().strip()
|
||
if title:
|
||
break
|
||
|
||
# 文件名前缀使用 figure 序号，保证即使标题重复也不会互相覆盖。
|
||
filename = f"figure_{figure_number:03d}_{safe_filename(title)}.png"
|
||
|
||
# 拼出完整保存路径。
|
||
path = output_dir / filename
|
||
|
||
# 真正保存图片。dpi 控制分辨率，bbox_inches='tight' 尽量裁掉多余白边。
|
||
fig.savefig(path, dpi=150, bbox_inches="tight")
|
||
|
||
# 记录保存路径。
|
||
saved_paths.append(path)
|
||
|
||
# 保存后关闭 figure，释放内存。
|
||
plt.close(fig)
|
||
|
||
# 打印保存数量和目录，方便运行脚本后立刻知道去哪里看图片。
|
||
print(f"Saved {len(saved_paths)} figures to {output_dir}")
|
||
|
||
# 返回路径列表，便于以后如果把 main 拆成函数时继续使用。
|
||
return saved_paths
|
||
|
||
|
||
def get_packet(p):
|
||
"""返回 Space Packet 原始字节。
|
||
|
||
`ccsds.extract_space_packets(..., get_timestamps=True)` 返回
|
||
`(packet_bytes, timestamp)` 二元组；不带时间戳时只返回 packet bytes。
|
||
这个辅助函数屏蔽两种返回格式差异。
|
||
"""
|
||
|
||
# 如果 p 是 tuple，说明它是 `(packet_bytes, timestamp)`；
|
||
# 如果不是 tuple，说明它本身就是 packet bytes。
|
||
return p[0] if type(p) is tuple else p
|
||
|
||
|
||
def packets_asarray(packets):
|
||
"""把 Space Packet 载荷转换为二维 uint8 数组。
|
||
|
||
每个 Space Packet 的前部是 CCSDS Space Packet Primary Header。
|
||
本函数跳过主头，只保留后续用户数据区，便于直接按字节或按数值类型
|
||
重解释为遥测通道矩阵。
|
||
"""
|
||
|
||
# 最终返回一个二维数组：行是 packet，列是 packet 载荷中的字节位置。
|
||
return np.array(
|
||
[
|
||
np.frombuffer(
|
||
# get_packet(p) 先拿到原始 packet bytes；
|
||
# 再跳过 Space Packet Primary Header，只保留用户数据区。
|
||
get_packet(p)[ccsds.SpacePacketPrimaryHeader.sizeof() :],
|
||
"uint8",
|
||
)
|
||
for p in packets
|
||
]
|
||
)
|
||
|
||
|
||
def plot_apids(apids, sc, vc):
|
||
"""按 APID 绘制 Space Packet 载荷热力图。
|
||
|
||
横轴对应载荷字节位置，纵轴对应同一 APID 下的 packet 序号。这个图常用于
|
||
观察哪些字段固定、哪些字段随时间变化，以及是否存在明显的分段结构。
|
||
"""
|
||
|
||
for apid in sorted(apids.keys()):
|
||
# 每个 APID 单独生成一张图。
|
||
plt.figure(figsize=(16, 16), facecolor="w")
|
||
|
||
# 把这个 APID 下所有 packet 载荷堆叠成二维数组。
|
||
ps = packets_asarray(apids[apid])
|
||
|
||
# imshow 把二维字节数组画成热力图，颜色代表字节值。
|
||
plt.imshow(ps, aspect=ps.shape[1] / ps.shape[0])
|
||
|
||
# 标题写清楚 APID、航天器 ID 和虚拟信道 ID，保存成文件后也容易辨认。
|
||
plt.title(f"APID {apid} Spacecraft {sc} Virtual channel {vc}")
|
||
|
||
|
||
def get_timestamps(aos):
|
||
"""从 AOS 帧插入区提取时间戳，并转换为 numpy.datetime64。"""
|
||
|
||
# a.insert_zone.timestamp 是 AOS 帧插入区中的原始时间戳字段。
|
||
return np.array(
|
||
[tianwen1_tm.parse_timestamp_datetime64(a.insert_zone.timestamp) for a in aos]
|
||
)
|
||
|
||
|
||
def get_packet_timestamps(packets):
|
||
"""从 `(packet, timestamp)` 形式的 Space Packet 列表提取时间戳。"""
|
||
|
||
# packets 中的 p[1] 是 extract_space_packets(..., get_timestamps=True)
|
||
# 附带返回的时间戳。
|
||
return np.array([tianwen1_tm.parse_timestamp_datetime64(p[1]) for p in packets])
|
||
|
||
|
||
def main():
|
||
"""按原 notebook 顺序执行完整分析流程。"""
|
||
|
||
# 切换到 Tianwen/ 数据目录后，np.fromfile("tianwen1_frames_20200730.u8")
|
||
# 这种相对路径才能正确找到数据文件。
|
||
os.chdir(DATA_DIR)
|
||
|
||
# -------------------------------------------------------------------------
|
||
# 读取原始 AOS 帧数据
|
||
# -------------------------------------------------------------------------
|
||
# 数据文件由连续的 220 字节帧组成。若文件末尾存在不足一帧的残留字节，
|
||
# 先截断到完整帧边界，再 reshape 成 `帧数 x 220` 的二维数组。
|
||
# 每个 AOS 帧固定 220 字节，这是后面 reshape 的依据。
|
||
frame_size = 220
|
||
|
||
# 从二进制文件中一次性读取所有字节，dtype='uint8' 表示每个元素是 0-255 的字节。
|
||
frames = np.fromfile("tianwen1_frames_20200730.u8", dtype="uint8")
|
||
|
||
# frames.size // frame_size 得到完整帧数量。
|
||
# 乘回 frame_size 后进行切片，可以丢弃文件末尾不满 220 字节的残余数据。
|
||
# reshape((-1, frame_size)) 把一维字节流变成二维表：每一行是一帧。
|
||
frames = frames[: frames.size // frame_size * frame_size].reshape((-1, frame_size))
|
||
|
||
# 打印总帧数，方便确认数据是否读取完整。
|
||
print("Number of frames:", frames.shape[0])
|
||
|
||
# -------------------------------------------------------------------------
|
||
# AOS headers
|
||
# -------------------------------------------------------------------------
|
||
# 解析所有 AOS 帧，得到 primary header、insert zone 等结构化字段。
|
||
# f 是一行 220 字节的 numpy 数组；ccsds.AOSFrame.parse 会按 CCSDS AOS
|
||
# 帧格式把它解释为带字段名的 Container。
|
||
aos = [ccsds.AOSFrame.parse(f) for f in frames]
|
||
|
||
# 统计 AOS transfer frame version number。主流值应集中在有效版本上，
|
||
# 零散异常值通常来自错误同步、误帧或数据噪声。
|
||
print(
|
||
"Transfer frame versions:",
|
||
# Counter 会返回类似 Counter({1: 111790, 2: 147}) 的统计结果。
|
||
collections.Counter([a.primary_header.transfer_frame_version_number for a in aos]),
|
||
)
|
||
|
||
# 统计 spacecraft ID，确认捕获数据中主要包含哪些航天器标识。
|
||
print(
|
||
"Spacecraft IDs:",
|
||
# spacecraft_id 用来区分不同航天器或不同数据源。
|
||
collections.Counter([a.primary_header.spacecraft_id for a in aos]),
|
||
)
|
||
|
||
# 分别观察 spacecraft 245 和 82 下出现的 virtual channel。
|
||
print(
|
||
"SC 245 virtual channels:",
|
||
collections.Counter(
|
||
[
|
||
# virtual_channel_id 是同一 spacecraft 下的数据通道编号。
|
||
a.primary_header.virtual_channel_id
|
||
for a in aos
|
||
# 这里只统计 spacecraft_id 为 245 的帧。
|
||
if a.primary_header.spacecraft_id == 245
|
||
]
|
||
),
|
||
)
|
||
print(
|
||
"SC 82 virtual channels:",
|
||
collections.Counter(
|
||
[
|
||
a.primary_header.virtual_channel_id
|
||
for a in aos
|
||
# 这里只统计 spacecraft_id 为 82 的帧。
|
||
if a.primary_header.spacecraft_id == 82
|
||
]
|
||
),
|
||
)
|
||
|
||
# 打印前 20 帧的 `(spacecraft_id, virtual_channel_id)` 序列，观察复用模式。
|
||
print(
|
||
"First 20 SC/VC pairs:",
|
||
[
|
||
# 每个 tuple 展示一帧的 spacecraft 和 virtual channel。
|
||
(a.primary_header.spacecraft_id, a.primary_header.virtual_channel_id)
|
||
for a in aos[:20]
|
||
],
|
||
)
|
||
|
||
# -------------------------------------------------------------------------
|
||
# Spacecraft 82 Virtual channel 1
|
||
# -------------------------------------------------------------------------
|
||
# 取出 SC 82 / VC 1 的帧，检查帧计数连续性和插入区字段。
|
||
print(
|
||
"SC 82 VC 1 first headers:",
|
||
# 先打印前 10 个 header，人工检查字段是否看起来合理。
|
||
[a.primary_header for a in aos if a.primary_header.spacecraft_id == 82][:10],
|
||
)
|
||
|
||
# 过滤出 spacecraft_id=82 且 virtual_channel_id=1 的所有 AOS 帧。
|
||
sc82_vc1 = [
|
||
a
|
||
for a in aos
|
||
if a.primary_header.spacecraft_id == 82
|
||
and a.primary_header.virtual_channel_id == 1
|
||
]
|
||
|
||
# 提取虚拟信道帧计数。正常情况下相邻帧计数应该连续递增。
|
||
fc = np.array([a.primary_header.virtual_channel_frame_count for a in sc82_vc1])
|
||
|
||
# 新建一张图，用来画帧丢失情况。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
|
||
# np.diff(fc) 是相邻帧计数差；如果连续递增，差值应为 1。
|
||
# 减 1 后，0 表示未丢帧，正数表示中间可能缺了多少帧。
|
||
plt.plot(np.diff(fc) - 1)
|
||
|
||
# 图标题和坐标轴标签会显示在图片上，也会参与生成保存文件名。
|
||
plt.title("Spacecraft 82 Virtual channel 1 frame loss")
|
||
plt.ylabel("Lost frames")
|
||
plt.xlabel("Frame number")
|
||
|
||
print(
|
||
"SC 82 VC 1 first insert-zone timestamps:",
|
||
# hex(...) 把整数时间戳打印成十六进制，方便和二进制协议字段对照。
|
||
[hex(a.insert_zone.timestamp) for a in sc82_vc1][:10],
|
||
)
|
||
print(
|
||
"SC 82 VC 1 insert-zone unknown values:",
|
||
# set 推导式会去重，用来查看 unknown 字段是否固定。
|
||
{hex(a.insert_zone.unknown) for a in sc82_vc1},
|
||
)
|
||
|
||
# 从该虚拟信道重组 Space Packet，并按 APID 分组。
|
||
# AOS 帧里承载的是 Space Packet 数据流；extract_space_packets 会把跨帧的
|
||
# Space Packet 重新拼出来。
|
||
sc82_vc1_packets = list(ccsds.extract_space_packets(sc82_vc1, 82, 1))
|
||
|
||
# 每个 Space Packet 开头都有 Primary Header，这里解析这些 header。
|
||
sc82_vc1_sp_headers = [
|
||
ccsds.SpacePacketPrimaryHeader.parse(p) for p in sc82_vc1_packets
|
||
]
|
||
|
||
# APID 是 Space Packet 的应用进程标识。统计 APID 可以知道有哪些数据类型。
|
||
sc82_vc1_apids = collections.Counter([p.APID for p in sc82_vc1_sp_headers])
|
||
print("SC 82 VC 1 APIDs:", sc82_vc1_apids)
|
||
|
||
# 构造字典：key 是 APID，value 是属于该 APID 的 packet 列表。
|
||
sc82_vc1_by_apid = {
|
||
apid: [
|
||
p
|
||
# zip 把 header 和 packet 一一配对，便于按 header.APID 筛选 packet。
|
||
for h, p in zip(sc82_vc1_sp_headers, sc82_vc1_packets)
|
||
if h.APID == apid
|
||
]
|
||
for apid in sc82_vc1_apids
|
||
}
|
||
|
||
# 为每个 APID 保存一张载荷热力图。
|
||
plot_apids(sc82_vc1_by_apid, 82, 1)
|
||
|
||
# -------------------------------------------------------------------------
|
||
# Spacecraft 245 Virtual channel 3
|
||
# -------------------------------------------------------------------------
|
||
# SC 245 / VC 3 似乎包含固定或填充类数据。这里检查帧计数和若干字节区间。
|
||
# 过滤出 spacecraft_id=245 且 virtual_channel_id=3 的帧。
|
||
sc245_vc3 = [
|
||
a
|
||
for a in aos
|
||
if a.primary_header.spacecraft_id == 245
|
||
and a.primary_header.virtual_channel_id == 3
|
||
]
|
||
|
||
# 提取 VC 3 的帧计数，用于检查是否连续。
|
||
fc = np.array([a.primary_header.virtual_channel_frame_count for a in sc245_vc3])
|
||
|
||
# 绘制 VC 3 的疑似丢帧情况。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(np.diff(fc) - 1)
|
||
plt.title("Spacecraft 245 Virtual channel 3 frame loss")
|
||
plt.ylabel("Lost frames")
|
||
plt.xlabel("Frame number")
|
||
|
||
print("SC 245 VC 3 first headers:", [a.primary_header for a in sc245_vc3[:10]])
|
||
|
||
# 这里需要原始字节帧，而不是解析后的 aos Container，所以从 frames 和 aos
|
||
# 同时遍历，筛选出 SC 245 / VC 3 对应的原始 220 字节帧。
|
||
sc245_vc3_frames = np.array(
|
||
[
|
||
f
|
||
for f, a in zip(frames, aos)
|
||
if a.primary_header.spacecraft_id == 245
|
||
and a.primary_header.virtual_channel_id == 3
|
||
]
|
||
)
|
||
|
||
# np.unique 查看某些字节范围里出现过哪些不同值。
|
||
# 如果只出现一个值，说明该字段可能是固定填充。
|
||
print("SC 245 VC 3 unique bytes 5-12:", np.unique(sc245_vc3_frames[:, 5:12]))
|
||
print("hex(170):", hex(170))
|
||
print("SC 245 VC 3 unique byte 12:", np.unique(sc245_vc3_frames[:, 12]))
|
||
print("hex(138):", hex(138))
|
||
print("SC 245 VC 3 unique bytes 13-end:", np.unique(sc245_vc3_frames[:, 13:]))
|
||
print("SC 245 VC 3 insert-zone unknown:", hex(sc245_vc3[0].insert_zone.unknown))
|
||
|
||
# -------------------------------------------------------------------------
|
||
# Spacecraft 245 Virtual channel 1
|
||
# -------------------------------------------------------------------------
|
||
# SC 245 / VC 1 是后续 APID 分析的主要数据源。
|
||
# 过滤出主要遥测通道 SC 245 / VC 1。
|
||
sc245_vc1 = [
|
||
a
|
||
for a in aos
|
||
if a.primary_header.spacecraft_id == 245
|
||
and a.primary_header.virtual_channel_id == 1
|
||
]
|
||
|
||
# 提取该虚拟信道的帧计数。
|
||
fc = np.array([a.primary_header.virtual_channel_frame_count for a in sc245_vc1])
|
||
|
||
# 绘制该通道的帧丢失情况。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(np.diff(fc) - 1)
|
||
plt.title("Spacecraft 245 Virtual channel 1 frame loss")
|
||
plt.ylabel("Lost frames")
|
||
plt.xlabel("Frame number")
|
||
|
||
print("SC 245 VC 1 first headers:", [a.primary_header for a in sc245_vc1[:10]])
|
||
|
||
# 从 AOS 插入区中取出时间戳，后续画“帧计数随时间变化”的散点图。
|
||
sc245_vc1_timestamps = get_timestamps(sc245_vc1)
|
||
|
||
# 画帧计数随时间的变化。横轴是真实时间，纵轴是 virtual channel frame count。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(sc245_vc1_timestamps, fc, ".")
|
||
plt.title("Spacecraft 245 Virtual channel 1")
|
||
plt.xlabel("Timestamp")
|
||
plt.ylabel("Virtual channel frame count")
|
||
plt.grid()
|
||
|
||
print(
|
||
"SC 245 VC 1 insert-zone unknown values:",
|
||
{hex(a.insert_zone.unknown) for a in sc245_vc1},
|
||
)
|
||
|
||
# 重组 SC 245 / VC 1 中的 Space Packet，并保留 AOS 帧插入区时间戳。
|
||
# get_timestamps=True 表示每个 packet 会附带其所在 AOS 帧的时间戳。
|
||
sc245_vc1_packets = list(
|
||
ccsds.extract_space_packets(sc245_vc1, 245, 1, get_timestamps=True)
|
||
)
|
||
|
||
# 因为 packet 现在是 `(packet_bytes, timestamp)`，所以解析 header 时使用 p[0]。
|
||
sc245_vc1_sp_headers = [
|
||
ccsds.SpacePacketPrimaryHeader.parse(p[0]) for p in sc245_vc1_packets
|
||
]
|
||
|
||
# 统计主遥测通道中每个 APID 出现的 packet 数量。
|
||
sc245_vc1_apids = collections.Counter([p.APID for p in sc245_vc1_sp_headers])
|
||
print("SC 245 VC 1 APIDs:", sc245_vc1_apids)
|
||
|
||
# 按 APID 对 packet 分桶，后续每个 APID 单独分析。
|
||
sc245_vc1_by_apid = {
|
||
apid: [
|
||
p
|
||
for h, p in zip(sc245_vc1_sp_headers, sc245_vc1_packets)
|
||
if h.APID == apid
|
||
]
|
||
for apid in sc245_vc1_apids
|
||
}
|
||
|
||
# 为 SC 245 / VC 1 的每个 APID 生成一张载荷热力图。
|
||
plot_apids(sc245_vc1_by_apid, 245, 1)
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 1280
|
||
# -------------------------------------------------------------------------
|
||
# 将 APID 1280 的载荷按大端 int16 解释，观察前 12 字节的遥测量。
|
||
# packets_asarray(...) 得到 uint8 字节矩阵；view("int16") 把每两个字节看成
|
||
# 一个 16 位整数；byteswap() 用来处理大端/小端字节序差异。
|
||
channels = packets_asarray(sc245_vc1_by_apid[1280]).view("int16").byteswap()
|
||
|
||
# 取出每个 packet 对应的 AOS 帧时间戳，作为绘图横轴。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[1280])
|
||
|
||
# 绘制 APID 1280 前 6 个 int16 通道，也就是原始载荷 bytes 0-12。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, :6])
|
||
plt.title("APID 1280 Spacecraft 245 Virtual channel 1 bytes 0-12 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(0, 12, 2)])
|
||
|
||
# 第二张图只看局部区间，便于放大观察某段时间内的细节变化。
|
||
plt.figure(figsize=(14, 6), facecolor="w")
|
||
|
||
# slice(850, 1200) 表示取第 850 到第 1199 个样本。
|
||
sel = slice(850, 1200)
|
||
plt.plot(timestamps[sel], channels[sel, :6])
|
||
plt.title("APID 1280 Spacecraft 245 Virtual channel 1 bytes 0-12 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(0, 12, 2)])
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 1281
|
||
# -------------------------------------------------------------------------
|
||
# APID 1281 同样按大端 int16 解释，分段查看 18-36 字节。
|
||
# channels 的行是 packet 序号，列是 int16 通道序号。
|
||
channels = packets_asarray(sc245_vc1_by_apid[1281]).view("int16").byteswap()
|
||
|
||
# timestamps 和 channels 的行一一对应。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[1281])
|
||
|
||
# 查看 int16 通道 9 到 14，对应字节范围 18-30。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, 9:15])
|
||
plt.title("APID 1281 Spacecraft 245 Virtual channel 1 bytes 18-30 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(9 * 2, 15 * 2, 2)])
|
||
|
||
# 放大 APID 1281 的一小段样本。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
sel = slice(400, 600)
|
||
plt.plot(timestamps[sel], channels[sel, 9:15])
|
||
plt.title("APID 1281 Spacecraft 245 Virtual channel 1 bytes 18-30 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(9 * 2, 15 * 2, 2)])
|
||
|
||
# 查看 int16 通道 15 到 17，对应字节范围 30-36。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, 15:18])
|
||
plt.title("APID 1281 Spacecraft 245 Virtual channel 1 bytes 30-36 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(2 * 15, 18 * 2, 2)])
|
||
|
||
# 同样对 bytes 30-36 的局部区间放大观察。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps[sel], channels[sel, 15:18])
|
||
plt.title("APID 1281 Spacecraft 245 Virtual channel 1 bytes 30-36 (int16)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(2 * 15, 18 * 2, 2)])
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 1288
|
||
# -------------------------------------------------------------------------
|
||
# APID 1288 的载荷偏移 80 后，前 24 字节按三个大端 float64 解释。
|
||
# notebook 中把最大时间间隔对应行置为 NaN，用于清理明显观测间断。
|
||
# [:, 80:] 表示丢掉载荷前 80 字节；[:, : 3 * 8] 表示再取 24 字节。
|
||
# 24 字节正好可以看成 3 个 float64。
|
||
channels = (
|
||
np.copy(packets_asarray(sc245_vc1_by_apid[1288])[:, 80:][:, : 3 * 8])
|
||
.view("float64")
|
||
.byteswap()
|
||
)
|
||
|
||
# APID 1288 的每个 packet 时间戳。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[1288])
|
||
|
||
# np.diff(timestamps) 找相邻样本的时间间隔；argmax 找最大间隔位置。
|
||
# 最大间隔可能代表观测中断，因此把该点设成 NaN，绘图时会断开。
|
||
channels[np.argmax(np.diff(timestamps)), :] = np.nan
|
||
|
||
# 绘制三个 float64 通道的原始变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels)
|
||
plt.title("APID 1288 Spacecraft 245 Virtual channel 1 bytes 80-104 (float64)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (float64)")
|
||
plt.legend([f"Bytes {j}-{j + 8}" for j in range(80, 104, 8)])
|
||
plt.grid()
|
||
|
||
# 对每个 float64 通道拟合 1 阶和 2 阶多项式，并绘制去趋势残差。
|
||
# 去趋势的目的：把总体平滑变化减掉，只看剩余的小幅波动。
|
||
for deg in [1, 2]:
|
||
# 每个多项式阶数单独画一张图。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
|
||
# channels.shape[1] 是通道数量，这里通常是 3。
|
||
for j in range(channels.shape[1]):
|
||
# 把 datetime64 转成“相对起点的秒数”，polyfit 需要普通数值横轴。
|
||
s = (timestamps - timestamps[0]) / np.timedelta64(1, "s")
|
||
|
||
# 过滤 NaN 点，否则 polyfit 无法拟合。
|
||
sel = ~np.isnan(channels[:, j])
|
||
|
||
# np.polyfit 拟合 deg 阶多项式，返回多项式系数。
|
||
p = np.polyfit(s[sel], channels[sel, j], deg)
|
||
|
||
# np.polyval 计算拟合曲线；原始值减拟合曲线就是残差。
|
||
plt.plot(timestamps, channels[:, j] - np.polyval(p, s))
|
||
plt.title("APID 1288 Spacecraft 245 Virtual channel 1 bytes 80-104 (float64)")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (float64)")
|
||
plt.legend(
|
||
[
|
||
f"Bytes {j}-{j + 8} minus degree {deg} polynomial"
|
||
for j in range(80, 104, 8)
|
||
]
|
||
)
|
||
plt.grid()
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 1287
|
||
# -------------------------------------------------------------------------
|
||
# APID 1287 的载荷偏移 66 后，前 24 字节按三个大端 float64 解释。
|
||
# 另外从 packet 内部字节 18 起取 6 字节时间字段，补两个高位零后按大端
|
||
# uint64 解析，与 AOS 插入区时间戳比较。
|
||
# 这里的处理方式和 APID 1288 类似，只是 float64 起始偏移变成 66。
|
||
channels = (
|
||
np.copy(packets_asarray(sc245_vc1_by_apid[1287])[:, 66:][:, : 3 * 8])
|
||
.view("float64")
|
||
.byteswap()
|
||
)
|
||
|
||
# AOS 帧插入区提供的时间戳。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[1287])
|
||
|
||
# Space Packet 自身载荷中也有一个 6 字节时间字段。
|
||
timestamps_packet = np.array(
|
||
[
|
||
tianwen1_tm.parse_timestamp_datetime64(
|
||
# b"\x00\x00" + 6 字节时间字段 = 8 字节，才能用 >Q 解析为 uint64。
|
||
struct.unpack(">Q", b"\x00\x00" + p[0][18:][:6])[0]
|
||
)
|
||
for p in sc245_vc1_by_apid[1287]
|
||
]
|
||
)
|
||
|
||
# 清理最大时间间隔处的样本，避免图中出现不合理连接。
|
||
channels[np.argmax(np.diff(timestamps)), :] = np.nan
|
||
|
||
# 比较 AOS 帧时间戳和 Space Packet 内部时间戳之间的差值。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
|
||
# 时间戳相减得到 timedelta64，再除以 1 秒单位，得到浮点秒数。
|
||
t_diff = (timestamps - timestamps_packet) / np.timedelta64(1, "s")
|
||
|
||
# 同样在最大间隔处断开曲线。
|
||
t_diff[np.argmax(np.diff(timestamps))] = np.nan
|
||
plt.plot(timestamps, t_diff)
|
||
plt.title("APID 1287 Spacecraft 245 Virtual channel 1 timestamp difference")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("AOS frame timestamp - Space packet timestamp")
|
||
|
||
# 绘制 APID 1287 三个 float64 通道，横轴使用 packet 内部时间戳。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps_packet, channels)
|
||
plt.title("APID 1287 Spacecraft 245 Virtual channel 1 bytes 66-90 (float64)")
|
||
plt.xlabel("Space Packet frame timestamp")
|
||
plt.ylabel("Raw value (float64)")
|
||
plt.legend([f"Bytes {j}-{j + 8}" for j in range(66, 90, 8)])
|
||
plt.grid()
|
||
|
||
# 分别对全时段和 2020-07-31 03:00 前的数据做 1 阶、2 阶去趋势分析。
|
||
# cut 是满足 timestamp <= 2020-07-31T03:00 的最后一个索引。
|
||
cut = np.where(timestamps_packet <= np.datetime64("2020-07-31T03:00"))[0][-1]
|
||
|
||
# 第一个 slice(None, None) 表示全时段；第二个 slice(0, cut) 表示前半段。
|
||
for selt in [slice(None, None), slice(0, cut)]:
|
||
for deg in [1, 2]:
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
for j in range(channels.shape[1]):
|
||
# 使用当前选中时间段的起点作为 0 秒。
|
||
s = (
|
||
timestamps_packet[selt] - timestamps_packet[selt][0]
|
||
) / np.timedelta64(1, "s")
|
||
|
||
# 只用非 NaN 数据点参与拟合。
|
||
sel = ~np.isnan(channels[selt, j])
|
||
|
||
# 对当前通道和当前时间段拟合多项式。
|
||
p = np.polyfit(s[sel], channels[selt][sel, j], deg)
|
||
|
||
# 绘制去趋势后的残差。
|
||
plt.plot(timestamps_packet[selt], channels[selt, j] - np.polyval(p, s))
|
||
plt.title("APID 1287 Spacecraft 245 Virtual channel 1 bytes 66-90 (float64)")
|
||
plt.xlabel("Space Packet timestamp")
|
||
plt.ylabel("Raw value (float64)")
|
||
plt.legend(
|
||
[
|
||
f"Bytes {j}-{j + 8} minus degree {deg} polynomial"
|
||
for j in range(66, 90, 8)
|
||
]
|
||
)
|
||
plt.grid()
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APIDs 2 and 3
|
||
# -------------------------------------------------------------------------
|
||
# 原 notebook 这里两个变量都从 APID 2 取数，保留原逻辑。若要比较 APID 2
|
||
# 与 APID 3，可把 `sc245_vc1_by_apid[2]` 的第二处改为 `[3]`。
|
||
# channels2 是 APID 2 的 uint8 载荷矩阵。
|
||
channels2 = packets_asarray(sc245_vc1_by_apid[2])
|
||
|
||
# timestamps2 是 APID 2 每个 packet 对应的 AOS 时间戳。
|
||
timestamps2 = get_packet_timestamps(sc245_vc1_by_apid[2])
|
||
|
||
# 注意：这里原 notebook 仍然取 APID 2，不是真正的 APID 3。
|
||
# 为了忠实保留原分析脚本，暂不改动。
|
||
channels3 = packets_asarray(sc245_vc1_by_apid[2])
|
||
|
||
# 同上，这里的 timestamps3 也来自 APID 2。
|
||
timestamps3 = get_packet_timestamps(sc245_vc1_by_apid[2])
|
||
|
||
# 对比 byte 0 的变化；'.-' 是点加线，'.' 是只画点。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps2, channels2[:, 0], ".-", label="APID 2 byte 0")
|
||
plt.plot(timestamps3, channels3[:, 0], ".", label="APID 3 byte 0")
|
||
plt.title("APIDs 2 and 3 Spacecraft 245 Virtual channel 1 byte 0")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.legend()
|
||
plt.grid()
|
||
|
||
# 对比 byte 2 的变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps2, channels2[:, 2], ".-", label="APID 2 byte 2")
|
||
plt.plot(timestamps3, channels3[:, 2], ".", label="APID 3 byte 2")
|
||
plt.title("APIDs 2 and 3 Spacecraft 245 Virtual channel 1 byte 2")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.legend()
|
||
plt.grid()
|
||
|
||
# 只取局部样本，放大查看 byte 4 的变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
|
||
# 从第 14350 个样本开始取 400 个样本。
|
||
sel = slice(14350, 14350 + 400)
|
||
plt.plot(timestamps2[sel], channels2[sel, 4], ".-", label="APID 2 byte 4")
|
||
plt.plot(timestamps3[sel], channels3[sel, 4], ".", label="APID 3 byte 4")
|
||
plt.title("APIDs 2 and 3 Spacecraft 245 Virtual channel 1 byte 4")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.legend()
|
||
plt.grid()
|
||
|
||
# 使用同一个局部范围，查看 byte 8 的变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps2[sel], channels2[sel, 8], ".-", label="APID 2 byte 8")
|
||
plt.plot(timestamps3[sel], channels3[sel, 8], ".", label="APID 3 byte 8")
|
||
plt.title("APIDs 2 and 3 Spacecraft 245 Virtual channel 1 byte 8")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.legend()
|
||
plt.grid()
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 454
|
||
# -------------------------------------------------------------------------
|
||
# APID 454 按 uint8 查看，重点观察 byte 0 以及 byte 4-7、18-20。
|
||
# 这里没有 view 成 int16/float64，说明按单字节原始值观察。
|
||
channels = packets_asarray(sc245_vc1_by_apid[454])
|
||
|
||
# APID 454 的 packet 时间戳。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[454])
|
||
|
||
# 放大一小段样本，查看 byte 0。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
sel = slice(1780, 1850)
|
||
plt.plot(timestamps[sel], channels[sel, 0], ".-")
|
||
plt.title("APID 454 Spacecraft 245 Virtual channel 1 byte 0")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.grid()
|
||
|
||
# 在同一张图中画 byte 4-6 和 byte 18-19，比较这些字段是否同步变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps[sel], channels[sel, 4:7], ".-")
|
||
plt.plot(timestamps[sel], channels[sel, 18:20], ".-")
|
||
plt.title("APID 454 Spacecraft 245 Virtual channel 1 bytes 4-7, 18-20")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (uint8)")
|
||
plt.legend([f"Byte {j}" for j in list(range(4, 7)) + list(range(18, 20))])
|
||
plt.grid()
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 462
|
||
# -------------------------------------------------------------------------
|
||
# APID 462 按大端 int16 解释，观察前 14 字节。
|
||
# 前 14 字节对应 7 个 int16 通道。
|
||
channels = packets_asarray(sc245_vc1_by_apid[462]).view("int16").byteswap()
|
||
|
||
# APID 462 的时间戳。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[462])
|
||
|
||
# 绘制 7 个 int16 通道的全时段变化。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, :7])
|
||
plt.title("APID 462 Spacecraft 245 Virtual channel 1 bytes 0-14")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.legend([f"Bytes {j} and {j + 1}" for j in range(0, 14, 2)])
|
||
plt.grid()
|
||
|
||
# -------------------------------------------------------------------------
|
||
# APID 1536
|
||
# -------------------------------------------------------------------------
|
||
# APID 1536 按大端 int16 解释，分别观察 bytes 18-20 和 20-22。
|
||
# int16 通道 9 对应 bytes 18-20，通道 10 对应 bytes 20-22。
|
||
channels = packets_asarray(sc245_vc1_by_apid[1536]).view("int16").byteswap()
|
||
|
||
# APID 1536 的时间戳。
|
||
timestamps = get_packet_timestamps(sc245_vc1_by_apid[1536])
|
||
|
||
# 绘制 int16 通道 9。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, 9])
|
||
plt.title("APID 1536 Spacecraft 245 Virtual channel 1 bytes 18-20")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.grid()
|
||
|
||
# 绘制 int16 通道 10。
|
||
plt.figure(figsize=(12, 6), facecolor="w")
|
||
plt.plot(timestamps, channels[:, 10])
|
||
plt.title("APID 1536 Spacecraft 245 Virtual channel 1 bytes 20-22")
|
||
plt.xlabel("AOS frame timestamp")
|
||
plt.ylabel("Raw value (int16)")
|
||
plt.grid()
|
||
|
||
# notebook 里依赖自动显示图像，但脚本在 WSL / 无显示环境下不应强制调用
|
||
# `plt.show()`，否则可能因为缺少图形界面后端而报错。这里统一用 savefig
|
||
# 保存所有 figure，方便在 WSL 中直接查看生成的 PNG 文件。
|
||
# OUTPUT_DIR 是 Tianwen/frame_analysis_figures。
|
||
save_all_figures(OUTPUT_DIR)
|
||
|
||
|
||
# 只有直接运行这个脚本时才会执行 main()。
|
||
# 如果以后从别的 Python 文件 import 这个脚本，main() 不会自动运行。
|
||
if __name__ == "__main__":
|
||
main()
|