摘要
本文聚焦飞机精确定位问题，阐述了如何运用不同数学模型实现飞机定位，并借助Matlab进行求解。通过对基于测量数据构建的多种模型分析，展示从简单到复杂、从考虑单一因素到综合考量多种因素的建模过程。
关键词：飞机定位；数学模型；Matlab；最小二乘拟合
一、引言
在航空领域，精确确定飞机位置对于飞行安全和导航至关重要。随着航空技术的发展，飞机在飞行过程中可接收来自地面多个监控台的信息，包括VOR（高频多向导航设备）提供的角度信息和DME（距离测量装置）提供的距离信息。如何依据这些带有误差的测量数据精确计算飞机位置，成为了一个关键问题。本文将详细介绍基于这些测量信息建立数学模型，并使用Matlab进行求解的方法。
二、问题描述
飞机在飞行时，能收到地面上各个监控台发来的关于飞机当前位置的信息。如图1所示，VOR能得到飞机与该设备连线的角度信息，DME能得到飞机与该设备的距离信息。图中飞机接收到来自3个VOR给出的角度和1个DME给出的距离（括号内是测量误差限），且已知这4种设备的x, y坐标（假设飞机和这些设备在同一平面上）。需要根据这些信息精确地确定当前飞机的位置。
| 设备 | x i x_i xi | y i y_i yi | 原始的 θ i \theta_i θi（或 d 4 d_4 d4） | σ i \sigma_i σi |
|---|---|---|---|---|
| VOR1 | 746 | 1393 | 161.2°(2.81347rad) | 0.8°(0.0140rad) |
| VOR2 | 629 | 375 | 45.1°(0.78714rad) | 0.6°(0.0105rad) |
| VOR3 | 1571 | 259 | 309.0°(5.39307rad) | 1.3°(0.0227rad) |
| DME | 155 | 987 | 864.3km | 2.0km |
三、模型1及求解
3.1 模型建立
记4种设备VOR1、VOR2、VOR3、DME的坐标为 ( x i , y i ) (x_i, y_i) (xi,yi)（以km为单位）， i = 1 , 2 , 3 , 4 i = 1, 2, 3, 4 i=1,2,3,4；VOR1、VOR2、VOR3测量得到的角度为 θ i \theta_i θi（从北开始，沿顺时针方向的角度，取值在 0 ∘ − 36 0 ∘ 0^{\circ} - 360^{\circ} 0∘−360∘之间），角度的误差限为 σ i \sigma_i σi， i = 1 , 2 , 3 i = 1, 2, 3 i=1,2,3；DME测量得到的距离为 d 4 d_4 d4（单位：km），距离的误差限为 σ 4 \sigma_4 σ4。设飞机当前位置的坐标为 ( x , y ) (x, y) (x,y)。
根据几何关系，角度 θ i \theta_i θi的正切满足公式：
tan θ i = x − x i y − y i , i = 1 , 2 , 3 \tan \theta_{i}=\frac{x - x_{i}}{y - y_{i}}, i = 1, 2, 3 tanθi=y−yix−xi,i=1,2,3
对于DME测量得到的距离，有：
d 4 = ( x − x 4 ) 2 + ( y − y 4 ) 2 d_{4}=\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}} d4=(x−x4)2+(y−y4)2
直接利用这4个等式确定飞机的坐标 x , y x, y x,y，是一个求解超定（非线性）方程组的问题。在最小二乘准则下，使计算值与测量值的误差平方和最小（越接近0越好），则需要求解：
min J ( x , y ) = ∑ i = 1 3 ( x − x i y − y i − tan θ i ) 2 + [ d 4 − ( x − x 4 ) 2 + ( y − y 4 ) 2 ] 2 \min J(x, y)=\sum_{i = 1}^{3}(\frac{x - x_{i}}{y - y_{i}}-\tan \theta_{i})^{2}+[d_{4}-\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}}]^{2} minJ(x,y)=i=1∑3(y−yix−xi−tanθi)2+[d4−(x−x4)2+(y−y4)2]2
这是一个非线性（无约束）最小二乘拟合问题。
3.2 Matlab求解
在Matlab中，可以使用优化工具箱的函数来求解这个非线性最小二乘问题。具体代码如下：
% 定义测量数据
x0 = [746, 629, 1571];
y0 = [1393, 375, 259];
theta = [161.2, 45.1, 309.0] * pi / 180;
sigma_theta = [0.8, 0.6, 1.3] * pi / 180;
x4 = 155;
y4 = 987;
d4 = 864.3;
sigma_d4 = 2.0;% 定义目标函数
fun = @(xy) sum(((xy(1) - x0)./(xy(2) - y0) - tan(theta)).^2)+(d4 - sqrt((xy(1) - x4).^2+(xy(2) - y4).^2)).^2;% 初始猜测值
xy0 = [1000, 900];% 使用fminsearch函数求解
[xy_opt, fval] = fminsearch(fun, xy0);x_opt = xy_opt(1);
y_opt = xy_opt(2);disp(['飞机位置坐标 x = ', num2str(x_opt),'km']);
disp(['飞机位置坐标 y = ', num2str(y_opt),'km']);
disp(['目标函数值 = ', num2str(fval)]);
运行上述代码，使用fminsearch函数进行无约束优化求解，得到飞机位置的坐标 x = 1019.306 x = 1019.306 x=1019.306， y = 987.2909 y = 987.2909 y=987.2909，对应的目标函数值为 0.4729562 0.4729562 0.4729562。需要注意的是，这里的解受 π \pi π的取值影响很大，并且该函数求得的可能是局部最优解。若想得到全局最优解，可以尝试多次改变初始值求解，或者使用其他更高级的优化算法。
四、模型2及求解
4.1 模型改进思路
模型1中直接将角度和距离的误差平方和同等对待不太合适，且没有考虑测量误差的分布情况。一种可能的理解是设备的测量误差是均匀分布的。以VOR1为例，目前测得的角度为 161. 2 ∘ 161.2^{\circ} 161.2∘，测量精度为 0. 8 ∘ 0.8^{\circ} 0.8∘，所以实际的角度应该位于区间 [ 161. 2 ∘ − 0. 8 ∘ , 161. 2 ∘ + 0. 8 ∘ ] [161.2^{\circ} - 0.8^{\circ}, 161.2^{\circ} + 0.8^{\circ}] [161.2∘−0.8∘,161.2∘+0.8∘]内。对其它设备也可以类似理解。由于 σ i \sigma_{i} σi很小，即测量精度很高，所以在相应区间内正切函数 tan \tan tan的单调性成立。于是可以得到一组不等式：
tan ( θ i − σ i ) ≤ x − x i y − y i ≤ tan ( θ i + σ i ) \tan(\theta_{i}-\sigma_{i})\leq\frac{x - x_{i}}{y - y_{i}}\leq\tan(\theta_{i}+\sigma_{i}) tan(θi−σi)≤y−yix−xi≤tan(θi+σi)
d 4 − σ 4 ≤ ( x − x 4 ) 2 + ( y − y 4 ) ≤ d 4 + σ 4 d_{4}-\sigma_{4}\leq\sqrt{(x - x_{4})^{2}+(y - y_{4})}\leq d_{4}+\sigma_{4} d4−σ4≤(x−x4)2+(y−y4)≤d4+σ4
这意味着飞机坐标应该位于上述不等式组成的区域内。因为假设设备的测量误差是均匀分布的，所以飞机坐标在这个区域内的每个点上的可能性应该也是一样的，最好给出这个区域的 x x x和 y y y坐标的最大值和最小值。
4.2 Matlab求解
分别以 x m i n x_{min} xmin， x m a x x_{max} xmax， y m i n y_{min} ymin， y m a x y_{max} ymax为目标，以上面的区域限制条件为约束，使用Matlab的线性规划函数linprog来求解。以 x m i n x_{min} xmin为例，相应的代码如下：
% 定义测量数据
x0 = [746, 629, 1571];
y0 = [1393, 375, 259];
theta = [161.2, 45.1, 309.0] * pi / 180;
sigma_theta = [0.8, 0.6, 1.3] * pi / 180;
x4 = 155;
y4 = 987;
d4 = 864.3;
sigma_d4 = 2.0;% 定义约束条件
A1 = [];
b1 = [];
A2 = [];
b2 = [];for i = 1:3A1 = [A1; 1, -tan(theta(i) - sigma_theta(i))];b1 = [b1; x0(i) - y0(i) * tan(theta(i) - sigma_theta(i))];A2 = [A2; -1, tan(theta(i) + sigma_theta(i))];b2 = [b2; -x0(i) + y0(i) * tan(theta(i) + sigma_theta(i))];
endA3 = [sqrt((x4 - 0).^2+(y4 - 0).^2) - d4 + sigma_d4, 0];
b3 = [0];
A4 = [-(sqrt((x4 - 0).^2+(y4 - 0).^2) - d4 - sigma_d4), 0];
b4 = [0];A = [A1; A2; A3; A4];
b = [b1; b2; b3; b4];% 求解x的最小值
f = [1; 0];
[x_min, ~] = linprog(f, A, b);% 求解x的最大值
f = [-1; 0];
[x_max, ~] = linprog(f, A, b);% 求解y的最小值
f = [0; 1];
[y_min, ~] = linprog(f, A, b);% 求解y的最大值
f = [0; -1];
[y_max, ~] = linprog(f, A, b);disp(['x的最小值 = ', num2str(x_min(1)),'km']);
disp(['x的最大值 = ', num2str(x_max(1)),'km']);
disp(['y的最小值 = ', num2str(y_min(1)),'km']);
disp(['y的最大值 = ', num2str(y_max(1)),'km']);
求得的 x x x的最小值为 974.8433 974.8433 974.8433，最大值为 982.2005 982.2005 982.2005； y y y的最小值为 717.1614 717.1614 717.1614，最大值为 733.1582 733.1582 733.1582。因此，最后得到的解是一个比较大的矩形区域，大致为 [ 975 , 982 ] × [ 717 , 733 ] [975, 982]×[717, 733] [975,982]×[717,733]。
五、模型3及求解
5.1 模型优化
模型2假设设备的测量误差是均匀分布的，这并不合理。一般来说，在多次测量中，应该假设设备的测量误差是正态分布的，而且均值为0。本例中给出的精度 σ i \sigma_{i} σi可以认为是测量误差的标准差。在这种理解下，用各自的误差限 σ i \sigma_{i} σi对测量误差进行无量纲化（也可以看成是一种加权法）处理是合理的，即求解如下的无约束优化问题更合理：
min E ( x , y ) = ∑ i = 1 3 ( α i − θ i σ i ) 2 + ( d 4 − ( x − x 4 ) 2 + ( y − y 4 ) 2 σ 4 ) 2 \min E(x, y)=\sum_{i = 1}^{3}(\frac{\alpha_{i}-\theta_{i}}{\sigma_{i}})^{2}+(\frac{d_{4}-\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}}}{\sigma_{4}})^{2} minE(x,y)=i=1∑3(σiαi−θi)2+(σ4d4−(x−x4)2+(y−y4)2)2
其中， tan α i = x − x i y − y i , i = 1 , 2 , 3 \tan \alpha_{i}=\frac{x - x_{i}}{y - y_{i}}, i = 1, 2, 3 tanαi=y−yix−xi,i=1,2,3。这仍然是一个非线性最小二乘拟合问题。
5.2 Matlab求解
在Matlab中，同样可以使用优化工具箱的函数来求解这个问题。代码如下：
% 定义测量数据
x0 = [746, 629, 1571];
y0 = [1393, 375, 259];
theta = [161.2, 45.1, 309.0] * pi / 180;
sigma_theta = [0.8, 0.6, 1.3] * pi / 180;
x4 = 155;
y4 = 987;
d4 = 864.3;
sigma_d4 = 2.0;% 定义目标函数
fun = @(xy) sum(((atan((xy(1) - x0)./(xy(2) - y0)) - theta)./sigma_theta).^2)+((d4 - sqrt((xy(1) - x4).^2+(xy(2) - y4).^2))./sigma_d4).^2;% 初始猜测值
xy0 = [1000, 900];% 使用fminsearch函数求解
[xy_opt, fval] = fminsearch(fun, xy0);x_opt = xy_opt(1);
y_opt = xy_opt(2);disp(['飞机位置坐标 x = ', num2str(x_opt),'km']);
disp(['飞机位置坐标 y = ', num2str(y_opt),'km']);
disp(['目标函数值 = ', num2str(fval)]);
启动Matlab的优化程序求解，得到全局最优解 x = 978.3071 x = 978.3071 x=978.3071， y = 723.9841 y = 723.9841 y=723.9841，对应的目标函数的值为 0.668035 0.668035 0.668035。这里得到的误差比模型1的大，这是因为模型1中使用的是绝对误差，而这里使用的是相对于精度 σ i \sigma_{i} σi的误差。对角度而言，分母 σ i \sigma_{i} σi很小，所以相对误差比绝对误差大，这是可以理解的。
六、总结
本文针对飞机精确定位问题，介绍了三种不同的数学模型及其Matlab求解方法。模型1基于简单的几何关系和最小二乘准则构建，能初步求解飞机位置，但未充分考虑测量误差特性；模型2考虑了测量误差的均匀分布，通过不等式约束确定飞机位置的取值范围；模型3则基于测量误差的正态分布假设，对误差进行无量纲化处理，使模型更加合理。
| 模型名称 | 模型核心思路 | 关键公式 | Matlab求解要点 | 求解结果 | 优缺点 |
|---|---|---|---|---|---|
| 模型1 | 基于几何关系构建方程，用最小二乘法使计算值与测量值的误差平方和最小，求解超定非线性方程组确定飞机位置 | tan θ i = x − x i y − y i \tan \theta_{i}=\frac{x - x_{i}}{y - y_{i}} tanθi=y−yix−xi d 4 = ( x − x 4 ) 2 + ( y − y 4 ) 2 d_{4}=\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}} d4=(x−x4)2+(y−y4)2 min J ( x , y ) = ∑ i = 1 3 ( x − x i y − y i − tan θ i ) 2 + [ d 4 − ( x − x 4 ) 2 + ( y − y 4 ) 2 ] 2 \min J(x, y)=\sum_{i = 1}^{3}(\frac{x - x_{i}}{y - y_{i}}-\tan \theta_{i})^{2}+[d_{4}-\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}}]^{2} minJ(x,y)=∑i=13(y−yix−xi−tanθi)2+[d4−(x−x4)2+(y−y4)2]2 | 定义测量数据和目标函数，以初始猜测值用fminsearch函数求解，可能得局部最优解 | x = 1019.306 x = 1019.306 x=1019.306， y = 987.2909 y = 987.2909 y=987.2909，目标函数值为 0.4729562 0.4729562 0.4729562 | 优点：计算简单直接。缺点：未考虑测量误差分布，解受 π \pi π取值影响大，易陷入局部最优 |
| 模型2 | 假设测量误差均匀分布，依据测量误差范围构建不等式组确定飞机坐标所在区域，求区域内 x x x、 y y y的最值 | tan ( θ i − σ i ) ≤ x − x i y − y i ≤ tan ( θ i + σ i ) \tan(\theta_{i}-\sigma_{i})\leq\frac{x - x_{i}}{y - y_{i}}\leq\tan(\theta_{i}+\sigma_{i}) tan(θi−σi)≤y−yix−xi≤tan(θi+σi) d 4 − σ 4 ≤ ( x − x 4 ) 2 + ( y − y 4 ) ≤ d 4 + σ 4 d_{4}-\sigma_{4}\leq\sqrt{(x - x_{4})^{2}+(y - y_{4})}\leq d_{4}+\sigma_{4} d4−σ4≤(x−x4)2+(y−y4)≤d4+σ4 | 定义测量数据和约束条件，分别以 x m i n x_{min} xmin、 x m a x x_{max} xmax、 y m i n y_{min} ymin、 y m a x y_{max} ymax为目标，用linprog函数求解 | x x x最小值为 974.8433 974.8433 974.8433，最大值为 982.2005 982.2005 982.2005； y y y最小值为 717.1614 717.1614 717.1614，最大值为 733.1582 733.1582 733.1582，解为矩形区域 [ 975 , 982 ] × [ 717 , 733 ] [975, 982]×[717, 733] [975,982]×[717,733] | 优点：考虑了测量误差范围。缺点：假设误差均匀分布不太符合实际，结果为较大矩形区域，定位不够精确 |
| 模型3 | 假设测量误差正态分布且均值为0，用误差限对测量误差无量纲化处理，构建无约束优化模型求解 | min E ( x , y ) = ∑ i = 1 3 ( α i − θ i σ i ) 2 + ( d 4 − ( x − x 4 ) 2 + ( y − y 4 ) 2 σ 4 ) 2 \min E(x, y)=\sum_{i = 1}^{3}(\frac{\alpha_{i}-\theta_{i}}{\sigma_{i}})^{2}+(\frac{d_{4}-\sqrt{(x - x_{4})^{2}+(y - y_{4})^{2}}}{\sigma_{4}})^{2} minE(x,y)=∑i=13(σiαi−θi)2+(σ4d4−(x−x4)2+(y−y4)2)2 tan α i = x − x i y − y i \tan \alpha_{i}=\frac{x - x_{i}}{y - y_{i}} tanαi=y−yix−xi | 定义测量数据和目标函数，以初始猜测值用fminsearch函数求解 | x = 978.3071 x = 978.3071 x=978.3071， y = 723.9841 y = 723.9841 y=723.9841，目标函数值为 0.668035 0.668035 0.668035 | 优点：基于更合理的误差分布假设，模型更科学。缺点：相对误差计算导致结果误差看似更大，理解和计算相对复杂 |