自聚焦算法

自聚焦算法

1.算法基本原理和推导

       自聚焦算法。假设一幅退化的雷达图像可以表示为$g(t_m,t_n)=I(t_m,t_n)+j\cdot Q(t_m,t_n)$。其中，考虑这一种典型的情况，图像的退化，是由于横向频率域中有误差相位引起的，有
$$g(t_m,t_n)=\int G(t_m,t_n)\cdot e^{j\cdot \varphi (f_n)}\cdot e^{j2\pi f_n{t}_n}d{f}_n\text{\hspace{1em}(1)}$$
作为最简单的一种情况，我们认为，误差函数$\varphi (f_n)$是一个二次函数，有
$$\varphi (f_n)=a\cdot f_n^2\text{\hspace{1em}(2)}$$
对于我们而言，找到对应的$a$就可以解决问题。
       实际过程中，我们采用图像熵来衡量图像的质量，熵值定义为
$$E=-\sum _{m=1}^M\sum _{n=1}^N(\frac{{|g(m,n)|}^2}{E_n}\cdot lo{g}_2\frac{{|g(m,n)|}^2}{E_n})\text{\hspace{1em}(3)}$$
其中，
$$E_n=\sum _{m=1}^M\sum _{n=1}^N|g(m,n){|}^2\text{\hspace{1em}(4)}$$
$g(m,n)$是对$g(t_m,t_n)$离散化的结果。
        对于函数的导数，我们可以采用近似
$$f^{{}^{\prime }}(x)=\underset{\Delta \to 0}{\lim }\frac{f(x+\Delta x)-f(x)}{\Delta x}\approx \frac{f(x+\Delta x)-f(x)}{\Delta x}\text{\hspace{1em}(5)}$$

2.级数反演法原理

        假设一函数$f(x)$在$x_0$处可以用泰勒级数展开为：
$$f(x)=\sum _{n=0}^{\infty }\frac{1}{n!}{f}^n(x_0)(x-x_0{)}^n$$
        则：
$$\begin{array}{rl}{f}^{{}^{\prime }}(x)=& \sum _{n=0}^{\infty }\frac{1}{n!}{f}^n(x_0)\cdot n\cdot (x-x_0{)}^{n-1}\\ =& \sum _{n=1}^{\infty }\frac{1}{(n-1)!}{f}^n(x_0)(x-x_0{)}^{n-1}\\ =& \sum _{n=0}^{\infty }\frac{1}{n!}{f}^{n+1}(x_0)(x-x_0{)}^n\end{array}$$
$$f^{{}^{\prime }}(x)-f^{{}^{\prime }}(x_0)=\sum _{n=1}^{\infty }{a}_n(x-x_0{)}^n{a}_n=\frac{1}{n!}{f}^{n+1}(x_0)\text{\hspace{1em}(6)}$$
        假设函数$f^{{}^{\prime }}(x)$在$x_0$的某一邻域内单调，则函数$f^{{}^{\prime }}(x)$在此邻域内存在反函数$f^{{}^{\prime }-1}(x)$。同样利用泰勒级数$f^{{}^{\prime }-1}(x)$在$f^{{}^{\prime }}(x_0)$展开为:
$$x-x_0=\sum _{m=1}^{\infty }{A}_m(f^{{}^{\prime }}(x)-f^{{}^{\prime }}(x_0){)}^m\text{\hspace{1em}(7)}$$
        令
$$Y=f^{{}^{\prime }}(x)-f^{{}^{\prime }}(x_0)y=x-x_0$$
       联立$(6)(7)$可得方程组:
$$\{\begin{array}{l}Y={\sum }_{n=1}^{\infty }{a}_n{y}^n\\ y={\sum }_{m=1}^{\infty }{A}_m{Y}^m\end{array}$$
       可得：
$$Y=\sum _{n=1}^{\infty }{a}_n(\sum _{m=1}^{\infty }{A}_m{Y}^m{)}^n$$
       对应项系数相等,可解得各系数$A_m$:
$$\{\begin{array}{l}1=a_1{A}_1\\ 0=a_1{A}_2+a_2{A}_1^2\\ 0=a_1{A}_3+a_2(A_1{A}_2+A_2{A}_1)+a_3{A}_1^3\\ 0=a_1{A}_4+a_2(A_1{A}_3+A_2^2+A_3{A}_1)+a_3(3\cdot A_1{A}_1{A}_2)+a_4{A}_1^4\\ \cdots \end{array}$$

3.实验涉及到数据和结果

       求得的最优值$a=860.2245$。恢复的图像中右边还是有点问题（可以看到它是整个图像循环向左平移了一部分，感兴趣的同学可以看一下matlab代码。）

4.附录：对应的matlab代码

%%%%%%%%%%%%%%%%%%%
%2020.09.16 李昕
%%%%%%%%%%%%%%%%%%%

clc, clear;
load image_bad
%展现原始图图像
figure();
image(abs(image_bad));
colormap(gray(256));

image_row_fft = fft(image_bad, [], 2);
[M,N] = size(image_row_fft);

%%利用枚举法画出a-E的曲线图,粗略确定a的范围
a_begin = 850;
a_end = 870;
delta = 0.05;
index = 1;
E_mat = [];
for a = a_begin:delta:a_end
    E_mat(index) = imageEntropy(a, N, image_row_fft, false);
    index = index + 1;
end
figure();
plot(a_begin:delta:a_end, E_mat);
% 利用枚举法结果图形画出a=860时恢复的图片
imageEntropy(860.05, N, image_row_fft, true);


%利用级数反演法求解最优图片熵值对应的a值，精确求解a的值

%从图形上看，它是一个多峰的函数，因此初始点的选取就很重要
%%%%初始值选取->优化求解结果
%%%%853->860.2245（求解结果对应极小值点） 
%%%%（865-870）->869.8797（求解结果对应极大值点） 
%%%%855->864.4170（求解结果对应极大值点）
%%%%当然有很多初始点求解的结果会发散

%虽然不同的初始值可以得到不同的结果，我们应该从中选取对应图片熵最小对应的a值
a = 853;
delta_a = 0.1;
cycle = 50;
for c = 1:cycle
    a_mat = zeros(1, 5);
    %求1-4次导
    for plus_delta_a = 0:4
        a_mat(plus_delta_a+1) = imageEntropy(a+plus_delta_a*delta_a, N, image_row_fft, false); 
    end
    d1 = diff(a_mat, 1);
    d2 = diff(a_mat, 2);
    d3 = diff(a_mat, 3);
    d4 = diff(a_mat, 4);
    %利用PDF中的联立的公式
    a0 = d1(1) / delta_a;
    a1 = d2(1) / (delta_a^2);
    a2 = d3(1) / (2 * delta_a^3);
    a3 = d4(1) / (6 * delta_a^4);
    A1 = 1 / a1;
    A2 = -(a2*A1^2)/a1;
    A3 = -(a2*(A1*A2+A2*A1)+a3*A1^3)/a1;
    a = a + (A1* (-a0)^1+A2*(-a0)^2+A3*(-a0)^3)
end
%显示最优a值对应的图像
imageEntropy(a, N, image_row_fft, true);


%计算特定a值下的图片熵以及是否显示图片
function E = imageEntropy(a, N, image_row_fft, display)
    n = 1:N;
    phi_mat = exp((n./N).^2 * j * a);
    image_recover = ifft(image_row_fft.* phi_mat, [], 2);
    if display
       figure();
       image(abs(image_recover));
       colormap((gray(256)));
    end
    image_square = abs(image_recover).^2;
    En = sum(image_square);
    E = - sum((image_square/En).*log2(image_square/En));
end