2022年国赛高教杯数学建模
B题 无人机遂行编队飞行中的纯方位无源定位
原题再现
无人机集群在遂行编队飞行时，为避免外界干扰，应尽可能保持电磁静默，少向外发射电磁波信号。为保持编队队形，拟采用纯方位无源定位的方法调整无人机的位置，即由编队中某几架无人机发射信号、其余无人机被动接收信号，从中提取出方向信息进行定位，来调整无人机的位置。编队中每架无人机均有固定编号，且在编队中与其他无人机的相对位置关系保持不变。接收信号的无人机所接收到的方向信息约定为：该无人机与任意两架发射信号无人机连线之间的夹角（如图 1 所示）。例如：编号为 FY01、FY02 及 FY03 的无人机发射信号，编号为FY04 的无人机接收到的方向信息是 𝛼1，𝛼2 和 𝛼3。
请建立数学模型，解决以下问题：
问题 1 编队由 10 架无人机组成，形成圆形编队，其中 9 架无人机（编号 FY01~FY09）均匀分布在某一圆周上，另 1 架无人机（编号 FY00）位于圆心（见图 2）。无人机基于自身感知的高度信息，均保持在同一个高度上飞行。
(1) 位于圆心的无人机（FY00）和编队中另 2 架无人机发射信号，其余位置略有偏差的无人机被动接收信号。当发射信号的无人机位置无偏差且编号已知时，建立被动接收信号无人机的定位模型。
(2) 某位置略有偏差的无人机接收到编号为 FY00 和 FY01 的无人机发射的信号，另接收到编队中若干编号未知的无人机发射的信号。若发射信号的无人机位置无偏差，除 FY00 和 FY01外，还需要几架无人机发射信号，才能实现无人机的有效定位？
(3) 按编队要求，1 架无人机位于圆心，另 9 架无人机均匀分布在半径为 100 m 的圆周上。当初始时刻无人机的位置略有偏差时，请给出合理的无人机位置调整方案，即通过多次调整，每次选择编号为 FY00 的无人机和圆周上最多 3 架无人机遂行发射信号，其余无人机根据接收到的方向信息，调整到理想位置（每次调整的时间忽略不计），使得 9 架无人机最终均匀分布在某个圆周上。利用表 1 给出的数据，仅根据接收到的方向信息来调整无人机的位置，请给出具体的调整方案。
问题 2 实际飞行中，无人机集群也可以是其他编队队形，例如锥形编队队形（见图 3，直线上相邻两架无人机的间距相等，如 50 m）。仍考虑纯方位无源定位的情形，设计无人机位置调整方案。
整体求解过程概述(摘要)
本文主要研究无人机纯方位无源定位，建立定位模型和调整方案模型完成有效定位被动接收信号无人机以及进行纯方位无源定位调整工作。
针对问题一(1)，首先根据3架无人机相对位置信息进行分类讨论，再结合正弦定理以及其他几何关系，确定六种定位模型，采用遍历算法，固定一架发射信号无人机，遍历另外两架无人机的相对位置关系，求解每种情况下所有位置上被动接收无人机的实际坐标，输出图解结果，发现定位效果良好，模型合理。
针对问题一(2)，在无人机知道自身编号与其余无人机的相对位置关系的假设条件下，确定位置编号发射信号无人机的相对位置，再推断其准确标号。随后结合问题一(1)定位模型计算被动接受无人机坐标，确定定位是否有效。利用仿真模拟生成的无人机实际坐标样本，验证得仅需通过增加一架无人机发射信号就能够有效定位接受信号无人机位置。
针对问题一(3)，分析给定初始坐标值，调整方案分为三步骤，步骤一:利用问题一(1)的定位模型，调整同一圆周上的三架无人机的相对位置接近理想位置关系。步骤二:以步骤一中三架无人机为参照系，利用问题一(1)的定位模型，调节余下无人机的角度坐标，使之接近理想相对位置，输出角度调节量；步骤三:同步骤二参照基准，建立几何关系模型，调节余下无人机的长度坐标每一步骤都对应各自独立的调整模型，模型求解利用迭代算法进行。通过输出调整后的坐标位置信息，对比初始位置信息确认此方案模型可行、实用。
针对问题二，将锥形编队分解成多个同心圆的圆形编队，利用问题一(3)中给出的调节方案调整各圆周上无人机，针对特殊点通过共线无人机制定相应的调整方案。
模型假设：
1、假设无人机发射信号与接收信号稳定不受干扰。
2、假设忽略被动接受无人机处理信号的时间。
3、假设每架无人机知道自身编号以及相对位置关系。
问题分析：
问题一(1)的分析
问题一首先根据无偏差的位置坐标确定圆周上两架发射信号无人机与单架接收信号无人机的相对位置关系，根据相对位置关系进行定位模型的分流，其中定位模型根据圆周上的三架无人机以及圆心无人机(FY00)之间的几何关系，利用正弦定理建立。最后在此定位模型的框架下，固定一架发射信号无人机位置，遍历第二架发射信号无人机和被动接收信号无人机的无偏位置坐标的所有可能，随机生成无偏坐标基础上的有偏接收无人机坐标值，求解出每一种可能情况下被动接收信号无人机的位置信息。
问题一(2)的分析
问题二首先确定来自编号FY01无人机发出的信号，此过程需要圆周上另外的两架无人机。通过两组无人机的三组角度信息确定FY01无人机发出的信息后，再从两架无人机中随机选择一架所提供的角度信息通过几何关系判断发射源与接收源的相对位置（此处假设被动接收信号无人机知道自身编号），根据相对位置的条件再进一步确定未知编号发射信号无人机的编号，从而利用已知的角度信息和坐标信息，代入问题一所得的定位模型，计算出被动接收信号无人机的实际坐标，建立偏差值的范围判定定位有效性，将实际坐标与无偏坐标作差并对比偏差值的范围，得出结论。
问题一(3)的分析
问题三首先分析给定的无人机初始位置表，确定基准圆周，本题考虑FY02(ρ2,θ2)、FY05(ρ5,θ5)、FY08(ρ8,θ8)所在半径98m的圆作为调整方案的理想圆周。调整方案大致分为三大步骤：一、θ2、θ5、θ8的调整，使得无人机FY02、FY05、FY08的相对位置接近理想值或达到理想值。二、由步骤一得到θ_2´、θ_5´、θ_8^´将其作为发射信号无人机，且为圆周上的参考系，调整余下圆周上所有无人机的角度坐标θk, ,k=1,3,﹍,9。三、经过前两步骤的调整，此时所有无人机均匀分布在以无人机FY00为圆心的同心圆上，以无人机FY02、FY05、FY08为参照系，利用几何关系调整其余无人机的长度坐标ρk,k=1,3,﹍,9，使得所有无人机均匀分布在半径为98m的圆上。
问题二的分析
问题二首先在锥形编队内找出一个无人机，使得以它作为圆心，划分多个同心圆，是尽可能多的无人机的理想坐标在圆周上，此无人机编号确定为FY05，随后分组调整无人机，每组包括的无人机为相互所成圆形编队的无人机，最后套用问题一（3）的方案模型，对各个圆形编队的无人机进行调整。唯一单独考虑的无人机编号为FY11和FY15，采用调整位置后的两组共线发射信号的无人机与FY05所成的角度信息的变化关系进去制定调整方案。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import warnings
#wanings.filterwarnings(‘ignore’)
plt.reParams[‘font.sans-serif’]=[‘SimHei’] #用来正常显示中文标签 plt.rcParams[‘axes.unicode_minus]=False #用来正常显示负号
class B_1_cls:
def_init_(self)
self.theta_rad=None
self.theta_deg=None
self.r=None
self.all_color=['b','e','g','k','m','y',' purple']
self.color=None
self.counter=0 #计时器，画图用
def_mid_k_funcself,theta_i, theta_j,theta_k_true):
“””
此时k在中间,偏差为+-1
“””
Theta_rad=None
Theta_deg=None
r = None
color=None
if theta_i-theta_k_true<180 and theta_k_true-theta_i<180:
alpha_i=(180-theta_k_true+theta_i)/2
alpha_i=np.arange(alpha_i-1,alpha_i+1,0.1)
alpha_j=(180+theta_k_true-theta_j)/2
alpha_j=np.arange(alpha_j-1,alpha_j+1,0.1)
#print(alpha_i)
#print(alpha_j)
#print(alpha_i.shape)
fenzi=np.sin(np.radians(alpha_j+theta_j))*np.sin(np.radians(alpha_i))+np.sin
np.radians(-alpha_i+theta_i))*np.sin(np.radians(alpha_j))
fenmu=np.cos(np.radians(alpha_j+theta_j))*np.sin(np.radians(alpha_i))+np.cos(np.radians(-alpha_i+theta_i))*np.sin(np.radians(alpha_i))
theta_rad=np.array(pd.Series(np.arctan(fenzi/fenmu)).map(lambdax:x if x>0 else x+np.pi))
theta_deg=np.rad2deg(theta_rad)
r=100*np.sin(np.radians(180-theta_j+theta_deg-alpha_i))/np.sin(np.radians(alpha_j))
#print("================")
#print(theta_deg)
#print(r)
#print("================")
#print(theta_deg shape)
#print(r.shape)
#print("theta为:{}".format(theta_deg))
#print("r 为:{}".format(r))
elif theta_j-theta_k_true> 180 and theta_k_true-theta_i<180:
alpha_i=(180-theta_k_true+theta_i)/2
alpha_i=np.arange(alpha_i-1,alpha_i+1,0.1)
alpha_j=(-180-theta_k_true+theta_j)/2
alpha_j=np.arange(alpha_j-1,alpha_j+1,0.1)
fenzi=np.sin(np.radians(-alpha_j+theta_j))*np.sin(np.radians(alpha_i))-np.sin(np.radians(-alpha_i+theta_i))*np.sin(np.radians(alpha_j))
fenmu=np.cos(np.radians(-alpha_j+theta_i))*np.sin(np.radians(alpha_i))-np.cos(np.radians(-alpha_i+theta_i))*np.sin(np.radians(alpha_j))
theta_rad=np.array(pd.Series(np.arctan(fenzi/fenmu)).map(lambda x:x if x>0 else x+np.pi))
theta_deg=np.rad2deg(theta_rad)
#print("======================")
#print(theta_deg)
#print("======================")
r=100*np.sin(np.radians(-180-theta_deg+theta_j-alpha_j))/np.sin(np.radians(alpha_j))
elif theta i-theta ktrue<180 and theta ktrue-thetai>180
alpha_i=(-180+theta ktrue-theta_i)/2
alpha_i=np.arange(alpha_i-1alpha_i+1,0.1)
alpha_j=(180+theta_k_true-theta_j)/2
alpha_j=np.arange(alpha_j-1alpha_j+1,0.1)
fenzi=np.sin(np.radians(alpha_i+theta i))*npsin(npradians(alpha_j))-np.sin(np.radians(alpha_j+theta_j))*np.sin(np.radians(alpha_i))
fenmu=np.cos(np.radians(alpha_i+theta_i))*np.sin(np.radians(alpha_j))np.cos(np.radians(alphaj+theta_j))*np.sin(np.radians(alpha_i))
theta_rad=nparray(pd.Series(nparctan(fenzifenmu))map(lambda x:x+np.pi if x>0 else x+2*np.pi))
theta_deg=np.rad2deg(theta_rad)
r=100*np.sin(np.radians(180+theta_deg-theta_j-alpha_j cr/np.sin(np.radians(alpha_i))
#print"===============”)
#print(theta_deg)
#print(r)
#print("=====================")
return theta_rad, theta_deg. r
def right or left kfunc(selftheta itheta jtheta k true)
“””
此时k在两边偏差为+-1
“””
theta_rad=None
theta deg=None
r=None
if theta_k_true-theta_j<180 and theta_k_true-theta_i<180
alpha_i=(180-theta_k_true+theta_i)/2
alpha_i=np.arange(alpha_i-1alpha_i+1,0.1)
alpha_j=(180-theta ktrue+theta_i)/2
alpha_j=np.arange(alpha_j-1alpha_j+1,0.1)
fenzi=np.sin(np.radians(alpha_j-theta_j))*np.sin(np.radians(alpha_i))-np.sin(np.radians(alpha_i-thetai))*np.sin(np.radians(alpha_i))
fenmu=np.cos(np.radians(alpha_i-theta_i))*np.sin(np.radians(alpha_i))-np.cos(np.radians(alpha_j-theta_i))*np.sin(np.radians(alpha_i))
theta_rad=np.array(pd.Series(np.arctan(fenzi/fenmu)).map(lambda x:x if x>0 else x+np.pi))
theta_deg=np.rad2deg(theta_rad)
r=100*np.sin(np.radians(-theta_j+theta_deg+alpha_j))/ np.sin(np.radians(alpha_j))
elif theta_k_true -theta_j<180 and theta_k_true-theta_i>180:
alpha_i=(-180+theta_k_true-theta_i)/2
alpha_i=nparange(alpha_i-1alpha_i+1,0.1)
alpha_j=(180-theta_k_true+theta_i)/2
alpha_j=np.arange(alpha_j-1alpha_j+1,0.1)
fenzi=np.sin(np.radians(-alpha_j+theta_j))*np.sin(np.radians(alpha_ i))+np.sin(np.radians(alpha_i+theta_i))*np.sin(np.radians(alpha_j))
fenmu=np.cos(np.radians(alpha_i+theta_i))*np.sin(np.radians(alpha _j))+np.cos(np.radians(alpha_j-theta_j))*np.sin(np.radians(alpha_i))
theta_rad=np.array(pd.Series(np.arctan(fenzi/fenmu)).map(lambda x:x+np.pi if x>0 else x+2*np.pi))
theta_deg=np.rad2deg(theta_rad)
r=100*np.sin(np.radians(theta_deg-theta_j+alpha_j)) / np.sin(np.radians(alpha_j))
elif theta_k_true-theta_j>180 and theta_k_true-theta_i>180:
alpha_i=(-180+theta_k_true-theta_i)/2
alpha_i=np.arange(alpha_i-1alpha_i+1,0.1)
alpha_j=(-180+theta_k_true-theta_j)/2
alpha_j=nparange(alpha_j-1alpha_j+1,0.1)
fenzi=npsin(np.radians(alphaj+thetaj))*np.sin(np.radians(alpha i))- np.sinnp.radians(alpha i+theta i))*np.sin(np.radians(alpha j))
fenmu=npcos(np.radians(alphaj+thetaj))*np.sin(np.radians(alpha i))-np.cos(np.radians(alpha_i+theta_i))*np.sin(np.radians(alpha j))
theta_rad=np.array(pd.Series(np.arctan(fenzi/fenmu)).map(lambda x:x+np.pi if x> 0 else x+ 2* np.pi))
theta deg=np.rad2deg(theta_rad)
r=100*np.sin(np.radians(-180+theta_deg-theta_i-alpha_i))/ np.sin(np.radians(alpha_j))
return theta_rad, theta_deg, r
def polar_plot(self):
Qi=0
Qj=np.arange(40,360,40)
Qk=np.arange(40,360,40)
for i, j in enumerate(Qj):
ax=plt.subplot(2,4,i+1,projection=’polar’)
self.theta_rad=np.array([])
self.theta_deg=np.array([])
self.r=np.array([])
self.color=nparray([])
self.counter=0
for k in Qk:
if k ==j:
continue
else:
if j> k:
theta_rad, theta_deg,temp_r=self. __mid_k_func(Qi,j,k)
self.theta_rad=np.hstack([selftheta radtheta rad])
self.theta deg=nphstack([selftheta degtheta deg]) selfr=np.hstack([selfr,temp_r]) self.color=np.hstack(
[self.color,np.array([selfallcolor[self.counter] for_in range(theta rad.shape[0])))
elif j< k:
theta rad theta degtempr=self. right or left k func(Qi,j.k)
self.theta rad=nphstack([selftheta radtheta rad)
selftheta deg=np.hstack([selftheta degthetadeg]) selfr=np.hstack([selfr,temp_r])
self.color=np.hstack(
[self.color,np.array ( [self.all color[selfcounter] for_in range(theta_radshape[0])])])
self.counter+=1
#print(self.theta_deg)
# print(selfr)
#print(“==========================”)
#print(“==========================”) #print(“==========================”)
#画图
# ax.scatter(self.theta_rad,selfrlabel=’定点i极坐标’, s=10,c=self.color ax.scatter(self.theta_rad,self.r,s=10,c=self.color)
ax.set_rlim(0,105)
ax.scatter(Qi,100,s=100,c=Y,marker=’*’)
if i== 0:
ax.scatter(np.radians(j),100,s=100,c=’r’,marker=’*’,labe=’发射信号无人机极坐标’)
ax.legend(bboxto anchor(551.35))
#ax legend(bbox to anchor(01.35))
ax.scatter(np.radians(j),10s=100c=,marker*)
#plt.title("极坐标图"，loc='left')
plt.show()
# @classmethod
Def b_2(self,j,k)
return self. __mid_k_func(O,j.k) if j>k else selfright or left k func(0,j,k)
if_name__==main_
b_1=B_1_cls0)
b_1.polar plot()
