基于贝叶斯准则的超分辨算法原理介绍

贝叶斯准则

阵列方位向输出结果可以表示为
$$y=Hx\text{+}n$$
式中：$H\in R^{L\times L}$表示卷积核循环矩阵，$y\in R^{L\times 1}$表示观测结果，$x\in R^{L\times 1}$表示目标散射信息，$n\in R^{L\times 1}$表示引入的加性噪声。由贝叶斯理论可知，在给定输出结果$y$的条件下，$x$的后验概率可以表示为
$$p(x|y)=\frac{p(y|x)p(x)}{p(y)}$$
式中：$p(x|y)$表示给定$y$的条件下$x$的后验概率，$p(x)$表示目标散射信息的先验概率分布，$p\left(y|x\right)$表示似然函数分布，$p(y)$为接收信息的先验概率分布。

最大后验概率估计的特点就是引入了原始信号的先验概率分布函数，可以获得关于$x$的最大后验概率分布函数$p(x|y)$，可以表示为
$$f=\arg \max p(x|y)\text{=arg max}\frac{p(y|x)p(x)}{p(y)}$$
上式中$p(y)$不是关于$x$的函数，其可以表示为
$$f=\arg \max p(y|x)p(x)$$
上述方法称为最大后验估计方法。一般情况下在无法获取$x$的先验信息时，通常认为$p\left(x\right)$是均匀分布，在这种情况下通过上式可以获得关于$x$的最大似然函数$p(y|x)$，这种方法称为最大似然估计方法。

最大似然准则

最大似然准则
$$y=Possion(Hx)+Possion(b)+Gaussian\left(n\right)$$
式中：$b$表示原始场景目标的环境噪声，$n$表示接收机存在的加性高斯噪声，根据信号的先验信息来选择合适的取值。接收信号的期望值是不含有高斯噪声的，可以表示为
$$\widehat{y_i}=\sum _{j=1}^L{H}_{ij}{x}_j+b_i$$
接收信号的各点元素的似然分布函数可以表示为
$$p(y_i|x)=\sum _{N=0}^{+\infty }\frac{e^{-{(Hx+b)}_i}(Hx+b{)}_i^N}{N!}\frac{1}{\sqrt{2\pi {\epsilon }^2}}{e}^{-\frac{{(x_i-N)}^2}{2{\epsilon }^2}}$$
设接收信号各点元素都相互独立，那么混合分布的联合概率密度函数可以表示为
$$p(y|x)=\prod _{i=1}^L\left(\sum _{N=0}^{+\infty }\frac{e^{-{(Hx+b)}_i}(Hx+b{)}_i^N}{N!}\frac{1}{\sqrt{2\pi {\epsilon }^2}}{e}^{-\frac{{(x_i-N)}^2}{2{\epsilon }^2}}\right)$$

为了实际中运算的方便，对上式取对数并且取反可以得到
$$-\ln \left[p(y|x)\right]\propto J(x)=\sum _{i=1}^L\left[{(Hx+b)}_i-\ln \left(\sum _{N=0}^{+\infty }\frac{(Hx+b{)}_i^N}{N!}\frac{1}{\sqrt{2\pi {\epsilon }^2}}{e}^{-\frac{{(x_i-N)}^2}{2{\epsilon }^2}}\right)\right]$$

求取原函数的最大值等价于求$J(x)$的最小值，对上式进行关于$x$求导可以得到
$$\frac{\partial \ln J(x)}{\partial x_i}=\sum _{i=1}^L{H}_{ij}-\sum _{i=1}^L{H}_{ij}\frac{{\varphi }_i(x)}{{\phi }_i(x)}=0$$
式中：
$${\varphi }_i(x)=\varphi ({(Hx+b)}_i,y_i)\text{=}\sum _{k=0}^{+\infty }\frac{(Hx+b{)}_i^k}{k!}{e}^{-\frac{{\left(k-y_i\right)}^2}{2{\epsilon }^2}}$$
$${\phi }_i(x)=\phi ({(Hx+b)}_i,y_i)\text{=}\sum _{k=0}^{+\infty }\frac{(Hx+b{)}_i^k}{k!}{e}^{-\frac{{\left(k\text{+}1-y_i\right)}^2}{2{\epsilon }^2}}$$
重新定义一个向量$h$，其中$h_j=\sum _{i=i}^L{H}_{ij}$，则上式可以表示为
$$\sum _{i=1}^L{H}_{ij}\frac{{\phi }_i(x)}{{\varphi }_i(x)}=h_j$$
进一步地
$$H^T\left(\frac{\phi (x)}{\varphi (x)}\right)=h$$
采用Picard迭代方法对上式进行求解可以得到
$$x^{k+1}=H^T\left(\frac{\phi (x^k)}{\varphi (x^k)}\right)\frac{x^k}{h}$$

式中：$k$表示迭代的次数，一般情况下令初始值$x^0=y$，同时对波束方向图进行归一化处理，即$h_j=\sum _{i=i}^L{H}_{ij}=1$。

由于上式属于非线性迭代过程，在每次迭代中都会产生新的频率分量，实现频谱的外推，进而实现角度超分辨。考虑到实际情况，$\frac{\phi (x)}{\varphi (x)}$严格为正，同时卷积核矩阵非零，因此若迭代过程中$x^0>0$，则每次迭代的结果都有$x^k>0$，所以重写可以得到
$$x^{k+1}=x^k+\frac{x^k}{h}\left(H^T\frac{\phi (x^k)}{\varphi (x^k)}-h\right)$$

对于上述迭代过程直接求解计算量较大，不利于实际的工程应用。由统计学知识可以得到，在一定的条件下，高斯分布和泊松分布是可以相互转化的，因此有两种方法：将高斯分布转化为泊松分布来求解；将泊松分布转化为高斯分布来求解。

对于将高斯分布转化为泊松分布求解来说，将方差为${\epsilon }^2$的高斯白噪声带入，对于第$i$个元素可以得到
$$y_i+{\epsilon }^2={(Hx+b)}_i+(n_i+{\epsilon }^2)$$
则迭代过程可以写作为
$$x^{k+1}=x^k\left[H^T\left(\frac{y+{\epsilon }^2}{H{x}^k+b+{\epsilon }^2}\right)\right]$$
对于将泊松分布转化为高斯分布来求解来说
$$\frac{{\phi }_i(x)}{{\varphi }_i(x)}=\exp \left[-\frac{1+2\left({(Hx+b)}_i-y_i\right)}{2\left({(Hx+b)}_i+{\epsilon }^2\right)}\right]\left[1+o(\frac{1}{\sqrt{y_i}})\right]$$

所以迭代的过程可以写作为
$$x^{k+1}=x^k{H}^T\exp \left[-\frac{1+2\left((H{x}^k+b)-y\right)}{2\left((H{x}^k+b)+{\epsilon }^2\right)}\right]$$
对上式中的指数项进行展开可以得到
$$\exp \left[-\frac{1+2\left((H{x}^k+b)-y\right)}{2\left((H{x}^k+b)+{\epsilon }^2\right)}\right]=1-\frac{1+2\left((H{x}^k+b)-y\right)}{2\left((H{x}^k+b)+{\epsilon }^2\right)}=\frac{y+{\epsilon }^2}{H{x}^k+b+{\epsilon }^2}$$

原始场景

接收场景

恢复场景

恢复的场景较接收的场景有了很大的改善。

最大后验概率准则

由前节叙述可知，关于$x$的最大后验估计为
$$f=\arg \max p(y|x)p(x)$$
若回波信号中的噪声服从复高斯分布，则对于每个采样点的似然函数可以表示为
$$p\left(y_n|x_n\right)=\frac{1}{\pi {\sigma }^2}{e}^{-\frac{1}{{\sigma }^2}{\left(y_n-H{x}_n\right)}^2}$$
式中：${\sigma }^2$表示复高斯分布函数的方差，若各回波采样点的噪声独立同分布，则联合似然分布函数为
$$p\left(y|x\right)=\frac{1}{\pi {\sigma }^{2MN}}{e}^{-\frac{1}{{\sigma }^2}{\Vert y-Hx\Vert }_2^2}$$
式中：$M,N$表示回波矩阵的行和列。如果用拉普拉斯分布表示目标的分布，则
$$p\left(x\right)=\prod _{m=1}^{MN}\frac{1}{2\mu }{e}^{\left(-\frac{|x_{{}_m}-b|}{\mu }\right)}$$
式中：$\mu >0$是拉普拉斯分布的尺度参数，$b=0$表示位置参数。那么，最大后验概率分布函数为
$$\begin{array}{rl}& f=p(x)p(y|x)\\ & =\prod _{m=1}^{MN}\frac{1}{2\mu }{e}^{\left(-\frac{|x_{{}_m}|}{\mu }\right)}\cdot \frac{1}{\pi {\sigma }^{2MN}}{e}^{-\frac{1}{{\sigma }^2}\Vert y-Hx\Vert {}_2^2}\end{array}$$
取对数并取反得
$$H\left(x\right)=-\ln f=\frac{1}{\mu }{\Vert x\Vert }_1+\frac{1}{{\sigma }^2}{\Vert y-Hx\Vert }_2^2-MN\ln \pi {\sigma }^2$$
式中：$1/\mu$表示正则化参数，由于${\Vert x\Vert }_1$在零点不可微，引入接近于0的非负实数$\epsilon$，则上式可以写为
$$H\left(x\right)\approx \frac{1}{\mu }\sum _{m=1}^{MN}\sqrt{{\left|x_m\right|}^2+\epsilon }+\frac{1}{{\sigma }^2}{\Vert y-Hx\Vert }_2^2-MN\ln \pi {\sigma }^2$$
对上式进行梯度运算可以得到
$$\nabla \left(H\left(x\right)\right)=H^{H}Hx-H^{H}y+\frac{{\sigma }^2}{\mu }W\left(x\right)x$$

式中：${\left(\bullet \right)}^H$表示矩阵的共轭，$W\left(x\right)=diag\left\{{\left({\left|x_m\right|}^2+\epsilon \right)}^{-\frac{1}{2}}\right\}$ 为对角矩阵。由于上式是关于$x$的非线性函数，无法利用$\nabla \left(H\left(x\right)\right)\text{=}0$直接求解，所以采用迭代的方法对上述求解，先得到其直接解为
$$x={\left(H^{H}H+\frac{{\sigma }^2}{\mu }W\left(x\right)\right)}^{-1}{H}^{H}y$$

采用最小二乘解$x={\left(H^{H}H\right)}^{-1}{H}^{H}y$作为迭代的初值，则迭代过程可以写作为
$$x_{k+1}={\left(H^{H}H+\frac{{\sigma }^2}{\mu }W\left(x_k\right)\right)}^{-1}{H}^{H}y$$

考虑到每次迭代噪声会对最终的结果造成影响，所以将每次迭代的噪声进行更新，即
$${{\sigma }_k}^2=\frac{1}{MN}{\Vert y-H{x}_k\Vert }_2^2$$
那么，最终的迭代过程为
$$x_{k+1}={\left(H^{H}H+\frac{{{\sigma }_k}^2}{\mu }W\left(x_k\right)\right)}^{-1}{H}^{H}y$$

接收场景

恢复场景

可以看出，先验信息的加入能够使得算法取得更好的效果。