CN105353375A - 一种镜像综合孔径辐射计反演方法 - Google Patents

Abstract

本发明公开了一种镜像综合孔径辐射计图像反演方法，该方法包括根据MIAS阵列输出互相关值向量R和MIAS阵列的冲击响应矩阵G获得天线输出互相关值向量R的概率测度初始值β和亮温图像不同像素点T_i的概率测度初始值α_i；计算后验概率的均值μ和后验概率的协方差Z，迭代更新β和α并判断天线输出互相关值向量R的概率测度更新值β^new和亮温T_i的概率测度更新值α^new是否同时收敛，若是则将最后一次迭代得到的μ＝βZGˊR作为亮温图像T的估计值进行输出，获得亮温图像分布T，若否，则返回进一步计算μ和Z。本发明能够有效的降低获得最优模型的计算复杂度，自动选取最优模型参数进行亮温反演，是一种新型的镜像综合孔径亮温重建方法，可以更高效地重建原始场景的亮温图像。

Description

一种镜像综合孔径辐射计反演方法

技术领域

本发明属于微波遥感及探测技术领域，更具体地，涉及一种镜像综合孔径辐射计(MIAS阵列)图像反演方法。

背景技术

镜像综合孔径辐射计就是一种由天线阵列构成的辐射计系统，所以也称为MIAS阵列。镜像综合孔径理论和镜像综合孔径微波辐射成像理论MIAS由华中科技大学的陈良兵博士提出。该理论是为了解决现有综合孔径系统的高分辨率与较大的系统复杂度之间的矛盾，采用较小的阵列加上微波反射面，从而可以在不增加天线阵元数目的情况下得到较高的余弦可见度采样频率，以此来提高系统的空间分辨率。镜像综合孔径理论提出了余弦可见度的概念，建立了亮温、余弦可见度和天线接收通道相关输出功率之间的关系，并给出了基本的亮温重建方法—反余弦变换，其中余弦可见度是指对亮温进行余弦变换所得物理量。此外，还给出了一种解决基线缺失的方法，即通过平移微波反射面，将多个位置处获得的采样频率进行组合，该方法在仿真中效果显著。

镜像综合孔径微波辐射成像方法结合了传统的二元干涉仪和海面Lloyd镜干涉，Lloyd镜干涉是物理光学中一个专有概念，表示通过海面反射信号与信号自身具有相关性的概念，两种方法的优势，微波反射面改变了目标场景辐射信号的传递路径，使得通道的相关输出也发生了变化，镜像综合孔径系统的两个通道之间只需进行同相相关，不需正交相关，即系统不需要IQ(正交相关)解调，如此一来便大大降低了信号处理的复杂度。

镜像综合孔径系统用于成像，同样面临着亮温图像重建的问题，同传统的综合孔径亮温图像重建一样，它也是一个病态的反问题，即由于误差等非理想因素的存在而导致问题的解可能不存在、不唯一或者不稳定。现有的正则化方法用于亮温重建，存在模型参数选取困难问题，反演效率难以进一步提高。

发明内容

本发明提供了一种镜像综合孔径辐射计图像反演方法，旨在解决现有的镜像综合孔径反演方法采用非最优的正则化参数导致反演结果的性能降低的问题。

为实现上述目的，本发明提供了一种镜像综合孔径辐射计图像反演方法，包括下述步骤：

一种镜像综合孔径辐射计图像反演方法，其特征在于，包括下述步骤：

(1)对于任意平面MIAS阵列，求出MIAS阵列输出互相关值向量R和天线阵列冲击响应矩阵G；其中，R的维度为该平面天线阵列的空间采样频率个数M，R的方差为G的行数为M，列数等于测量的亮温图像T的像素点数N；测量的亮温图像T各像素点i的亮温的方差为

(2)迭代计算R的概率测度更新值β_new和各像素点i处亮温的概率测度值α＝[α₁,α₂,…，α_N]；用β表示MIAS阵列每个基线对应的互相关值的方差的倒数，赋其初值为用α_i表示像素点i的亮温T_i的方差的倒数，赋其初值为

(2.1)分别计算方差矩阵Z和亮温后验概率的均值μ：

Z＝(A+βG'G)^-1，(I)

μ＝βZG'R，(II)

其中A＝diag(α)为对角阵；Diag(*)是Matlab中用于构造一个对角矩阵的函数，其输入参数为向量，返回一个方阵，其对角线元素为输入向量，不在对角线上元素全为0。

(2.2)令γ_i＝1-α_iZ_ii，将(2.1)计算出的Z和μ代入下式进行迭代计算，

${\mathrm{\α}}_{i,new}=\frac{{\mathrm{\γ}}_i}{{\mathrm{\μ}}_i^2},---\left(III\right)$

${\mathrm{\β}}_{new}=\frac{M-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\γ}}_i}{\left|\right|R-G\mathrm{\μ}|{|}^2},---\left(IV\right)$

其中，Z_ii表示矩阵Z的第i个对角元；||*||表示对*取F-范数；

(2.3)判断所述天线输出互相关值的概率测度更新值β_new和亮温T的概率测度更新值α_i，new是否同时收敛，是则转步骤(3)；否则转子步骤(2.4)；

(2.4)判定迭代次数是否小于预定的迭代上限次数K，是则将β、α_i值分别更新为β_new、α_i,new，转子步骤(2.1)继续迭代计算；否则判定不收敛，迭代结束，算法中止，将最后一次迭代得到β_new、α_i,new作为迭代终值；

(3)根据迭代终值β_new以及α_i，new，按式(I)、(II)计算亮温后验概率的均值μ，即亮温分布估计值T。

步骤(1)所述的采样频率，是由阵列中天线的相对位置决定的，这里假设天线1的位置为(x₁,y₁)，天线2的位置为(x₂,y₂)，镜像综合孔径系统工作的中心频率对应的波长为λ，则上述采样频率定义为那么对于一维镜像综合孔径阵列而言，每一对天线可以形成两个采样频率，对二维而言，则可以形成四个采样频率；

进一步的，所述的反演方法，其特征在于，在步骤(1)之前，还包括下述步骤：

以镜像综合孔径辐射计系统工作频率半波长的整数倍为单位，平移阵列反射面；对一维阵列，沿阵列垂直的方向进行平移；对二维阵列，反射面在水平和垂直两个维度上分别进行移动；阵列反射面移动次数，以在所形成的采样频率均匀布满最小频率与最大频率之间为原则。这样可以进一步弥补缺失基线，提升图像重建质量。均匀布满，就是如果最小采样频率为f，最大采样频率为N*f，N为正整数，那么不移动，在1～N之间就会有些采样频率是没有的，称为缺失，会引入不可弥补的成像误差；通过移

动，就可以保证1～N之间的所有采样频点都得到，称为均匀布满。

进一步的，所述的反演方法，其特征在于，在子步骤(2.3)中，通过判别是否|β_new-β|²≤e和||α_new-α||²≤e，确定β_new和α_new是否同时收敛，其中e为迭代计算精度控制量，大小在10^-5～10^-3之间，具体值根据收敛速度和计算精度需求进行权衡。

进一步的，所述的反演方法，其特征在于，在子步骤(2.4)中迭代上限次数K根据需要设定，取值范围是100-300。

本发明提出的方法，利用了贝叶斯模型和公式，也可称为贝叶斯统计反演方法。本发明提出反演方法是基于如下原理的：

设从亮温到可见度传递过程中的所有误差之和为ε，可以将其看做总体系统噪声，则镜像综合孔径阵列的系统接收模型可描述为：R＝GT+ε，由于噪声ε_i近似满足均值为0，方差为σ²的高斯分布，且各通道噪声之间满足独立同分布的条件，其联合概率密度函数为：p(ε)＝N(0，σ²I)，由此可得R的条件概率分布p(R|T)服从正态分布N(GT,σ²I)，其中I表示单位矩阵，维度与R的元素个数相等，得出p(R|T)表达式为：

$p\left(R\right|T)={\left(\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\right)}^M{e}^{-\frac{\left|\right|R-GT|{|}^2}{2{\mathrm{\σ}}^2}}---\left(1\right)$

M表示向量R的元素个数；σ为向量R中各元素的方差，这里设R中各元素的方差是相同的；

依据亮温T_i的概率分布p(T_i|α_i)服从正态分布且相互独立，得到T的联合概率分布为：

$p\left(T\right|\mathrm{\α})=\underset{i=1}{\overset{N}{\Π}}p\left(T_i\right|{\mathrm{\α}}_i^{-1})={\left(\frac{1}{\sqrt{2\mathrm{\π}}}\right)}^N\underset{i=1}{\overset{N}{\Π}}\sqrt{{\mathrm{\α}}_i}{e}^{-\frac{{\mathrm{\α}}_i{T}_i^2}{2}}---\left(2\right),$

式中T_i表示像素点i处的亮温，α_i为不同像素点i处的亮温的概率测度值，代表亮温的方差；α＝[α₁,α₂,…，α_N]，N表示亮温的像素点个数；

根据贝叶斯公式得到p(T|R,α,β)(其推导过程参考：MendelsonS,SmolaAJ.AdvancedLecturesonMachineLearning[M].SpringerBerlinHeidelberg,2003)

$p\left(T\right|R,\mathrm{\α},\mathrm{\β})={\left(\frac{1}{\sqrt{2\mathrm{\π}}}\right)}^M\frac{1}{\sqrt{\left|Z\right|}}{e}^{-\frac{{(T-\mathrm{\μ})}^{,}{Z}^{-1}(T-\mathrm{\μ})}{2}}---\left(3\right),$

上式中，μ为亮温后验概率的均值，方差Z＝(A+βG'G)^-1，μ＝βZG'R；β＝σ^-2是天线输出互相关值的概率测度初始值，A＝diag(α)是由α构成的对角阵；由公式(2)和(3)可以得到天线输出互相关向量R关于超参数α和σ²的边缘似然函数：

$p\left(R\right|\mathrm{\α},{\mathrm{\σ}}^2)={\left(\frac{1}{\sqrt{2\mathrm{\π}}}\right)}^M\frac{1}{\sqrt{|{\mathrm{\σ}}^{2}I+{\mathrm{GA}}^{-1}{G}^{,}|}}{e}^{-\frac{R^{,}{({\mathrm{\σ}}^{2}I+{\mathrm{GA}}^{-1}{G}^{,})}^{-1}R}{2}}---\left(4\right),$

公式(4)给出了各天线互相关值之间的联合分布函数，将公式(4)作为贝叶斯统计反演模型中的证据函数，用以推导待估计参数的表达式。对公式(4)两边取对数，求导取极值，得到α_i以及σ²的迭代公式(具体推导过程参见：杨宏.基于机器学习的SAIR图像反演和目标检测方法研究[D].华中科技大学,2012.)：

${\mathrm{\α}}_{i,new}=\frac{{\mathrm{\γ}}_{i,old}}{{\mathrm{\μ}}_{i,old}^2}---\left(5\right),$

${\mathrm{\β}}_{new}=\frac{M-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\γ}}_{i,old}}{\left|\right|R-{\mathrm{G\μ}}_{old}|{|}^2}---\left(6\right),$

迭代第一步，将所述MIAS阵列输出互相关值矢量R的概率测度初始值β₁、不同像素点i处亮温的概率测度初始值α_i,1和所述亮温后验概率的均值μ以及方差Z代入迭代公式(5)、(6)，经过迭代运算，可得到天线输出互相关值矢量R的概率测度更新值β_new和不同像素点i处亮温的概率测度更新值α_i，new＝[α_1,new,α_2new,…，α_N，new]；γ_i＝1-α_iZ_ii，Z_ii表示矩阵Z的第i个对角元，||*||表示对*取F-范数；如果β_new和α_i，new同时收敛，则按式Z＝(A+βG'G)^-1，μ＝βZG'R，用β_new和T_i的迭代终值计算出亮温后验概率的均值μ，即亮温分布估计值T。

本发明从贝叶斯模型选取的角度阐明镜像综合孔径辐射计亮温反演问题；通过对亮温以及天线输出互相关值作统计假设，得出亮温到天线互相关输出的统计传播特性。根据贝叶斯原理推导出亮温的后验概率公式，结合极大似然估计方法，提出迭代计算公式，用于估计反演参数。本发明能够有效的降低获得最优模型的计算复杂度，自动选取最优模型参数进行亮温反演，是一种新型的镜像综合孔径亮温重建方法。

附图说明

图1是本发明实施例提供的一种镜像综合孔径辐射计图像反演方法的实现流程图；

图2是误差因素对本方法反演及反余弦变换反演算法误差均方根的影响结果图；

图3是飓风场景下Tikhonov正则化反演结果和贝叶斯统计反演结果对比图；

图4是一维镜像综合孔径阵列与反射面组合示意图。

图5是具体实施方式中的一维镜像综合孔径阵列反演亮温图像结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例提供的镜像综合孔径微波辐射计图像反演方法基于贝叶斯回归估计，针对镜像综合孔径成像的实时应用，能够有效的降低其计算复杂度，是一种可以自动的选取最优模型的镜像综合孔径图像反演方法。本实施例以一个16单元一维镜像综合孔径阵列MIAS为例，包括一个一维最小冗余线阵加一个镜反射面结构，阵列单元接收机噪声温度T_rec＝500k，系统带宽B＝200MHz，积分时间τ＝1秒，亮温图像平均亮温T_A＝250k，这些数据对R、G，以及β和α_i的初值和初值有影响，但是R和G又不是仅仅依赖于这些参数，主要还是由阵列和被测场景亮温T决定。

天线阵列以及镜反射面位置组合示意图如图4所示。反演过程如下(如图1所示)：

第一步，输入R、G；R、G是已知的输入参数，可根据实际的天线阵列参数采用通用的方法计算得到，R、G也可以通过以下步骤取得：

1.输入阵列参数：天线数目、阵元最小间距、阵列类型、阵列空间采样频率；2.输入反射面参数：反射面起始位置、反射面位置组合；3.输入MIAS系统参数：工作频率、接收机带宽、噪声系数、积分时间；通过阵列信号互相关得到天线阵列输出互相关值向量R；通过外部点响应测量或者天线阵列内部分系统测量得到冲击响应矩阵G，并将每个反射面位置处获得的G矩阵进行组合，获取组合后的G矩阵；

第二步，根据MIAS阵列参数T_rec、T_A、B和τ给出R的概率测度初始估计值待反演的亮温图像T的概率测度赋予初始估计值本实施例中迭代计算精度控制量e取值为10^-4，迭代上限次数K取值为200；

第三步，根据权利要求书中公式(I)(II)(III)(IV)迭代计算β和α_i的更新值β_new和α_i，new，本实施例中，β_new和α_i，new同时收敛，迭代过程终止。

第四步，根据迭代过程结束时的β_new以及α_i，new和权利要求书中的公式(II)计算亮温后验概率的均值μ，得到要求的亮温分布值T，如图5所示。

为了更进一步的说明本发明性能，现通过仿真对比实验中的反余弦变换反演法和本实施例所用的统计反演方法，仿真参数如表1所示：

表1

图2显示了系统误差因素对本方法反演及现有的反余弦变换反演算法误差均方根的影响仿真结果对比，可见本发明提出的方法能够有效的减少反演图像和真实场景的偏差，受误差因素影响更小，因此本方法反演性能优于反余弦变换反演法。

为了进一步通过仿真对比分析传统数值化正则化反演方法同本方法的性能。我们将对本方法、Tikhonov正则化反演法的反演性能进行分析。图3显示了飓风场景的Tikhonov正则化反演和本方法反演仿真结果对比，我们可以看出本方法反演可以取得同传统数值化正则化反演方的最优化的解性能相同的反演结果。不同于传统Tikhonov数值化正则化反演方法需要人为的获得最优解，本反演方法可以自动的获得最优值，得到偏差和方差最小的反演图像，具有更快的求解速度和无需人工干预的优点。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种镜像综合孔径辐射计图像反演方法，其特征在于，包括下述步骤：

(1)对于任意平面MIAS阵列，求出MIAS阵列输出互相关值向量R和天线阵列冲击响应矩阵G；其中，R的维度为该平面天线阵列的空间采样频率个数M，R的方差为G的行数为M，列数等于测量的亮温图像T的像素点数N；测量的亮温图像T各像素点i的亮温方差为

(2)迭代计算R的概率测度更新值β_new和各像素点i处亮温的概率测度值α＝[α₁,α₂,…，α_N]；用β表示MIAS阵列每个基线对应的互相关值的方差的倒数，赋其初值为用α_i表示像素点i的亮温T_i的方差的倒数，赋其初值为

(2.1)分别计算方差矩阵Z和亮温后验概率的均值μ：

Z＝(A+βG'G)^-1，(I)

μ＝βZG’R,(II)

其中A＝diag(α)为对角阵；

(2.2)令γ_i＝1-α_iZ_ii，将(2.1)计算出的Z和μ代入下式进行迭代计算，

${\mathrm{\α}}_{i,new}=\frac{{\mathrm{\γ}}_i}{{\mathrm{\μ}}_i^2},---\left(III\right)$

${\mathrm{\β}}_{new}=\frac{M-\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\γ}}_i}{\left|\right|R-G\mathrm{\μ}|{|}^2},---\left(IV\right)$

其中，Z_ii表示矩阵Z的第i个对角元；||*||表示对*取F-范数；

(2.3)判断所述天线输出互相关值的概率测度更新值β_new和亮温T的概率测度更新值α_i，new是否同时收敛，是则转步骤(3)；否则转子步骤(2.4)；

(2.4)判定迭代次数是否小于预定的迭代上限次数K，是则将β、α_i值分别更新为β_new、α_i,new，转子步骤(2.1)继续迭代计算；否则判定不收敛，迭代结束，算法中止，将最后一次迭代得到β_new、α_i,new作为迭代终值；

(3)根据迭代终值β_new以及α_i，new，按式(I)、(II)计算亮温后验概率的均值μ，即亮温分布估计值T。

2.根据权利要求1所述的反演方法，其特征在于，在步骤(1)之前，还包括下述步骤：

以镜像综合孔径辐射计系统工作半波长的整数倍为单位，平移阵列反射面；对一维阵列，沿阵列轴线方向进行平移；对二维阵列，反射面在水平和垂直两个方向上分别进行移动；阵列反射面移动次数，以所形成的采样频率均匀布满最小频率与最大频率之间为原则。

3.如权利要求1所述的反演方法，其特征在于，在子步骤(2.3)中，通过判别是否|β_new-β|²≤e和||α_new-α||²≤e，确定β_new和α_new是否同时收敛，其中e为迭代计算精度控制量，大小在10^-5～10^-3之间，具体值根据收敛速度和计算精度需求进行权衡。

4.如权利要求1所述的反演方法，其特征在于，在子步骤(2.4)中迭代上限次数K根据需要设定，取值范围是100-300。