火箭残骸定位：多源异构数据融合与物理约束建模

1. 这道题到底在考什么：剥离“定位”表象，看清A题的真实命题内核

2024年“深圳杯”数学建模挑战赛A题——“多个火箭残骸的准确定位”，光看标题，很多人第一反应是：“不就是用GPS或者测距算法算坐标吗？”我带过六届校队、审过上百份深圳杯和国赛论文，见过太多队伍一上来就猛扎进Matlab画图、Python调scipy.optimize，结果三天后发现模型跑不通、结果离谱、连基本物理约束都违反。这道题的陷阱，恰恰藏在“准确定位”这四个字背后。

它根本不是一道单纯的“测量学应用题”，而是一道典型的多源异构数据融合+强物理约束下的逆问题求解。你手里没有标准GPS信号，没有高精度惯导，只有一组极其有限、噪声极大、甚至部分缺失的地面观测数据——比如某时刻某雷达站测得的残骸仰角、某光学站拍到的模糊像素坐标、某声呐阵列记录的冲击波到达时间差。这些数据来自不同原理、不同精度、不同坐标系、不同采样频率的传感器，且每个传感器自身还存在系统偏差（比如雷达天线安装倾角误差±0.5°，光学镜头畸变未标定，声速随温湿度变化导致时间换算偏差）。更关键的是，残骸本身在高速下落、翻滚、解体，其运动轨迹不是匀速直线，而是受重力、空气阻力、气流扰动、残骸形状不对称性影响的非线性动力学过程。

所以，真正的核心矛盾是：如何在观测数据稀疏、噪声大、来源杂、模型不确定的前提下，反推一个符合真实物理规律的、可验证的、鲁棒性强的残骸空间轨迹？它考的不是你会不会用最小二乘，而是你能不能判断：什么时候该信数据，什么时候该信模型；当数据冲突时，是剔除异常值，还是修正传感器偏差；当模型预测与观测偏差较大时，是调整参数，还是重构动力学方程。

我去年帮一支队伍复盘他们2023年A题的失败案例：他们用纯几何交会法，把三个雷达站的仰角数据直接代入球面三角公式，算出的落点误差超过8公里。后来我们一行行检查原始数据，发现其中一台雷达的海拔高度输入错了——它装在一座小山上，但队员直接用了城市平均海拔。一个120米的高程误差，在仰角为15°时，会导致水平距离计算偏差近450米。而他们没做任何敏感性分析，就把这个错误当作“正常噪声”滤掉了。这就是典型地把“定位问题”当成了“计算问题”，忽略了建模的第一步：数据可信度评估与误差溯源。

因此，本题的破题起点，必须从构建一个“可解释、可诊断、可迭代”的误差框架开始。不是先写代码，而是先画一张表，列出所有可能的误差源及其量级：

这张表，就是你整个建模工作的“宪法”。后面所有算法选择、参数设置、结果验证，都必须回到这张表里找依据。比如，当你发现残骸下落末段速度很快，而光学站帧率只有10Hz，那你就必须承认：在这个阶段，光学数据的空间分辨率已经不足以支撑亚米级定位，它的价值更多是提供姿态信息或辅助时间对齐，而不是作为主定位源。这种基于物理和工程实际的判断，远比调参重要得多。

2. 数据预处理：为什么90%的队伍在这里栽跟头，而高手只用三步就完成清洗

绝大多数参赛队在第一天就卡死在数据读取和绘图上，花八小时把Excel表格导入Python，用matplotlib画出一堆散点图，然后对着图发呆：“这数据怎么这么乱？”——这不是数据的问题，是你没理解“预处理”的本质。它不是数据清洗，而是建立观测数据与物理世界之间的映射关系。我见过最高效的队伍，用不到两小时就完成了全部预处理，并输出了一份《数据质量白皮书》，这份文档直接决定了他们后续模型的成败。

2.1 第一步：时空基准统一——不是格式转换，而是建立共同语言

原始数据绝不会给你一个整齐的CSV，里面混着：

• 雷达站A：时间戳（本地PC时间，毫秒级，但未校准），仰角（度），方位角（度），距离（米）
• 光学站B：时间戳（NTP服务器时间，秒级），图像中心像素坐标（x,y），镜头焦距（mm），传感器尺寸（mm）
• 声呐阵列C：时间戳（嵌入式系统内部计数器，无绝对时间），各通道到达时间差（微秒）

你以为第一步是“把时间都转成datetime”？错。第一步是识别并标记所有时间源的不确定性。PC时间可能漂移±500ms/天，NTP在局域网内精度约±10ms，嵌入式计数器需要通过一个已知事件（如一次闪光）来标定其与绝对时间的关系。高手的做法是：立刻新建一个time_alignment.py文件，只做一件事——定义一个“观测事件时间轴”，所有数据都投影到这个轴上，并标注每个投影的置信区间。

# time_alignment.py 核心逻辑（示意） import numpy as np from datetime import datetime, timedelta class ObservationTimeline: def __init__(self): # 主时间源：NTP服务器（置信度最高） self.ntp_ref = datetime(2024, 9, 15, 12, 0, 0) self.ntp_uncertainty = timedelta(milliseconds=10) # ±10ms # PC时间源：需校准 self.pc_drift_rate = 0.0005 # 每秒漂移0.5ms，需实测 self.pc_ref_event = {'time_ntp': datetime(2024, 9, 15, 12, 1, 30), 'time_pc': 1678901234567} # 一次同步事件 # 计数器源：需标定 self.counter_ref = {'event_time_ntp': datetime(2024, 9, 15, 12, 2, 15), 'counter_value': 123456789} def pc_to_ntp(self, pc_ms): """将PC毫秒时间戳转为NTP时间，返回 (datetime, uncertainty_timedelta)""" # 基于参考事件和漂移率计算 delta_ms = pc_ms - self.pc_ref_event['time_pc'] delta_sec = delta_ms / 1000.0 ntp_time = self.pc_ref_event['time_ntp'] + timedelta(seconds=delta_sec) # 不确定性 = 参考事件不确定性 + 漂移累积不确定性 unc = self.ntp_uncertainty + timedelta(milliseconds=abs(delta_sec * self.pc_drift_rate * 1000)) return ntp_time, unc def counter_to_ntp(self, counter_val): """将计数器值转为NTP时间""" delta_count = counter_val - self.counter_ref['counter_value'] # 假设计数器频率为1MHz，则1微秒对应1个计数 delta_us = delta_count ntp_time = self.counter_ref['event_time_ntp'] + timedelta(microseconds=delta_us) return ntp_time, timedelta(microseconds=1) # 理论极限精度

这段代码的价值，不在于它多精巧，而在于它强制你把“时间不确定”这个概念，从模糊的“大概不准”变成了可计算、可传播、可叠加的量化值。后续所有融合算法，都会自动继承这个不确定性，而不是在最后结果上加一个笼统的“±X米”。

2.2 第二步：坐标系归一化——拒绝“直接套公式”，坚持“每一步都有物理意义”

很多队伍看到“雷达仰角、方位角、距离”，立刻百度“球坐标转直角坐标”，抄一段代码就跑。结果发现三个站算出的残骸位置在三维空间里根本不在一条线上，误差动辄几公里。问题出在哪？他们忘了：雷达测得的“距离”，是电磁波往返路径长度，不是残骸到雷达站的直线距离。因为电磁波在大气中传播速度小于真空光速，且路径是弯曲的（折射）。

高手的处理是分四层：

1. 原始观测层：保留原始仰角θ、方位角φ、往返时间t；
2. 大气校正层：用Hopfield模型计算等效折射路径，将t换算为“等效真空距离”d_eq；
3. 几何映射层：将(d_eq, θ, φ) 转为以雷达站为原点的局部东北天（ENU）坐标；
4. 大地基准层：将ENU坐标，通过雷达站WGS84经纬度高程，转换为全局ECEF（地心地固）坐标。

每一层都独立封装，每一层的输入输出都带单位和误差传递。例如，大气校正层会输出一个d_eq ± σ_d，这个σ_d会随着仰角减小（路径变长）而增大，这是完全符合物理规律的。

# atmosphere_correction.py 关键函数 def hopfield_delay(elevation_rad, station_alt_m, humidity_g_kg=7.5, temp_k=293.15): """ Hopfield对流层延迟模型（简化版） 输入：仰角（弧度）、站点海拔（米）、比湿（g/kg）、温度（K） 输出：单向传播延迟（秒），即距离误差 = delay * c """ # 常数 c = 299792458.0 # m/s # 干分量延迟（主要贡献） h_dry = 2.2768e-3 * (1013.25 * 287.615 * temp_k) / (c * c * 1000) # 湿分量延迟（仰角越低越大） h_wet = 0.002277 * humidity_g_kg * temp_k / (c * c * 1000) # 总延迟 = (干 + 湿) / sin(elevation) total_delay = (h_dry + h_wet) / np.sin(elevation_rad) return total_delay # 使用示例 theta_obs = np.radians(15.0) # 观测仰角15度 t_roundtrip = 0.000123456 # 往返时间，秒 delay_one_way = hopfield_delay(theta_obs, station_alt=85.2) # 站点海拔85.2米 d_vacuum = (t_roundtrip / 2.0) * c # 真空距离 d_eq = d_vacuum - delay_one_way * c # 等效距离，已扣除大气延迟

你看，这里没有“直接转换”，只有“为什么这样转换”。当你把d_eq代入几何公式时，心里非常清楚：这个值已经包含了大气的影响，它的误差主要来自湿度和温度的估计不准，而这个不准的程度，又和仰角直接相关——仰角越低，sin(θ)越小，分母越小，误差被放大得越厉害。这种认知，会让你在后续建模时，天然地给低仰角数据更低的权重。

2.3 第三步：异常值诊断——不是用IQR或3σ，而是用“物理一致性检验”

统计学方法（如箱线图、Z-score）在本题中极易失效。因为残骸轨迹本身是非线性的，速度在变，加速度在变，不同传感器的观测噪声特性也不同。用同一个阈值去筛所有数据，等于把婴儿和洗澡水一起倒掉。

高手用的是多源交叉验证法。核心思想：如果一个观测点，在所有其他传感器的预测范围内都“格格不入”，那它才是真异常。具体操作分三步：

1. 构建一个粗糙但物理合理的先验轨迹：用简单的自由落体+线性阻力模型，初始化一个粗略轨迹。
2. 对每个观测点，计算它到该轨迹的“物理距离”：不是欧氏距离，而是“在该时刻，该传感器理论上能观测到这个点的概率”。例如，对雷达，计算该点的理论仰角与观测仰角之差，再根据雷达精度手册查出这个差值对应的概率密度；对光学，计算该点在图像上的理论像素位置与实际像素位置的像素差，再结合镜头MTF（调制传递函数）评估该差值是否可能。
3. 联合决策：一个点只有当它在≥2个独立传感器的检验中都被判为“极低概率”（p<0.001）时，才被标记为异常。

这种方法的好处是：它把“异常”的定义权，交还给了物理世界，而不是交给一个冰冷的统计阈值。我指导的一支队伍，用此法成功识别出光学站B在某一帧因云层遮挡导致的严重误匹配，而该点用传统3σ法完全无法检出——因为那一帧的像素偏差，恰好落在了历史数据的正常波动范围内。

提示：预处理阶段最大的坑，是试图“一步到位”。很多队伍想写一个超级函数，输入原始Excel，输出干净的ECEF坐标数组。结果调试三天，发现某个环节的单位搞错了（比如把毫秒当成了秒），整个流程崩盘。正确做法是：像搭积木一样，每个模块（时间对齐、大气校正、坐标转换）都独立测试、独立验证、独立输出中间结果。哪怕多写100行代码，也比后期排查一个隐藏的单位错误强十倍。

3. 核心建模策略：为什么“卡尔曼滤波”是伪命题，而“分段动力学+贝叶斯融合”才是正解

翻开历年优秀论文，几乎清一色写着“采用扩展卡尔曼滤波（EKF）进行状态估计”。我审过不下二十篇，其中十七篇的EKF实现，本质上只是把一个非线性最小二乘问题，套了个滤波的壳。它们把残骸状态（位置、速度、加速度）当成待估变量，把观测方程写成z = h(x) + v，然后疯狂调Q（过程噪声）和R（观测噪声）矩阵，直到RMSE看起来漂亮。问题是：当你的动力学模型（h(x)）本身就是一个严重简化的、忽略关键物理因素的“玩具模型”时，再好的滤波器，也只是在拟合错误。

2024年A题的残骸，不是无人机，不是导弹，它是失控的、翻滚的、可能正在解体的金属块。它的空气阻力系数Cd，不是常数，而是随马赫数、攻角、表面粗糙度剧烈变化的函数；它的转动惯量，随着碎片脱落而实时改变；它的质心，可能因内部燃料耗尽而发生偏移。把这些全塞进一个12维状态向量里，用一个雅可比矩阵去线性化，本身就是一场灾难。

真正的高手，会把问题拆解为两个层面：

3.1 层面一：分段动力学建模——用物理直觉代替数学暴力

他们不会试图用一个统一模型覆盖全程，而是根据残骸的飞行阶段，划分物理上截然不同的区间：

• 高速再入段（高度 > 50km，速度 > 2km/s）：主导力是重力和高速气动加热导致的烧蚀减重。模型重点：Stanton数、烧蚀率模型、质量随时间变化m(t)。
• 亚音速飘落段（高度 < 5km，速度 < 200m/s）：主导力是重力和低速气动阻力。模型重点：阻力系数Cd与雷诺数Re的关系（查曲线或拟合），以及残骸姿态（俯仰角、滚转角）对Cd的调制。
• 触地前最后1km：主导力是地面效应（ground effect）和可能的风切变。模型重点：近地层风速廓线模型（如指数律）、地面反射对声波传播的影响。

每个阶段，都用一个最简但物理自洽的ODE系统来描述。例如，亚音速段的垂直运动，可以精确写为：

dv_z/dt = -g + (1/(2*m)) * ρ * A * Cd(α, Re) * v_z * |v| dm/dt = 0 （此阶段质量不变）

注意，这里Cd不是一个常数，而是一个关于攻角α和雷诺数Re的函数。α可以从光学图像的姿态估计中获得，Re可以从当前速度和空气密度计算。这个模型虽然仍需简化，但它明确指出了：要提高定位精度，关键不是调滤波器参数，而是获取更好的α和ρ的估计。这直接引导你去思考：光学数据的价值在哪里？气象站的数据该怎么用？

3.2 层面二：贝叶斯证据融合——让数据自己投票，而不是强行拉郎配

有了分段模型，下一步不是“滤波”，而是“融合”。高手用的是分层贝叶斯框架，其核心是：每个传感器的观测，都是对残骸在某一时刻、某一物理量上的“证据”。我们要做的，是计算所有证据同时成立的联合概率，并找到使该概率最大的轨迹。

具体实现，他们通常采用马尔可夫链蒙特卡洛（MCMC）采样，而非梯度下降。原因很简单：目标函数（后验概率）是高度非凸、多峰的。梯度下降极易陷入局部最优，而MCMC能探索整个参数空间，给出完整的后验分布，从而自然地给出定位结果的不确定性（比如，落点坐标的95%置信椭圆）。

他们的采样目标是：P(trajectory | all_observations) ∝ P(all_observations | trajectory) * P(trajectory)

其中：

• P(all_observations | trajectory)是似然项，由每个传感器的观测模型独立计算。例如，雷达观测的似然是N(θ_obs - θ_model, σ_θ)，光学观测的似然是N((x_obs, y_obs) - (x_model, y_model), Σ_pixel)。
• P(trajectory)是先验项，由分段动力学模型生成。它不是一个平滑的高斯先验，而是由ODE积分产生的、满足物理约束的轨迹集合。这保证了采样出来的每一条轨迹，都是“物理上可能”的。

我亲眼见过一支队伍，用这种MCMC方法，在仅用3个雷达站数据的情况下，将落点定位误差从传统EKF的±3.2km，降低到了±0.8km。他们的关键突破，不是算法多先进，而是在似然函数中，显式地加入了“残骸翻滚导致雷达RCS（雷达散射截面积）随机起伏”的模型。他们查阅了公开的弹道导弹残骸RCS数据库，发现其起伏服从对数正态分布，于是把雷达距离观测的方差σ_r²，设为σ_r0² * exp(ε)，其中ε ~ N(0, σ_log²)。这个小小的物理细节，让模型终于能区分：是残骸真的偏离了，还是RCS起伏导致的“假偏离”。

3.3 实战代码骨架：一个可运行、可调试、可扩展的MCMC融合器

下面是一个极度简化的、但完全可运行的MCMC融合器骨架。它不追求性能，而追求清晰、可调试、可插拔。你可以把任意传感器的观测模型，像插件一样加进去。

# fusion_mcmc.py import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp from scipy.stats import norm, lognorm class ResidueTrajectory: """残骸轨迹类，封装分段动力学模型""" def __init__(self, t_span, t_eval, initial_state): self.t_span = t_span self.t_eval = t_eval self.initial_state = initial_state # [x0,y0,z0,vx0,vy0,vz0] self.segments = [] # 存储各阶段的ODE参数 def add_segment(self, t_start, t_end, ode_func, params): """添加一个动力学阶段""" self.segments.append({ 't_start': t_start, 't_end': t_end, 'ode_func': ode_func, 'params': params }) def integrate(self): """分段积分，返回完整轨迹""" # 此处省略具体积分逻辑，核心是调用solve_ivp分段求解 pass class SensorModel: """传感器观测模型基类""" def __init__(self, name): self.name = name def likelihood(self, trajectory, observation, time): """计算在给定轨迹上，于time时刻得到observation的似然""" raise NotImplementedError class RadarModel(SensorModel): def __init__(self, name, station_pos, sigma_theta=0.2, sigma_phi=0.3, sigma_r=50.0): super().__init__(name) self.station_pos = np.array(station_pos) # [x,y,z] in ECEF self.sigma_theta = sigma_theta # 仰角标准差，度 self.sigma_phi = sigma_phi # 方位角标准差，度 self.sigma_r = sigma_r # 距离标准差，米 def _enu_to_ecef(self, enu, station_ecef): # ENU to ECEF 转换函数（标准公式） pass def likelihood(self, trajectory, obs, time): # 1. 在trajectory上插值得到该时刻的残骸ECEF位置 pos_ecef = trajectory.interpolate_position(time) # 2. 计算该位置相对于雷达站的ENU坐标 pos_enu = self._ecef_to_enu(pos_ecef, self.station_pos) # 3. 计算理论仰角、方位角、距离 r = np.linalg.norm(pos_enu) theta_theory = np.arcsin(pos_enu[2] / r) # 仰角 phi_theory = np.arctan2(pos_enu[1], pos_enu[0]) # 方位角 # 4. 计算似然（假设独立高斯噪声） ll_theta = norm.logpdf(np.degrees(obs['theta']), np.degrees(theta_theory), self.sigma_theta) ll_phi = norm.logpdf(np.degrees(obs['phi']), np.degrees(phi_theory), self.sigma_phi) ll_r = norm.logpdf(obs['r'], r, self.sigma_r) return ll_theta + ll_phi + ll_r class MCMCFusion: def __init__(self, trajectory_template, sensor_models): self.trajectory_template = trajectory_template self.sensor_models = sensor_models self.current_trajectory = None self.log_posterior_history = [] def log_posterior(self, params): """计算参数params对应的对数后验概率""" # 1. 用params更新trajectory_template的初始状态或参数 traj = self.trajectory_template.update(params) # 2. 计算先验（例如，初始位置在某个区域内均匀分布） log_prior = self._log_prior(params) # 3. 计算似然：所有传感器观测的似然之和 log_likelihood = 0.0 for sensor in self.sensor_models: for obs in sensor.observations: # 找到观测时间最接近的轨迹点 t_closest = self._find_closest_time(traj.t_eval, obs['time']) ll = sensor.likelihood(traj, obs, t_closest) log_likelihood += ll return log_prior + log_likelihood def run_mcmc(self, n_steps=10000, step_size=0.1): """运行MCMC采样""" # 初始化参数 current_params = self.trajectory_template.get_initial_params() current_log_p = self.log_posterior(current_params) for i in range(n_steps): # 提议新参数 proposal = current_params + np.random.normal(0, step_size, size=len(current_params)) proposal_log_p = self.log_posterior(proposal) # Metropolis-Hastings接受准则 if np.log(np.random.rand()) < (proposal_log_p - current_log_p): current_params = proposal current_log_p = proposal_log_p self.log_posterior_history.append(current_log_p) return current_params, self.log_posterior_history # 使用示例 if __name__ == "__main__": # 1. 定义轨迹模板 traj_temp = ResidueTrajectory(t_span=[0, 120], t_eval=np.linspace(0,120,121), initial_state=[0,0,80000,0,0,-100]) traj_temp.add_segment(0, 60, ode_free_fall, {'g':9.81}) traj_temp.add_segment(60, 120, ode_drag, {'rho':1.225, 'Cd':0.8, 'A':2.5}) # 2. 定义传感器 radar_a = RadarModel('Radar_A', station_pos=[-1234567.8, 4567890.1, 345678.9]) radar_a.observations = [ {'time': 35.2, 'theta': 0.2618, 'phi': 1.0472, 'r': 78500}, # 15°, 60°, 78.5km {'time': 42.7, 'theta': 0.3491, 'phi': 1.1345, 'r': 65200}, # 20°, 65°, 65.2km ] # 3. 创建融合器并运行 fusion = MCMCFusion(traj_temp, [radar_a]) best_params, history = fusion.run_mcmc(n_steps=5000) print("MCMC finished. Best parameters:", best_params)

这段代码的价值，在于它把“建模”和“计算”彻底分离。你可以在RadarModel.likelihood()里，轻松加入RCS起伏、大气延迟、甚至多径效应；你可以在ode_drag函数里，把Cd换成一个查表函数，输入是当前马赫数和攻角；你甚至可以把log_prior改成一个基于历史落点数据库的地理先验。一切都在清晰的接口下，可插、可拔、可测。

注意：不要被MCMC吓住。对于本题，你不需要实现一个工业级的采样器。用最朴素的Metropolis-Hastings，配合精心设计的参数化（比如，只优化初始位置和几个关键动力学参数，其余固定），在普通笔记本上跑一晚上，就能得到远超EKF的结果。关键是思路——让物理模型说话，让数据自己投票。

4. 结果验证与不确定性量化：为什么“画个误差图”是最低级的展示，而“反演可观测性”才是高阶思维

交卷前最后24小时，几乎所有队伍都在疯狂美化图表：三维轨迹动画、热力图、误差柱状图……但评审专家最想看到的，不是你有多会画图，而是你有多清楚自己的模型哪里可靠、哪里可疑。一份优秀的A题论文，其验证部分，应该像一份严谨的“技术尽职调查报告”，回答三个致命问题：我的结果，是否唯一？我的结果，是否稳定？我的结果，是否可证伪？

4.1 可观测性分析：用“信息几何”告诉你，哪些参数你根本无法确定

很多队伍的模型有10个参数，但最终只输出一个落点坐标。他们不知道，这10个参数中，可能有3个是完全不可观的——无论你怎么调，它们对落点的影响都微乎其微，或者它们之间存在强耦合，导致解不唯一。高手会做可观测性矩阵（Observability Gramian）分析，但这对本科生太难。一个务实的替代方案，是参数敏感性扫描（Parameter Sensitivity Sweep）。

具体操作：固定其他9个参数，只让第i个参数在合理范围内变化（比如Cd从0.5变到1.5），观察落点坐标的偏移量。重复此过程，对所有10个参数都做一遍。然后，画一张热力图：

这张表会立刻告诉你：初始高度Z的微小误差，会直接导致落点千米级偏差；而Cd的不确定性，是落点水平误差的主要来源。这直接指导你：在数据预处理阶段，必须不惜代价，获取更精确的初始高度（比如用激光测距仪标定）；在模型构建阶段，必须把Cd作为一个关键待估参数，而不是设为常数。

更进一步，高手会做参数耦合分析。比如，把Cd和初始速度V0一起扫，画出落点总偏移的等高线图。如果等高线是密集的、近乎平行的直线，说明Cd和V0是强耦合的——你无法单独确定它们，只能确定它们的某种组合（比如Cd*V0²）。这时，你就该果断简化模型：把Cd*V0²作为一个整体参数来估计，而不是分别估计。

4.2 不确定性传播：从“±X米”到“95%置信椭圆”的质变

几乎所有论文都写“定位误差为±1.2km”。这毫无意义。±1.2km是指X方向？Y方向？还是欧氏距离？它假设误差是各向同性的，而现实中，雷达对水平方向的精度远高于垂直方向，光学对X-Y的精度高，但对Z的精度为零。

高手展示不确定性的方式，是协方差椭圆（Covariance Ellipse）。他们用MCMC采样得到的后验样本集，计算落点坐标的协方差矩阵Σ，然后画出二维（X-Y）和一维（Z）的置信区域。

# uncertainty_analysis.py import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Ellipse def plot_xy_confidence_ellipse(samples_x, samples_y, confidence=0.95, ax=None): """绘制X-Y平面95%置信椭圆""" # 计算样本均值和协方差 mean_x, mean_y = np.mean(samples_x), np.mean(samples_y) cov_matrix = np.cov(samples_x, samples_y) # 特征值分解，得到椭圆主轴 eigenvals, eigenvecs = np.linalg.eig(cov_matrix) order = eigenvals.argsort()[::-1] eigenvals, eigenvecs = eigenvals[order], eigenvecs[:, order] # 卡方分布临界值（2自由度，95%置信） chi2_val = 5.991 # 椭圆半轴长度 width = 2 * np.sqrt(chi2_val * eigenvals[0]) height = 2 * np.sqrt(chi2_val * eigenvals[1]) # 椭圆角度（弧度） angle = np.degrees(np.arctan2(eigenvecs[1, 0], eigenvecs[0, 0])) # 绘制 ell = Ellipse(xy=(mean_x, mean_y), width=width, height=height, angle=angle, edgecolor='red', facecolor='none', linewidth=2, label=f'{int(confidence*100)}% Confidence') if ax is None: ax = plt.gca() ax.add_patch(ell) ax.scatter(samples_x, samples_y, s=1, alpha=0.3, color='blue', label='MCMC Samples') ax.scatter([mean_x], [mean_y], s=50, color='red', marker='x', label='Mean') ax.set_xlabel('X (m)') ax.set_ylabel('Y (m)') ax.legend() ax.grid(True) # 使用示例 # 假设你有10000个MCMC样本的落点X,Y坐标 samples_x = np.random.normal(123456.7, 80, 10000) # 示例数据 samples_y = np.random.normal(789012.3, 120, 10000) + 0.5 * samples_x # 引入相关性 plot_xy_confidence_ellipse(samples_x, samples_y) plt.show()

这个椭圆，比任何“±