%使用卡爾曼濾波方法對飛行航班進行軌跡預測
%數據來源：FlightAware（https://zh.flightaware.com）
%航????班：CES9937?????寧波櫟社國際機場飛往成都雙流國際機場
%飛行時間：2018-07-17??07:07-09:55
%說明：取起飛后前20組數據作為實驗數據。對時間點進行近似取值，假設每隔30s上報一次數據
clear;
clc;
%采樣點的個數
N=228;
%測試數據：緯度
latitude=load（‘C:\workspace\matlab?workspace\CES9937\latitude.txt‘）;
%真實維度值
lat=latitude;
%卡爾曼濾波處理的狀態，即估計值
lat_kf=zeros（1N）;
%測報值
lat_z=zeros（1N）;
P=zeros（1N）;
%初始緯度值
lat（1）=29.8131;
%初始值的協方差
P（1）=0.09;
%初始測報值
lat_z（1）=29.8027;
%初始估計狀態。假設和初始測報值相同
lat_kf（1）=lat_z（1）;
%噪聲方差
%系統噪聲方差
Q=0.1;
%測量噪聲方差
R=0.001;
%方差決定噪聲大小
W=sqrt（Q）*randn（1N）;
V=sqrt（R）*randn（1N）;
%系統矩陣
F=1;
G=1;
H=1;
%本系統狀態為1維
I=eye（1）;
%模擬緯度測報，并濾波
for?k=2:N
????%隨時間推移，飛行緯度逐漸變化
????%k時刻的真是緯度值是測報儀器不知道的，測報值可能是無限接近于真實值，但并不是真實值
????%lat（k）=F*lat（k-1）+G*W（k-1）;
????%緯度在k時刻的測報值
????lat_z（k）=H*lat（k）+V（k）;
????%kalman濾波
????%有了k時刻的測報值lat_z（k）和k-1時刻的狀態，那么就可以進行濾波了
????%狀態預測
????lat_pre=F*lat_kf（k-1）;
????%協方差預測
????P_pre=F*P（k-1）*F‘+Q;
????%計算卡爾曼增益
????Kg=P_pre*inv（H*P_pre*H‘+R）;
????%新息
????e=lat_z（k）-H*lat_pre;
????%狀態更新
????lat_kf（k）=lat_pre+Kg*e;
????%協方差更新
????P（k）=（I-Kg*H）*P_pre;
end
%計算誤差
%測量值與真實值之間的偏差
Err_Messure=zeros（1N）;
%kalman估計與真實值之間的偏差
Err_Kalman=zeros（1N）;
for?k=1:N
????Err_Messure（k）=abs（lat_z（k）-lat（k））;
????Err_Kalman（k）=abs（lat_kf（k）-lat（k））;
end
t=1:N;
%濾波效果圖
figure
plot（tlat‘g-‘tlat_z‘b-‘tlat_kf‘r-‘）;
legend（‘真實值‘‘觀測值‘‘kalman濾波值‘）;
xlabel（‘測報時間點‘）;
ylabel（‘緯度值‘）;
%誤差分析圖
figure
plot（tErr_Messure‘b-‘tErr_Kalman‘r-‘）;
legend（‘測報偏差‘‘kalman濾波偏差‘）;
xlabel（‘測報時間點‘）;
ylabel（‘緯度偏差值‘）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????153993??2018-07-24?16:28??CES9937.docx
?????文件????????1825??2018-07-17?15:07??latitude.txt
?????文件????????2051??2018-07-17?15:15??longitude.txt
?????文件????????1988??2018-07-26?16:59??QA_KF_lat.m
?????文件?????????410??2018-07-17?11:16??飛行數據-20組.txt