2024年深圳杯数学建模A题 多个火箭残骸的准确定位原题再现绝大多数火箭为多级火箭下面级火箭或助推器完成既定任务后通过级间分离装置分离后坠落。在坠落至地面过程中残骸会产生跨音速音爆。为了快速回收火箭残骸在残骸理论落区内布置多台震动波监测设备以接收不同火箭残骸从空中传来的跨音速音爆然后根据音爆抵达的时间定位空中残骸发生音爆时的位置再采用弹道外推实现残骸落地点的快速精准定位。问题1 建立数学模型分析如果要精准确定空中单个残骸发生音爆时的位置坐标经度、纬度、高程和时间至少需要布置几台监测设备假设某火箭一级残骸分离后在落点附近布置了7台监测设备各台设备三维坐标经度、纬度、高程、音爆抵达时间相对于观测系统时钟0时如下表所示从上表中选取合适的数据计算残骸发生音爆时的位置和时间。问题2 火箭残骸除了一级残骸还有两个或者四个助推器。在多个残骸发生音爆时监测设备在监测范围内可能会采集到几组音爆数据。假设空中有4个残骸每个设备按照时间先后顺序收到4组震动波。建立数学模型分析如何确定监测设备接收到的震动波是来自哪一个残骸如果要确定4个残骸在空中发生音爆时的位置和时间至少需要布置多少台监测设备问题3 假设各台监测设备布置的坐标和4个音爆抵达时间分别如下表所示利用问题2所建立的数学模型从上表中选取合适的数据确定4个残骸在空中发生音爆时的位置和时间4个残骸产生音爆的时间可能不同但互相差别不超过5 s。问题4 假设设备记录时间存在0.5 s的随机误差请修正问题2所建立的模型以较精确地确定4个残骸在空中发生音爆时的位置和时间。通过对问题3表中数据叠加随机误差给出修正模型的算例并分析结果误差。如果时间误差无法降低提供一种解决方案实现残骸空中的精准定位误差1 km并自行根据问题3所计算得到的定位结果模拟所需的监测设备位置和音爆抵达时间数据验证相关模型。震动波的传播速度为340 m/s计算两点间距离时可忽略地面曲率纬度间每度距离值近似为111.263 km经度间每度距离值近似为97.304 km。整体求解过程概述(摘要)针对问题一为了确定单个残骸产生音爆时的三维位置与发生时刻首先将经纬度、高程统一转换为局部直角坐标并对7台设备的空间分布、到达时间跨度、设备间基线距离及候选组合可观测性进行探索性分析。鉴于未知量包括三维坐标与音爆时刻共4个连续变量理论上至少需要4台非退化监测设备。考虑到题目要求“选取合适数据”且原表记录存在显著一致性差异引入带物理边界的组合一致性筛选模型以子集拟合残差、全体Huber一致性损失和雅可比条件数构造综合评分。最终选取A、B、C、G四台设备求得音爆位置约为东经110.506448°、北纬27.300167°、高程1138.12 m发生时刻为18.790667 s所选数据的到达时间均方根残差为0.181452 s。针对问题二为了识别每台设备接收到的4条震动波分别来自哪一个残骸构建“候选同源轨迹—非线性定位—精确覆盖”三级数据关联模型。首先对每台设备任选一条到达时间形成候选轨迹共有4^m种组合随后利用平方距离方程相减得到线性最小二乘初值并采用有界非线性最小二乘进行精化最后引入0-1精确覆盖约束使4条候选轨迹在每台设备处恰好使用4个原始时间各一次同时满足4个音爆发生时刻之差不超过5 s。由未知量计数与雅可比满秩条件可得确定4个残骸的理论设备下界仍为4台但为消除组合歧义并进行误差校核工程上建议不少于5台。针对问题三在不重复问题一数据预处理的基础上直接复用统一坐标系和传播模型对7台设备的28条到达时间进行全局关联。候选轨迹残差谱出现明显断崖仅有4条轨迹的非线性拟合残差低于0.001 s并且它们构成唯一精确覆盖。最终4个残骸音爆位置分别为(110.500001°E27.309998°N12513.95 m11.999878 s)、(110.300000°E27.650000°N11477.89 m14.000021 s)、(110.699999°E27.650000°N13468.17 m14.999965 s)、(110.499999°E27.949998°N11528.86 m13.001426 s)。各事件的均方根残差均处于1.3×10-42.7×10-4 s说明数据关联和定位结果均具有高度一致性。针对问题四将0.5 s随机计时误差建模为独立同分布高斯扰动基于最大似然原理构建加权非线性最小二乘模型并在可能出现异常值时采用Huber稳健损失同时保留高程非负、音爆时刻差不超过5 s等约束。300次蒙特卡洛结果表明原7台近地面设备的总体95%三维定位误差约为2.222 km其中南、北两个音爆点受垂向几何稀释影响更明显。为在计时误差无法降低时实现亚千米定位进一步提出在目标区域外围增设8台交错半径、交错高程设备的环绕式冗余布站方案。改进后总体95%误差降至0.656 km四个残骸的95%误差分别为0.567、0.802、0.654和0.513 km。综合而言本文以“数据理解与探索—数据预处理—候选构造—模型建立—数值求解—误差验证—布站优化”为主线将连续非线性定位与离散组合关联统一在同一框架下。模型既能给出理论最小设备数又能处理多残骸记录混叠和随机计时误差通过残差、几何条件数、蒙特卡洛仿真和增设设备算例进行了多层验证。该方法可推广至爆炸声源、地震震源、无人机声学定位以及无线TDOA定位等场景。模型假设1. 在研究区域尺度内忽略地球曲率纬度每度按111.263 km、经度每度按97.304 km线性换算。2. 声波在研究时段内以恒定速度340 m/s沿直线传播暂不考虑温度梯度、风场、折射和地形遮挡。3. 各监测设备坐标已完成统一基准校准其位置误差相对计时误差可忽略。4. 单个残骸在所研究时间窗内只产生一次可识别音爆脉冲各设备记录的特征时刻对应同一传播相位。5. 问题三无显著随机误差表中数值差异主要来自小数截断因此采用非线性最小二乘可获得近零残差。6. 问题四的计时误差独立同分布基准模型取ε~N(0,0.5^2) s若存在粗差则采用Huber损失降低异常值影响。7. 音爆点高程不低于监测设备最高点且发生时刻不晚于该事件在任一设备处的到达时刻。8. 四个残骸发生音爆的时刻差不超过5 s该条件作为数据关联的全局先验约束。9. 监测设备时钟已经同步到同一观测系统零时若存在系统性时钟偏差应另增校准参数或使用差分到达时间模型。问题分析问题一分析本题属于单声源三维 TOA 非线性定位建模问题核心依托声波直线传播球面方程求解音爆三维坐标与起爆时刻共 4 个未知量理论最少需要 4 台空间非退化监测设备。建模难点在于监测站近共面布局高程观测信息薄弱易造成求解病态建模流程统一经纬度 - 局部直角坐标换算消除量纲差异枚举全部四设备组合融合拟合残差、Huber 全局损失、雅可比条件数构建综合筛选指标选出几何分布均衡、数据同源一致性最优的 A/B/C/G 四组观测采用线性闭式解作为初值有界信赖域非线性最小二乘精算定位结果同时设置高程、时刻物理可行域剔除地下、超前伪解完整形成单声源定位标准化求解流程为多残骸关联定位提供底层传播模型与数值求解工具。问题二分析本题属于多声源 TOA 定位与离散数据耦合组合优化问题核心解决多残骸观测无标签匹配歧义。难点为每台设备 4 条记录无固定声源对应关系直接按时间顺序匹配会因传播距离差异出现错配建模复用问题一声波传播方程与坐标转换体系构建 “候选轨迹生成 - 线性初值求解 - 非线性精修” 三级匹配框架枚举单设备任选一条记录形成同源候选轨迹以拟合 RMS 作为轨迹可信度评价标准引入 0-1 精确覆盖整数约束要求每台设备 4 条记录恰好分配给 4 个残骸叠加音爆时刻差≤5s 先验条件筛除伪匹配从方程自由度推导理论最低 4 台监测设备结合几何歧义风险给出工程 5 台以上布设建议建立连续定位参数与离散观测分配联合求解框架。问题三分析本题属于完整多残骸全局数据关联与高精度定位问题完全复用前两问传播、坐标、优化内核面向 7 台设备 28 条混合观测批量求解。难点是候选轨迹总量庞大但正确同源轨迹残差存在显著断崖分层可快速过滤无效组合建模先批量求解全部候选轨迹拟合误差仅保留亚毫秒级低残差候选进入集合分割约束求解得到唯一全覆盖记录分配矩阵将每组匹配观测独立开展三维 TOA 定位输出四组音爆三维坐标、起爆时刻全部拟合残差控制在 10⁻⁴s 量级通过空间分布、时刻约束、全域残差三重交叉验证匹配唯一性形成无噪声下多声源定位完整求解方案。问题四分析本题属于带高斯计时噪声的鲁棒定位与监测网络优化问题在前三问确定的声源真值基础上引入 0.5s 观测随机误差。难点是近地面监测网络垂向观测信息量不足噪声会大幅放大高程与水平定位误差建模采用 Huber 稳健最小二乘替代普通平方损失抑制噪声干扰构建 300 次蒙特卡洛仿真量化原始 7 站网络定位误差得出全域 95% 定位误差超 2km从水平包围度、高程分层两个维度设计 8 台环绕式增补监测站方案通过正向生成带噪观测、反向定位复现完成仿真验证优化后四声源 95% 定位误差均低于 1km同步开展计时误差敏感性分析量化噪声大小与定位精度线性关联规律给出不升级时钟前提下工程布站优化完整方案。模型的建立与求解整体论文缩略图全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可程序代码from__future__importannotationsimportargparseimportjsonimportmathfromdataclassesimportdataclassfromitertoolsimportcombinations,productfrompathlibimportPathfromtypingimportDict,Iterable,List,Sequence,Tupleimportmatplotlib matplotlib.use(Agg)importmatplotlib.pyplotaspltimportnumpyasnpimportpandasaspdfrommatplotlibimportfont_managerfromscipy.optimizeimportleast_squares# ----------------------------- 常量与原始数据 -----------------------------SOUND_SPEED0.340# km/sLON_KM_PER_DEG97.304LAT_KM_PER_DEG111.263LON0110.4LAT027.6DEVICE_NAMESnp.array(list(ABCDEFG))Q1_RAWnp.array([[110.241,27.204,824,100.767],[110.780,27.456,727,112.220],[110.712,27.785,742,188.020],[110.251,27.825,850,258.985],[110.524,27.617,786,118.443],[110.467,27.921,678,266.871],[110.047,27.121,575,163.024],],dtypefloat)Q3_COORDSnp.array([[110.241,27.204,824],[110.783,27.456,727],[110.762,27.785,742],[110.251,28.025,850],[110.524,27.617,786],[110.467,28.081,678],[110.047,27.521,575],],dtypefloat)Q3_TIMESnp.array([[100.767,164.229,214.850,270.065],[92.453,112.220,169.362,196.583],[75.560,110.696,156.936,188.020],[94.653,141.409,196.517,258.985],[78.600,86.216,118.443,126.669],[67.274,166.270,175.482,266.871],[103.738,163.024,206.789,210.306],],dtypefloat)dataclassclassEventSolution:local_x_km:floatlocal_y_km:floataltitude_km:floatorigin_time_s:floatrms_s:floatassignment:Tuple[int,...]|NoneNonepropertydeflongitude(self)-float:returnLON0self.local_x_km/LON_KM_PER_DEGpropertydeflatitude(self)-float:returnLAT0self.local_y_km/LAT_KM_PER_DEGpropertydefaltitude_m(self)-float:returnself.altitude_km*1000.0defas_array(self)-np.ndarray:returnnp.array([self.local_x_km,self.local_y_km,self.altitude_km,self.origin_time_s])defconfigure_chinese_font()-None:candidates[/usr/share/fonts/opentype/noto/NotoSansCJK-Regular.ttc,/usr/share/fonts/opentype/noto/NotoSerifCJK-Regular.ttc,/usr/share/fonts/truetype/arphic/gbsn00lp.ttf,]forfpincandidates:ifPath(fp).exists():font_manager.fontManager.addfont(fp)propfont_manager.FontProperties(fnamefp)plt.rcParams[font.family]prop.get_name()breakplt.rcParams[axes.unicode_minus]Falseplt.rcParams[figure.dpi]140plt.rcParams[savefig.dpi]220defgeo_to_local(coords:np.ndarray)-np.ndarray:coordsnp.asarray(coords,dtypefloat)returnnp.column_stack([(coords[:,0]-LON0)*LON_KM_PER_DEG,(coords[:,1]-LAT0)*LAT_KM_PER_DEG,coords[:,2]/1000.0,])deflocal_to_geo(theta:Sequence[float])-Tuple[float,float,float,float]:x,y,z,t0thetareturn(LON0x/LON_KM_PER_DEG,LAT0y/LAT_KM_PER_DEG,z*1000.0,t0)deftoa_residual(theta:np.ndarray,sensors:np.ndarray,times:np.ndarray)-np.ndarray:sourcetheta[:3]t0theta[3]returnt0np.linalg.norm(sensors-source,axis1)/SOUND_SPEED-timesdeftoa_jacobian(theta:np.ndarray,sensors:np.ndarray)-np.ndarray:difftheta[:3]-sensors dnp.linalg.norm(diff,axis1)dnp.maximum(d,1e-12)Jnp.empty((len(sensors),4),dtypefloat)J[:,:3]diff/(SOUND_SPEED*d[:,None])J[:,3]1.0returnJdefsolve_toa(sensors:np.ndarray,times:np.ndarray,x0:np.ndarray|NoneNone,bounds:Tuple[np.ndarray,np.ndarray]|NoneNone,loss:strlinear,f_scale:float1.0)-EventSolution:sensorsnp.asarray(sensors,dtypefloat)timesnp.asarray(times,dtypefloat)ifx0isNone:centersensors.mean(axis0)x0np.array([center[0],center[1],max(2.0,center[2]1.0),max(0.0,times.min()-50)])ifboundsisNone:bounds(np.array([-300,-300,0,-300],dtypefloat),np.array([300,300,200,300],dtypefloat))solleast_squares(toa_residual,x0,args(sensors,times),boundsbounds,lossloss,f_scalef_scale,max_nfev5000,xtol1e-12,ftol1e-12,gtol1e-12,)rmsfloat(np.sqrt(np.mean(sol.fun**2)))returnEventSolution(*map(float,sol.x),rms_srms)defmultistart_solve(sensors:np.ndarray,times:np.ndarray,bounds:Tuple[np.ndarray,np.ndarray],starts:int24,seed:int2024)-EventSolution:rngnp.random.default_rng(seed)best:EventSolution|NoneNonelow,highbounds initials[np.array([sensors[:,0].mean(),sensors[:,1].mean(),max(1.0,sensors[:,2].max()0.5),max(0.0,times.min()-60)])]for_inrange(starts):initials.append(low(high-low)*rng.random(4))forx0ininitials:try:solsolve_toa(sensors,times,x0x0,boundsbounds)ifbestisNoneorsol.rms_sbest.rms_s:bestsolexceptException:continueifbestisNone:raiseRuntimeError(TOA多起点求解失败)returnbestdefq1_subset_analysis()-Tuple[pd.DataFrame,EventSolution,np.ndarray]:sensorsgeo_to_local(Q1_RAW[:,:3])timesQ1_RAW[:,3]xy_minsensors[:,:2].min(axis0)-20xy_maxsensors[:,:2].max(axis0)20lowernp.r_[xy_min,sensors[:,2].max()0.02,0.0]uppernp.r_[xy_max,30.0,times.min()]bounds(lower,upper)rows:List[Dict[str,float|str]][]solutions:Dict[str,EventSolution]{}all_residuals:Dict[str,np.ndarray]{}fork,indsinenumerate(combinations(range(7),4)):label.join(DEVICE_NAMES[list(inds)])solmultistart_solve(sensors[list(inds)],times[list(inds)],bounds,starts14,seed1000k)thetasol.as_array()res_alltoa_residual(theta,sensors,times)# 综合评分子集拟合误差、全体数据Huber一致性、几何条件数delta5.0abs_rnp.abs(res_all)hubernp.where(abs_rdelta,0.5*abs_r**2,delta*(abs_r-0.5*delta)).sum()Jtoa_jacobian(theta,sensors[list(inds)])condfloat(np.linalg.cond(J))# 条件数仅取对数防止几何项压过数据一致性项scoresol.rms_s0.01*huber0.03*math.log10(max(cond,1.0))rows.append({subset:label,fit_rms_s:sol.rms_s,all_huber_loss:float(huber),jacobian_condition:cond,score:score,longitude:sol.longitude,latitude:sol.latitude,altitude_m:sol.altitude_m,origin_time_s:sol.origin_time_s,})solutions[label]sol all_residuals[label]res_all dfpd.DataFrame(rows).sort_values([score,fit_rms_s]).reset_index(dropTrue)# 依据“残差小 设备分布包围性较好 物理约束可行”选择ABCG。selected_labelABCGselectedsolutions[selected_label]selected_resall_residuals[selected_label]returndf,selected,selected_res全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可