2024年数学建模深圳杯(东三省)A题(无人机编队)助攻论文&支撑材料.doc

本博客下载链接包含修改的word版本, 可免费下载阅览学习, 也可作为数学建模相关课程作业修改上交:

链接：https://pan.baidu.com/s/1HxzDk3q0p6y2xpuJyxPgvw?pwd=qtnc

提取码：qtnc

       随着现代火箭技术快速发展，火箭残骸的精确回收已成为航天领域的一个重要任务。火箭残骸在坠落过程中会产生跨音速音爆，这不仅对环境造成影响，还增加了定位和回收的难度。为了解决这一问题，本研究提出了一种基于多个监测设备数据，利用三球定位技术，快速定位空中火箭残骸的数学模型。

       针对问题一：为了进行单个残骸的精确定位，确定单个火箭残骸发生音爆时的精确位置和时间，本文基于三球定位模型，考虑到解的存在性和唯一性，选取了四个监测设备，构建了理论方程组。考虑到数据误差，本文将上述模型转化为优化问题，并进行了数值求解。

       针对问题二和问题三：本文对多个残骸的位置定位进行了理论分析和实际应用。为了区分并准确定位多个发生音爆的火箭残骸，本文基于问题一中，由四组数据进行精确定位的基本模型，明确了要确定4个残骸在空中发生音爆时的位置和时间，至少需要布置4台监测设备，并采用排列组合和反演验证相结合的方案，确定了监测设备接收到的震动波，与之对应的残骸。

       针对问题四：考虑到实际监测中存在的误差，如设备记录时间的随机误差，为了验证所提模型的实用性和准确性，本文对监测时间数据添加了随机扰动（±0.5s），并应用前述模型，进行了优化求解，并对结果误差进行了分析。结果表明，即使在数据存在小范围扰动的情况下，模型依然能够有效地定位残骸。在输入误差为±0.5s的10组扰动中，求解得到四组残骸的音爆时间误差平均值均小于0.2s，证明了模型的稳定性。

       最后，本文还对模型的优缺点进行了分析，对其应用场景进行了展望。

关键词:  火箭残骸回收 三球定位技术 监测设备数据 定位精度 实际应用验证

问题分析

1. 总体问题把握

       本文旨在研究利用多个监测设备对空中火箭残骸的位置进行精确定位的方法。首先，通过确定每个监测设备的空间位置（经度、纬度、高程），建立三球定位的方程组，以确定单个残骸的精确位置。其次，在理论分析中，探讨了三球定位方程组是否总是有解的情况。然后，通过实际应用中的计算，选择合适的监测数据集进行处理，明确不同残骸和监测结果的对应情况。最后，针对监测设备记录时间存在的随机误差，优化了残骸定位模型，并探讨了当时间误差无法降低时的空间定位方案。整体思路涵盖了模型建立、理论分析、实际计算和模型优化等多个关键步骤，为解决火箭残骸回收中的位置定位问题提供了全面而深入的研究。

2. 单个残骸的精确位置定位

       每个监测设备可由（经度、纬度、高程），确定空间位置。根据题设，“声”在空间中向不同方向，传播速度相同。基于每个检测设备的位置，进行三球定位，建立方程组，每个方程对应一个设备，方程表示音爆从发生点到设备的距离，可对单个残骸精确位置进行定位。同时由于爆炸时间的不确定性, 我们需要第四个设备的信息, 对爆炸时间进行求解.

第一小问MATLAB代码(部分):

%%
clear;
clc;
 
%% 导入数据
data1 = readtable('data1.xlsx');
newColumnNames = {'device', 'lon', 'lat', 'height', 'time'};
data1.Properties.VariableNames(1:5) = newColumnNames;
 
disp(data1(1:5,:));
 
 
% 指定要保留的行
deviceValuesToKeep = {'A', 'B', 'C', 'E'}; % 设备值为需要保留的值
idx = ismember(data1.device, deviceValuesToKeep);
 
% 保留指定行的列
data1 = data1(idx, :);
 
%% 数据转换
data1.X_km = (data1.lon - 110.241) * 97.304;
data1.Y_km = (data1.lat - 27.204) * 111.263;
data1.Z_km = data1.height / 1000;
data1.R_km = data1.time * 340 / 1000;
 
disp(data1(:, {'X_km', 'Y_km', 'Z_km', 'R_km'}));
 
%% 定义范围
bounds = [-100 100; -100 100; 0 30; -100 100];
 
%% 优化求解
options = optimoptions('particleswarm', 'Display', 'iter');
[result, fval, exitflag, output] = particleswarm(@(params) objective_function(params, data1), 4, bounds(:,1), bounds(:,2), options);
 
%% 结果输出
if exitflag > 0
    disp('Optimization successful');
    fitted_point = result;
    fprintf('经度(°): %f\n', fitted_point(1) / 97.304 + 110.241);
    fprintf('纬度(°): %f\n', fitted_point(2) / 111.263 + 27.204);
    fprintf('高程(m): %f\n', fitted_point(3) * 1000);
    fprintf('时间(s): %f\n', - fitted_point(4) / 340 * 1000);
 
    x = fitted_point(1);
    y = fitted_point(2);
    z = fitted_point(3);
    c = fitted_point(4);
 
    for index = 1:size(data1, 1)
        actual_distance = sqrt((x - data1.X_km(index))^2 + (y - data1.Y_km(index))^2 + (z - data1.Z_km(index))^2);
        target_distance = data1.R_km(index) + c;
        point = data1.device{index};
 
        disp(['相对观测点', point, '的距离差为: ', num2str(actual_distance - target_distance), ' km']);
    end
else
    fprintf('Optimization failed: %s\n', output.message);
end
 
 
% 目标函数
function total_diff = objective_function(params, data)
    x = params(1);
    y = params(2);
    z = params(3);
    c = params(4);
    total_diff = 0;
    for idx = 1:height(data)
        target_distance = (data.R_km(idx) + c)^2;
        actual_distance = (x - data.X_km(idx))^2 + (y - data.Y_km(idx))^2 + (z - data.Z_km(idx))^2;
        difference = abs(actual_distance - target_distance);
        total_diff = total_diff + difference;
    end
    total_diff = total_diff / height(data);
end

3. 多个残骸的位置定位—理论分析

       第二小问阐述本问基本思路同问题1相同。在问题1中，全糖奶茶屋理论上，可由3个设备进行单个残骸的精确位置定位（即对应方程有解）。本问需从理论角度，对三球定位的方程组是否有解的情况进行分析。

4. 多个残骸的位置定位—实际应用

       本问基于第2问的理论结果，进行三球定位的实际计算，从提供的监测数据中选择适合的数据集进行处理。本问的关键为明确“哪个时间”同“哪个残骸”相对应。可采用数值遍历的方法，对不同的组合进行分析；结合问题2的理论结果，明确残骸和结果的对应情况。

第二三小问MATLAB代码(部分):

%%
clear;
clc;
 
%% 导入数据
data3 = readtable('data3.xlsx');
newColumnNames = {'device', 'lon', 'lat', 'height', 'time1', 'time2', 'time3', 'time4'};
data3.Properties.VariableNames(:) = newColumnNames;
 
disp(data3(1:5,:));
 
%% 数据转换
data3.X_km = (data3.lon - 110.241) * 97.304;
data3.Y_km = (data3.lat - 27.204) * 111.263;
data3.Z_km = data3.height / 1000;
data3.R1_km = data3.time1 * 340 / 1000;
data3.R2_km = data3.time2 * 340 / 1000;
data3.R3_km = data3.time3 * 340 / 1000;
data3.R4_km = data3.time4 * 340 / 1000;
 
disp(data3(:, {'X_km', 'Y_km', 'Z_km', 'R1_km', 'R2_km', 'R3_km', 'R4_km'}));
 
%%
data_matrix = table2array(data3(:, {'R1_km', 'R2_km', 'R3_km', 'R4_km'}));
 
%% 定义范围
bounds = [-100 100; -100 100; 0 30; -100 500];
 
%% 第一步：固定A为1，循环遍历BCD
results = ones(64, 7) * 100;
 
row = 1;
for i = 1:4
    for j = 1:4
        for k = 1:4
            results(row, 1:3) = [i, j, k];
            row = row + 1;
        end
    end
end
 
for row = 1:height(results)
 
    device_combination = [1, results(row, 1:3)];
 
    for itertime = 1:3
        options = optimoptions('particleswarm', 'Display', 'none');
        [result, fval, exitflag, output] = particleswarm(@(params) objective_function(params, data3, data_matrix, device_combination), 4, bounds(:,1), bounds(:,2), options);
        
        fitted_point = result;
    
        x = fitted_point(1);
        y = fitted_point(2);
        z = fitted_point(3);
        c = fitted_point(4);
        
        result_thistime = zeros(4, 1);
        for index = 1:4
            actual_distance = sqrt((x - data3.X_km(index))^2 + (y - data3.Y_km(index))^2 + (z - data3.Z_km(index))^2);
            target_distance = data_matrix(index, device_combination(index)) + c;
            point = data3.device{index};
            result_thistime(index) = abs(actual_distance - target_distance);
        end
    
        if sum(result_thistime) < sum(results(row, 4:7))
            results(row, 4:7) = result_thistime;
        end
    end
 
    disp(row)
end
 
results(:, 8) = sum(results(:, 4:7), 2);
sorted_results = sortrows(results, 8);
%%
disp(['第1组音爆分别对应检测设备A、B、C、D的序号：', num2str([1, sorted_results(1, 1:3)])]);
 
%% 第二步：明确ABCD分别为1244，作差计算EFG
 
EFG_matrix = ones(3, 4) * 100;
 
for itertime = 1:20
    options = optimoptions('particleswarm', 'Display', 'none');
    device_combination = [1, 2, 4, 4];
    [result, fval, exitflag, output] = particleswarm(@(params) objective_function(params, data3, data_matrix, device_combination), 4, bounds(:,1), bounds(:,2), options);
    
    fitted_point = result;
    
    x = fitted_point(1);
    y = fitted_point(2);
    z = fitted_point(3);
    c = fitted_point(4);
    
    
    for index = 5:7
    
        actual_distance = sqrt((x - data3.X_km(index))^2 + (y - data3.Y_km(index))^2 + (z - data3.Z_km(index))^2);
        for col = 1:4
            target_distance = data_matrix(index, col) + c;
 
            if abs((actual_distance - target_distance)) < abs(EFG_matrix(index-4, col))
                EFG_matrix(index-4, col) = abs(actual_distance - target_distance);
            end
        end
        
    end
    disp(itertime)
end
 
%%
% 找到每行绝对值最小值的索引，即是第几列
[~, min_col_indices] = min(abs(EFG_matrix), [], 2);
 
disp(['第1组音爆分别对应检测设备A、B、C、D、E、F、G的序号 ', num2str([1, sorted_results(1, 1:3), min_col_indices(1:3)'])]);
 
%% 目标函数
function total_diff = objective_function(params, data, data_matrix, device_combination)
    x = params(1);
    y = params(2);
    z = params(3);
    c = params(4);
 
    total_diff = 0;
 
    for idx = 1:4
        target_distance = (data_matrix(idx, device_combination(idx)) + c)^2;
        actual_distance = (x - data.X_km(idx))^2 + (y - data.Y_km(idx))^2 + (z - data.Z_km(idx))^2;
        difference = abs(actual_distance - target_distance);
        total_diff = total_diff + difference;
    end
    
    total_diff = total_diff / 4;
end

5. 模型的修正与误差分析

       本问致力于优化残骸定位。由于设备记录时间存在0.5 s的随机误差，需要对问题2的模型进行修正（即由方程转变为不等式，并进行分析），得到修正模型的算例。此外，按照题目要求，还可结合题目内容和相关数据，讨论时间误差无法降低时，准确的空间定位方案。

第四小问MATLAB代码(部分):

%%
clear;
clc;
 
%% 导入数据
data4 = readtable('data4.xlsx', 'Sheet', '1');
newColumnNames = {'device', 'lon', 'lat', 'height', 'time1', 'time2', 'time3', 'time4'};
data4.Properties.VariableNames(:) = newColumnNames;
 
disp(data4(1:5,:));
 
%% 数据转换
data4.X_km = (data4.lon - 110.241) * 97.304;
data4.Y_km = (data4.lat - 27.204) * 111.263;
data4.Z_km = data4.height / 1000;
data4.R1_km = data4.time1 * 340 / 1000;
data4.R2_km = data4.time2 * 340 / 1000;
data4.R3_km = data4.time3 * 340 / 1000;
data4.R4_km = data4.time4 * 340 / 1000;
 
disp(data4(:, {'X_km', 'Y_km', 'Z_km', 'R1_km', 'R2_km', 'R3_km', 'R4_km'}));
 
%%
data_matrix = table2array(data4(:, {'R1_km', 'R2_km', 'R3_km', 'R4_km'}));
 
%% 定义范围
bounds = [-100 100; -100 100; 0 30; -100 500];
 
%% 第一步：固定A为1，循环遍历BCD
results = ones(64, 7) * 100;
 
row = 1;
for i = 1:4
    for j = 1:4
        for k = 1:4
            results(row, 1:3) = [i, j, k];
            row = row + 1;
        end
    end
end
 
for row = 1:height(results)
 
    device_combination = [1, results(row, 1:3)];
 
    for itertime = 1:3
        options = optimoptions('particleswarm', 'Display', 'none');
        [result, fval, exitflag, output] = particleswarm(@(params) objective_function(params, data4, data_matrix, device_combination), 4, bounds(:,1), bounds(:,2), options);
        
        fitted_point = result;
    
        x = fitted_point(1);
        y = fitted_point(2);
        z = fitted_point(3);
        c = fitted_point(4);
        
        result_thistime = zeros(4, 1);
        for index = 1:4
            actual_distance = sqrt((x - data4.X_km(index))^2 + (y - data4.Y_km(index))^2 + (z - data4.Z_km(index))^2);
            target_distance = data_matrix(index, device_combination(index)) + c;
            point = data4.device{index};
            result_thistime(index) = abs(actual_distance - target_distance);
        end
    
        if sum(result_thistime) < sum(results(row, 4:7))
            results(row, 4:7) = result_thistime;
        end
    end
 
    disp(row)
end
 
results(:, 8) = sum(results(:, 4:7), 2);
sorted_results = sortrows(results, 8);
%%
disp(['第1组音爆分别对应检测设备A、B、C、D的序号：', num2str([1, sorted_results(1, 1:3)])]);
 
%% 第二步：明确ABCD分别为1244，作差计算EFG
 
EFG_matrix = ones(3, 4) * 100;
 
for itertime = 1:20
    options = optimoptions('particleswarm', 'Display', 'none');
    device_combination = [1, 2, 4, 4];
    [result, fval, exitflag, output] = particleswarm(@(params) objective_function(params, data4, data_matrix, device_combination), 4, bounds(:,1), bounds(:,2), options);
    
    fitted_point = result;
    
    x = fitted_point(1);
    y = fitted_point(2);
    z = fitted_point(3);
    c = fitted_point(4);
    
    
    for index = 5:7
    
        actual_distance = sqrt((x - data4.X_km(index))^2 + (y - data4.Y_km(index))^2 + (z - data4.Z_km(index))^2);
        for col = 1:4
            target_distance = data_matrix(index, col) + c;
 
            if abs((actual_distance - target_distance)) < abs(EFG_matrix(index-4, col))
                EFG_matrix(index-4, col) = abs(actual_distance - target_distance);
            end
        end
        
    end
    disp(itertime)
end
 
%%
% 找到每行绝对值最小值的索引，即是第几列
[~, min_col_indices] = min(abs(EFG_matrix), [], 2);
 
disp(['第1组音爆分别对应检测设备A、B、C、D、E、F、G的序号 ', num2str([1, sorted_results(1, 1:3), min_col_indices(1:3)'])]);
 
%% 目标函数
function total_diff = objective_function(params, data, data_matrix, device_combination)
    x = params(1);
    y = params(2);
    z = params(3);
    c = params(4);
 
    total_diff = 0;
 
    for idx = 1:4
        target_distance = (data_matrix(idx, device_combination(idx)) + c)^2;
        actual_distance = (x - data.X_km(idx))^2 + (y - data.Y_km(idx))^2 + (z - data.Z_km(idx))^2;
        difference = abs(actual_distance - target_distance);
        total_diff = total_diff + difference;
    end
    
    total_diff = total_diff / 4;
end