【matlab】雷達成像系列 之 BP（BackProjection，後向投影) 成像演算法
一、什麼是BP演算法？
由來：BP演算法最初是McCorkle受計算機層析技術的啟發推導而來，所謂的計算機層析：就是CT（Computer Tomograpy），這是在醫院中再普遍不過的技術了。
BP演算法的原理：BP演算法參考了“時延-疊加”的思想，在雷達應用中，其對雷達接收天線接收到的回波訊號進行距離向匹配率，獲取回波資料中包含的相幅資訊，再通過IFFT進行逆傅立葉變換，獲取收發天線組合的時延，最後累積訊號相干相加得到目標函式。
1.1 BP演算法
回波訊號與參與電訊號進行匹配濾波後，獲得的距離壓縮輸出訊號可表示為：
src(h)=Aexp{-j4pf0R(h)/c}
其中，A表示幅度；f0表示載波頻率；R( 表示雷達至點目標的距離h)/
二：BP成像演算法在SAR中的應用
2.1 優缺點分析：
BP演算法有一獨到的有點，其可以應用於多種架構的雷達天線而不受陣列形式的限制，這也使得其在雷達應用中顯得舉足輕重。
另一方面，BP演算法也有著巨大的缺陷，其運算資料量比較大，存在冗餘的現象，這直接導致了其計算效率低下（在作者模擬BP演算法對點目標進行成像的時候就發現需要等待的時間特別長，後文中將具體詳述，附圖+原始碼
2.2 FAST-BP演算法的特點
在雙站合成孔徑雷達（Synthetic Aperture Radar,SAR）中,BP演算法通過將雷達回波資料反向投影到成像區域的各個畫素，畫素值通過計算雷達回波在雷達天線和影象畫素之間的距離的延時來進行成像。
目前相關提升BP演算法速度的研究成果不是很顯著，但是在相關文獻上還是可以看到一些相應的快速BP演算法，這些演算法的大致有已下幾種套路：
1） 深究演算法結構，從理論層面對公式進行簡化，優化運算中可以改進的變換，保質保量。
2） 做出一些精確度上的犧牲，通過允許系統上誤差來對經典BP演算法進行修改或近似，比如西電李浩林研究的“機載SAR快速後向投影演算法”，國防科技大學電工院發表的《多級多分辨快速後向投影成像演算法》文章。
3） 基於子模組或域分割的簡化演算法，就是通過大化小的思想把成像區域分成一塊塊，從總體上降低演算法結構的複雜度。
這個作者沒有進行深究，讀者有興趣可自行查閱相關資料進行思考。歡迎評論區拍磚切磋！
三、單站SAR的BP演算法模擬
》》在這一段作者將盡全力對matlab原始碼進行附：解釋
3.1 單站SAR
首先作者先以單站SAR為例，這個資料比較齊全，解釋起來容易的多。參考教材《合成孔徑雷達與成像——演算法與實現》一書。（主要還是作者太懶了，不想自己畫，而且CSDN編輯器確實難用了點 
執行步驟：
例項：目標場景
好了，大夥就將就看哈，如上圖所示。
工作引數：成像區域網格:2^8&2^10（非嚴格尺寸，意思是對的。註明：圖取自網路，示意用）
運算時間：~4min~(欠火候)
成像區域範圍:A到左邊緣20m,C到右邊緣20m；
》》》》》上邊緣縱座標50m，下邊緣縱座標-50m
3.2 圖及matlab原始碼（關鍵部分詳述，有問題請在評論區諮詢）
！ ! ! ! -_-hahah
clear all; %清除工作空間的內容
%MATLAB有個基本的工作空間，用base標識，此外，當開啟一個函式m檔案時，可能會產生很多工作空間。
%每一個函式對應一個工作空間。例如，一個圖形使用者介面程式test，可能有test、gui_mainfcn、
%pushbutton1_callback等工作空間。這時，如果呼叫clear命令時，需要注意了：%
clc;%清除命令視窗內容</pre><pre class="sql" name="code">% clc命令是用來清除命令視窗的內容，這點不用多說。
% 不管開啟多少個應用程式，命令視窗只有一個，所以clc無論是
% 在指令碼m檔案或者函式m檔案呼叫時，clc命令都會清除命令視窗的內容。clear;
%fmin=？; %請讀者自行賦值%fmax=？; %請讀者自行賦值
c=3e8; M=901;
kmin=2*fmin/c; %波數
B=2*(fmax-fmin)/c; %頻寬
f_step=fmax-fmin/M-1; %步進頻率，步進頻的好處在於雷達發射此形式訊號可以得到目標在各頻點的幅度和相位資訊R0=10; %希望在評論區看到讀者的答案！！
R_v=1.5; %希望在評論區看到讀者的答案
rx1=-0.4; %希望在評論區看到讀者的答案
rx2=36; %希望在評論區看到讀者的答案
rx3=0; %希望在評論區看到讀者的答案
rx4=0.9; %希望在評論區看到讀者的答案
rx5=0.3; %希望在評論區看到讀者的答案
N=161;
angle_step=0.1; %角度最小分量
for n=1:N;
sita=(n-(N+1)/2)*angle_step*pi/180; %sita值表示：
R1(1,n)=sqrt(R0^2+rx1^2-2*R0*rx1*cos(sita+pi/2)); %求一個神奇公式的平方根，
R2(1,n)=sqrt(R0^2+rx2^2-2*R0*rx2*cos(sita+pi/2));
R3(1,n)=sqrt(R0^2+rx3^2-2*R0*rx3*cos(sita+pi/4));
R4(1,n)=sqrt(R0^2+rx4^2-2*R0*rx4*cos(sita));
R5(1,n)=sqrt(R0^2+rx5^2-2*R0*rx5*cos(sita+pi/4));
end
for m=1:M
k=2*pi*(8e9+(m-1)*f_step)/c;
% E(1:N,m)=exp(-2j*k*(R1(1,:)))./R0+exp(-2j*k*(R2(1,:))).
%/R0+exp(-2j*k*(R3(1,:)))./R0+exp(-2j*k*(R4(1,:)))./R0+exp(-2j*k*(R5(1,:)))./R0;
% E(1:N,m)=exp(-2j*k*(R1(1,:)))./R1(1,:);
E(1:N,m)=exp(-2j*k*(R1(1,:)))./R1(1,:)+exp(-2j*k*(R2(1,:)))./R0;
end
for n=1:N
E1(n,:)=ifft(E(n,:)); %求逆傅立葉變換
end
Lt=(0:900)/B;
plot(Lt,abs(E1(81,:)));
F=4;
MF=F*(M-1)+1;
N1=round((R0-R_v)*B*F);
N2=round((R0+R_v)*B*F);
D_N=N2-N1+1;
for n=1:N % per det_sita
for m=1:M % per det_f
G(1,m)=(B*(m-1)./(M-1)+kmin).*E(n,m)*( 0.54-0.46.*cos(2*pi*m/M) ); %fHamming;對矩陣的每行用fHamming進行加權。
end
P_sita(n,:)=ifft(G(1,1:M),MF);
P_sita1(n,1:D_N)=P_sita(n,N1:N2); % 擷取目標區
end
dx=0.005;
dy=dx;
DN=floor(sqrt(2)*R_v/dx);
Fig=zeros(DN,DN);
h=waitbar(0,'計算中....'); %由於BP演算法計算效率低下，所以特意加了進度條一方讀者長時間看不到出圖誤以為程式出錯
for xn=1:DN
xxn=(xn-(DN+1)/2)*dx;
for yn=1:DN
yyn=(yn-(DN+1)/2)*dy;
for n=1:N
sita=(n-(N+1)/2)*angle_step*pi/180;
ham=0.54-0.46*cos(2*(n-1)*pi/(N-1)); %漢明窗
L0=sqrt(R0.^2+xxn.^2+yyn.^2+2*R0*(yyn*cos(sita)-xxn*sin(sita)))-R0; %中場
% L0=yyn*cos(sita)-xxn*sin(sita); % 遠場
LB=L0*B*(MF/M)+D_N/2+1;
d_L=LB-floor(LB);
P_L=P_sita1(n,floor(LB))+d_L*(P_sita1(n,floor(LB)+1)-P_sita1(n,floor(LB)));
Fig(xn,yn)=Fig(xn,yn)+P_L*exp(i*2*pi*kmin*L0)*ham;
end
end
waitbar(xn/DN,h);
end
close(h)
x=(1-(DN+1)/2)*dx:dx:(DN-(DN+1)/2)*dx;
y=x;
figure; surf(x,y,abs(Fig));shading flat; %shading 是用來處理色彩效果的，分以下三種：
%no shading 一般的預設模式 即shading faceted
%shading flat 在faceted的基礎上去掉圖上的網格線
%shading interp 在flat的基礎上進行色彩的插值處理，使色彩平滑過渡
xlabel('X');ylabel('Y');
感謝各位讀者的支援，作者會努力提高部落格水準，歡迎讀者對錯誤或有待改進的地方提出建議和意見 ——》》PS《《——話說現在寫部落格花的時間越來越多了，好費時間！小碩，努力積累一點一滴。