Y. Chen, J. Sun, Y. Lin, G. Gui and H. Sari, "Hybrid n-Inception-LSTM-based aircraft coordinate prediction method for secure air traffic," IEEE Transactions on Intelligent Transportation Systems, vol. 23, no. 3, pp. 2773-2783, Mar. 2022, doi: 10.1109/TITS.2021.3095129.（Impact factors: 9.6）
航空大数据专题已经更新完毕，欢迎读者交流。上方为已发表论文，含有该专题的更多细节。
第一章：项目背景及数据集分析
第二章：神经网络输入/输出数据集制作
第四章：项目资源汇总及开源
目录
1、csv转换为matlab可用数组（此步很重要，影响后面训练集制作）
背景
了解ADS-B报文系统工作原理的都知道该系统是一种开放型的系统，只要有接收机任何人都能接收ADS-B报文，因此存在数据被篡改的风险。于是想出了一种仅使用ADS-B报文的信号强度（而非报文具体内容）通过神经网络预测飞机坐标的方法。该方法通过牺牲一定的预测准确性来换取预测可靠性（因为信号强度不易被篡改，但携带信息有限）。可用于飞行器搜救（飞行器遇难前ADS-B报文系统很有可能已被篡改，由ADS-B系统获得的定位可能是虚假定位）等领域。
项目来源
项目来自于OpenSky发布在AIcrowd上的任务：Aircraft Localization Competition（有梯子更快） （该竞赛已经完赛，公开了一些优秀的飞行器定位方案可以学习）
本文的目标与该竞赛的目标不同，本文仅使用该竞赛提供的数据集。
数据集
本文共涉及两个数据集：主数据集 和 传感器数据集。
数据集标签解释
点击此处查看官方给出的标签解释原文。
图1 主数据集
id：表示该条数据在该数据集中的收集顺序编号。
timeAtServer(时间戳)：以秒为单位，指OpenSky的服务器接收信息时间。 在每个1h主数据集中（官方提供了8个主数据集，每个数据集的搜集时间均为1小时），此时间戳从0开始，截止至3600秒，此时间戳具有毫秒精度（在传播延迟、处理延迟和网络延迟之后确定的）。
Aircraft(航班号)：每个主数据集中每个飞行器唯一的标识符（如：主数据集1中Aircraft 为123的飞行器与数据集2中Aircraft 为123的飞行器不是相同的飞行器，即不同数据集的数据无法合并）。
Latitude/longitude/geoAltitude(飞行器纬度/经度/地理高度)：由飞行器报告的当前位置坐标。纬度和经度以十进制度数（WGS84）表示，geoAltitude以米为单位。这些位置信息的准确性通常是未知的（有些飞行器会发送来自其内部系统而非GPS传感器的位置信息，在这种情况下，位置信息可能会漂移，甚至可能会相差几百米，但是这种情况很少见），但是大多数飞行器都以适当的精度（10米的精度）报告其位置。
baroAltitude(飞行器气压高度)：飞行器报告的气压高度（以米为单位）。气压高度取决于天气，可能与几何高度相差数百米，但是可以从“已知”飞行器上获悉差异，然后可以使用此信息估算几何高度或估算边界。
numMeasurements(应答传感器数量)：捕获到该飞行器信号的地面传感器数量。
Measurements(传感器测量数据)：每个捕获到该飞行器信号的地面传感器的测量数据。测量数据由包含三元组的JSON数组字符串表示。每个三元组包含：报告测量数据的传感器的全局唯一序列号、信号到达接收器的时间戳（以纳秒为单位）和接收信号强度指示器（RSSI）。RSSI的确切定义取决于接收器的类型，以dB为单位。纳秒级时间戳的精准度也取决于接收器的类型。对于型号为dump1090的传感器，这些时间戳通常不同步，并且仅具有12MHz的分辨率，不同步的时间戳会导致这些时间戳存在漂移（有时会很严重）。对于Radarcape和GRX1090传感器，这些时间戳与GPS同步（GPS同步意味着它们会不断重新同步以补偿时钟漂移），分辨率约为40-60MHz。
图2 OpenSky服务器系统
需要注意的是，主数据集涉及两个时间戳，分别是OpenSky服务器时间戳timeAtServer和Measurements中的传感器接收信号时间戳。如图2所示系统，飞行器间隔发出ADS-B报文，在接收范围内的接收机能够接收到此报文，最后通过网络上传至OpenSky服务器，由服务器汇总每个时间点各个接收机收到的对应航班的报文。OpenSky服务器和每个接收机（传感器）都有各自的时间戳，且有各自的时间基准，又或存在网络延迟导致时钟漂移，因此它们的时间戳不是同一个时间戳，无法合并。
图3 传感器数据集
Serial(序列号)：与图1主数据集中Measurements（传感器测量数据）栏数据相对应，是传感器的全局唯一序列号。
Latitude/longitude/height(地面传感器的纬度/经度/高度)：纬度和经度以十进制度（WGS84）为单位，高度以米为单位。此信息的准确性/不确定性未知，对于dump1090类型的传感器，这些位置由用户自己输入，大多数用户使用诸如Google Maps或智能手机之类的服务来确定天线的位置，此过程会出现误差，误差通常在数米范围内，但也有一些用户出于隐私原因，故意报告错误的位置，以伪装其确切位置，在这种情况下，位置甚至可能会出现几公里的误差。对于Radarcape或GRX1090类型的传感器，其位置更准确，因为它们由集成的GPS接收器自动报告，并且GPS天线通常位于靠近ADS-B天线的位置。
Type(传感器类型)：OpenSky当前支持基于开源S模式软件解码器dump1090和商用设备Radarcape，SBS-3和GRX1090的传感器。设备类型决定了测量噪声分布和时间戳精度（测量噪声分布和时间戳精度还取决于许多其他因素，例如：温度，内部/室外部署等）。
Good(接收数据的准确性)：True表示该传感器数据绝对可靠，False表示该传感器数据的可靠性存疑（具体误差未知）。
数据集分析
官方已经对数据进行了一定的分析，本文根据官方分析，使用matlab对官方分析进行复现，并做相关理解。
1、csv转换为matlab可用数组（此步很重要，影响后面训练集制作）
官方下载的数据集是csv格式，导入到matlab，为方便后续计算整合需要将已有csv转换成double类型矩阵。csv直接导入matlab会得到table类型矩阵，如图1，该主数据集中，A-H列数据可以直接转换成double类型数据，I列数据为string类型，需要对其分割并转换成double类型，如：
[[470,982753933,38],[499,6026974083.33333,30]] 要转成6列double类型数据
| 470 | 982753933 | 38 | 499 | 6026974083.33333 | 30 |
为方便管理数据，还需要将图1和图3中的数据整理到同一个矩阵中，重点需要整合图3中的每个传感器的坐标。如：[[470,982753933,38],[499,6026974083.33333,30]]，将传感器对应坐标整合入后得12列double类型数据
| 470 | 982753933 | 38 | 470号传感器纬度 | 470号传感器经度 | 470号传感器高度 | 499 | 6026974083.33333 | 30 | 499号纬度 | 499号经度 | 499号高度 |
为方便在传感器数据集中查找传感器对应坐标，先将图3传感器数据集中的A-D列转换成double类型矩阵（E、F列暂时用不着，先不管），存于变量sensors_double中，供后续使用。
整合图1和图3数据集的matlab代码如下：
load('sensors_double.mat');%加载double类型的传感器矩阵
temp=set_1;%set_1为直接从csv文件中读取的table类型矩阵，即图1中的table
[x,~]=size(temp);
measurements=9;%表示measurements所在列在第9列
allinfo=table2array(temp(:,1:8));%A-H行直接转换
out_i=0;
fprintf('初始化成功！\n')
for i = 1:x
if(~mod(i,1000))
fprintf('已完成 %d 行！\n',i);
end
str = temp(i,measurements);%取出i行9列string
str = table2array(str);
str = str{1:1};
str = str(3:length(str)-2);
S = regexp(str, '],[', 'split');
S = regexp(S, ',', 'split');
[~,m]=size(S);%传感器数据分割
k = measurements;%处理后的数据存放的起始位置（列）
out_i = out_i+1;%处理后的数据存放的起始位置（行）
for t = 1:m
[n,~]=find(sensors_double(:,1)==str2double(cell2mat(S{1,t}(1,1))));
sensors_id=n;%传感器在传感器数据集中的对应所在行
for l = 1:3
allinfo(out_i,k)=str2double(cell2mat(S{1,t}(1,l)));
k=k+1;
end
allinfo(out_i,k)=sensors_double(sensors_id,2);%加上该传感器的地理坐标
allinfo(out_i,k+1)=sensors_double(sensors_id,3);
allinfo(out_i,k+2)=sensors_double(sensors_id,4);
k=k+3;
end
end
整合后的数据被存入变量allinfo中，由于有多个主数据集，考虑到内存限制，可以对每个主数据集分开处理。（即allinfo是图1和图3中数据的整合）
2、分析服务器与传感器两大时间戳关系
研究时间戳关系只需取出其中一个航班进行观察即可，将矩阵allinfo按第3列升序排序（图1中的C列，Aircraft航班号），截取其中一个航班的所有数据，存入矩阵flightA_allinfo，供后续使用。（补充说明：原始数据默认按服务器时间戳升序排列，此时的航班号是乱的，若按航班号升序排列，会把相同的航班放一起，同时每个航班内按服务器时间戳升序排列，方便后续处理）
此时只分析传感器时间戳和服务器时间戳的关系，因此继续对矩阵flightA_allinfo处理，删除无用数据列，仅保留两大时间戳，存入变量flightA_timestamps，截取矩阵flightA_timestamps部分数据，如图4所示。
图4 矩阵flightA_timestamps部分数据情况
每列数据表示意义：
| 1 | 2 | 3 | 4 | 5 | ... | |
| 数据采集点1 | 服务器 时间戳 | 接收到该数据的第一个传感器id | 接收到该数据的第一个传感器时间戳 | 接收到该数据的第二个传感器id | 接收到该数据的第二个传感器时间戳 | ... |
| 数据采集点2 | ... | |||||
| ... |
如图4所示，第一个采集点共有4个传感器接收到数据，第三个采集点仅有2个传感器接收到数据，后面的数据列用0补充。因此矩阵flightA_timestamps的列数取决于所有采集点中能接收到数据的传感器的最大个数。
matlab实现服务器与传感器时间戳关系作图：
load('./point_color.mat');%随机生成的颜色矩阵，用于画图取色
load('./flightA_timestamps.mat');
temp=flightA_timestamps;
[n,m]=size(temp);
figure;
xlabel("Server's Timestamp (s)");
ylabel("Receiver's Timestamp (ns)");
set(gca,'FontName','Times New Roman','FontSize',20);
grid on;
hold on;
for i = 1:n
if(~mod(i,1000))
fprintf('已完成 %d 行！\n',i);
end
for t = 3:2:m
if temp(i,t)==0
break;%碰到0表示改行后面已无有效传感器数据，可以直接跳入下一行
end
%描点，不同传感器有不同的颜色
plot(temp(i,1),temp(i,t),'.','Color',[point_color(temp(i,t-1),1),point_color(temp(i,t-1),2),point_color(temp(i,t-1),3)],'LineWidth',1.5);
end
end
图5 传感器-服务器时间戳关系曲线
如图5所示，纵轴为各个传感器的接收信号时间戳（单位：纳秒），横轴为timeAtServer (OpenSky服务器接收时间戳，单位：秒)，该曲线体现了以航班A为例的各个传感器接收信号时间戳与OpenSky服务器接收信号时间戳的关系，图中不同颜色表示不同的传感器，从图5中能够发现有三个传感器偏移了其他传感器的曲线，这三个传感器存在时间戳的损坏或者不同步等问题。
通过矩阵flightA_timestamps继续研究时间戳差值，即替换描点代码为：
plot(temp(i,1),temp(i,t)/10e8-temp(i,1),'.','Color',[point_color(temp(i,t-1),1),point_color(temp(i,t-1),2),point_color(temp(i,t-1),3)],'LineWidth',1.5);
图6 传感器-服务器时间戳差值关系曲线
将图5中每个传感器的接收信号时间戳与服务器接收时间戳相减得到图6，该图能够更加直观地看出有三个传感器曲线偏移了其他传感器曲线。
在图6中还能看出，绝大多数的传感器时间戳与服务器时间戳的时间基准不同，即服务器时间戳为0时，传感器时间戳不一定为0，多数传感器时间戳与服务器时间戳在数值上相差0~100秒。再次印证了传感器时间戳相互之间独立，且与服务器时间戳也独立，它们有着各自的时间基准，无法直接合并使用。
3、分析传感器接收信号强度与接收距离关系
依旧利用矩阵flightA_allinfo，观察各个传感器接收到信号强度与接收距离的关系。matlab代码如下：（补充说明：该代码仅做数据可视化，用于定位“无效”传感器，但未做剔除，即没有把“无效”的传感器删除。若要删除，先记录下“无效”传感器id，后面做数据集时对这些“无效”传感器做跳过处理）
temp=flightA_allinfo;
time=temp(:,2);
temp(:,1:3)=[];
temp(:,3)=[];
temp(:,4)=[];
[~,m]=size(temp);
temp(:,5:6:m-4)=[];%删除无用数据
[~,m]=size(temp);
for i = 1:5:m-2
temp(:,i:i+2)=lla2ecef(temp(:,i:i+2));%将经纬度坐标转换成笛卡尔坐标
end
for i = 6:5:m-2
temp(:,i)=((temp(:,i)-temp(:,1)).^2+(temp(:,i+1)-temp(:,2)).^2+(temp(:,i+2)-temp(:,3)).^2).^0.5;%计算当前飞机坐标与传感器之间的距离（三维空间距离）
end
load('./point_color.mat');
[m,n]=size(temp);
figure;
xlabel("Distance (m)");
ylabel('Signal Strength (dB)');
set(gca,'FontName','Times New Roman','FontSize',20);
grid on;
hold on;
for i = 1:m
fprintf('i=%d\n',i);
for t = 4:5:n-4
if temp(i,t)==0%碰到0表示该行后面已无有效传感器数据，可以直接跳入下一行
break;
end
%单独观察315号传感器
%t+1为强度，t+2为距离
% if temp(i,t)==315
% plot(temp(i,t+2),temp(i,t+1),'.','Color',[point_color(temp(i,t),1),point_color(temp(i,t),2),point_color(temp(i,t),3)]);
% break;
% end
%显示所有传感器的强度
plot(temp(i,t+2),temp(i,t+1),'.','Color',[point_color(temp(i,t),1),point_color(temp(i,t),2),point_color(temp(i,t),3)]);
end
end如图7所示，该图体现了以航班A为例的各个传感器接收信号强度与接收距离之间的关系（不同的传感器用不同的颜色表示），纵轴表示传感器接收信号的强度（单位：dB），横轴表示传感器与飞行器之间的三维空间距离（单位：米）。从图7的整体分布能够看出，传感器与飞行器之间的三维空间距离越小，传感器接收强度越大，同时能够发现有些点不随距离变化而变化，这些传感器也许存在损坏。
图7 传感器接收信号强度与接收距离关系
4、分析采集到的数据集有效性
数据集的有效性主要包括：传感器数据有效性和飞行轨迹有效性。
（1）传感器数据有效性
传感器的筛选主要用于提高训练集或测试集输入数据的可靠性。本文优先使用标签Good(接收数据的准确性)为Ture的数据，由于标签为True的数据较少，为了尽可能利用更多的数据，对于剩余标签为False的数据，结合传感器的类型再进一步筛选。
根据数据集标签Type(传感器类型)的描述，可知各类传感器精准度排序为：GRX1090 = Radarcape > SBS-3 > dump1090。
Radarcape、GRX1090或SBS-3类型的传感器相比于dump1090类型的传感器精度更高，且自身坐标更准确，所以优先采用Radarcape、GRX1090或SBS-3的数据，最后再考虑dump1090的数据。
（2）飞行轨迹有效性
飞行器轨迹筛选主要用于提高训练集输出数据的可靠性。由于飞行器报告的当前位置坐标可能存在一定误差，本文对训练集的所有航班轨迹进行可视化人工筛选。如图8，例举了损坏的轨迹、信号中断的轨迹和良好的轨迹。损坏的轨迹和信号中断的轨迹显然需要从数据集中剔除，否则会影响后续训练性能。
（a）损坏的轨迹
（b）信号中断的轨迹
（c）良好的轨迹
图8 所采集的常见飞行轨迹例举
轨迹筛选过程完全由人工主观完成，没有统一的标准，适当保留一些无效的轨迹也能提高后续训练模型的鲁棒性。此过程很费时间，官方提供的数据集中存在大量无量轨迹，需要耐心筛选，源自于真实世界的数据确实会存在种种干扰，也是我们需要面对的一大挑战。
整个筛选过程在矩阵allinfo上进行，matlab代码如下（键盘输入d删除当前航班，a跳过，s删除上一个航班，b退出筛选）：
% 航班计数
% 先将allinfo按航班号升序排序，再进行如下操作
temp=allinfo;
clear allinfo;%allinfo占用内存太大，需及时清理
train_count=[];
train_id=[];
train_index=[];
[m,~]=size(temp);
reg=temp(1,3);%保存第一个航班号
count=0;
start=1;
for i = 1:m
if(~mod(i,10000))
fprintf('已完成 %d 行！\n',i);
end
if temp(i,3)==reg
count=count+1;
else
train_id=[train_id start];%每个航班号的起始索引
train_count=[train_count count];%每个航班号的个数
train_index=[train_index reg];%每个航班的航班号（经过多次筛选，航班号不一定连续）
start=start+count;
reg=temp(i,3);
count=1;
end
end
train_id=[train_id start];%每个航班号的起始索引
train_count=[train_count count];%每个航班号的个数
train_index=[train_index reg];%每个航班的航班号
% 筛选航班 （直接把图片关了也可以退出）
% 键盘输入d删除当前航班，a跳过, s删除上一个航班，b退出筛选
[~,m]=size(train_id);
a=input('是否要清空train_delete?（y/n）','s');
% 输入y/n
if a=='y'
train_delete=[];
% train_delete存放需要删除的航班号起始索引（根据实际情况判断是否要清空，初次筛选需要清空作为初始化矩阵）
end
input_='c';
set(gcf,'CurrentCharacter','c');
for i = 1:m
fprintf('%d、当前航班数据采集点个数 %d ,',i,train_count(i));
scatter3(temp(train_id(i):train_id(i)+train_count(i)-1,4),temp(train_id(i):train_id(i)+train_count(i)-1,5),temp(train_id(i):train_id(i)+train_count(i)-1,7),'r','.');
xlabel('Latitude /°');
ylabel('Longitude /°');
zlabel('Altitude /m');
grid on;
set(gca,'FontName','Times New Roman','FontSize',20);
while 1
waitforbuttonpress;%等待键盘输入
input_=get(gcf,'CurrentCharacter');
if input_=='d'
train_delete=[train_delete temp(train_id(i),3)];
fprintf('按键 %c \n',input_);
set(gcf,'CurrentCharacter','c');
break;
end
if input_=='s'
train_delete=[train_delete temp(train_id(i-1),3)];
fprintf('按键 %c \n',input_);
set(gcf,'CurrentCharacter','c');
break;
end
if input_=='a'
fprintf('按键 %c \n',input_);
set(gcf,'CurrentCharacter','c');
break;
end
if input_=='b'
fprintf('按键 %c \n',input_);
set(gcf,'CurrentCharacter','c');
break;
end
end
if input_=='b'
break;
end
end
close all;
%% 删除存放在变量train_delete中的对应航班数据
% aircraft_id=temp(:,3);
% [~,m]=size(train_delete);
% for i = 1:m
% fprintf('i= %d \n',i);
% index=find(aircraft_id==train_delete(i));
% temp(index(1):index(end),:)=nan;
% end
% temp(all(isnan(temp(:,1)),2),:)=[];% 最终矩阵temp即为筛选后的数据补充说明：代码上半部分是标记需要剔除的航班号，存于变量train_delete中，真正删除的操作需要打开最下面的注释，最终得到的temp便是剔除“无效”航班后的数据，重命名作为新的allinfo，供后面制作数据集使用。
筛选过程如下：
本章总结
通过对数据集的分析，发现这些由传感器真实采集的数据受多重因素干扰，会存在时间戳不同步、时钟漂移、传感器损坏、飞行轨迹损坏或信号中断等问题，为后续制作训练集造成很大的困难。
1445
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
