背景

了解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类型数据

为方便管理数据，还需要将图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('./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

整合后的数据被存入变量allinfo中，由于有多个主数据集，考虑到内存限制，可以对每个主数据集分开处理。（即allinfo是图1和图3中数据的整合）

2、分析服务器与传感器两大时间戳关系

研究时间戳关系只需取出其中一个航班进行观察即可，将矩阵allinfo按第3列升序排序（图1中的C列，Aircraft航班号），截取其中一个航班的所有数据，存入矩阵flightA_allinfo，供后续使用。（补充说明：原始数据默认按服务器时间戳升序排列，此时的航班号是乱的，若按航班号升序排列，会把相同的航班放一起，同时每个航班内按服务器时间戳升序排列，方便后续处理）

此时只分析传感器时间戳和服务器时间戳的关系，因此继续对矩阵flightA_allinfo处理，删除无用数据列，仅保留两大时间戳，存入变量flightA_timestamps，截取矩阵flightA_timestamps部分数据，如图4所示。

图4 矩阵flightA_timestamps部分数据情况

每列数据表示意义：

如图4所示，第一个采集点共有4个传感器接收到数据，第三个采集点仅有2个传感器接收到数据，后面的数据列用0补充。因此矩阵flightA_timestamps的列数取决于所有采集点中能接收到数据的传感器的最大个数。

matlab实现服务器与传感器时间戳关系作图：

图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，供后面制作数据集使用。

筛选过程如下：

本章总结

通过对数据集的分析，发现这些由传感器真实采集的数据受多重因素干扰，会存在时间戳不同步、时钟漂移、传感器损坏、飞行轨迹损坏或信号中断等问题，为后续制作训练集造成很大的困难。

原文链接：https://blog.csdn.net/qq_39291503/article/details/117673374