航空大数据 |
您所在的位置：网站首页 › ais报文 › 航空大数据 |
航空大数据
|
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） 航空大数据专题已经更新完毕，欢迎读者交流。上方为已发表论文，含有该专题的更多细节。 第一章：项目背景及数据集分析 第二章：神经网络输入/输出数据集制作 第三章：N-Inception-LSTM新型网络 第四章：项目资源汇总及开源 目录 背景 项目来源 数据集 数据集标签解释 数据集分析 1、csv转换为matlab可用数组（此步很重要，影响后面训练集制作） 2、分析服务器与传感器两大时间戳关系 3、分析传感器接收信号强度与接收距离关系 4、分析采集到的数据集有效性 本章总结 背景了解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类型数据 470982753933384996026974083.3333330为方便管理数据，还需要将图1和图3中的数据整理到同一个矩阵中，重点需要整合图3中的每个传感器的坐标。如：[[470,982753933,38],[499,6026974083.33333,30]]，将传感器对应坐标整合入后得12列double类型数据 47098275393338470号传感器纬度470号传感器经度470号传感器高度4996026974083.3333330499号纬度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部分数据情况 每列数据表示意义： 12345...数据采集点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，供后面制作数据集使用。 筛选过程如下：
通过对数据集的分析，发现这些由传感器真实采集的数据受多重因素干扰，会存在时间戳不同步、时钟漂移、传感器损坏、飞行轨迹损坏或信号中断等问题，为后续制作训练集造成很大的困难。 |
【本文地址】
今日新闻 |
点击排行 |
|
推荐新闻 |
|
图片新闻 |
|
专题文章 |