CN103513256B - 一种极地探冰雷达数据层位提取方法 - Google Patents

Abstract

本发明公开了一种极地探冰雷达数据层位提取方法，包括步骤：步骤1，算法初始化，输入待处理的图像，同时给出原始图像恢复、图像层位提取、规整化系数以及层位提取规整化因子的初始值；步骤2，判断迭代终止条件是否满足，若满足，则输出原始图像恢复以及图像层位提取结果，若不满足，则进入步骤3的迭代过程；步骤3，将原始图像恢复和图像层位提取分开进行，对于原始图像恢复，采用梯度框架进行处理；对于图像层位提取，采用总变分最小化结合梯度框架进行处理；在迭代过程的最后，采用共轭梯度算法更新层位提取规整化因子的值；其中步骤3每迭代一次，需要返回步骤2进行判断，直至满足迭代终止条件。利用该方法能够得到冰川的准确层位。

Description

一种极地探冰雷达数据层位提取方法

技术领域

本发明属于极地探冰雷达成像技术领域，特别涉及一种极地雷达数据层位提取方法。该方法可用于极地冰下层位提取，也可应用于探地雷达的工程探测、水文地质探测、生态环境探测，以及火星冰下探测等数据处理方面。

背景技术

极地探冰雷达利用一个发射天线发射高频宽频带的电磁脉冲，并通过一个接收天线接收来自(地下)介质层面的反射波。当电磁波在介质中传播时，其路径、电磁场强度及波形、相位等随所穿越介质的电磁特性及几何形态而变化。因此，通过检测回波时间、幅度、相位等参量，能够解算出目标深度、介质特性及结构等信息。在数据可视化上，可运用数字图像的恢复与重建技术，对冰下目标进行成像处理，以期达到对冰下目标真实和直观的显示效果。

极地探冰雷达应用于极地探测等情况下，需要对冰川的厚度和沉积过程进行探测与分析，因此，得到冰川的准确层位显得十分重要。所以，冰川层位的有效提取是目标识别和参数提取等一系列技术成功应用的根本保证。目前对于冰川层位提取采用的方法主要是从成像方面进行，如衍射层析成像、探地雷达波前成像算法和探地雷达时域反向投影成像算法等，这些方法能够得到冰川层位提取图像，但是得到的层位提取图像存在噪声，使得层位并不清晰。

发明内容

基于现有技术存在的缺陷，为了提高层位提取的有效性，本发明基于鲁棒的主分量分解-总变分原理，提出了一种极地探冰雷达数据层位提取方法，该方法能够消除得到的层位提取图像噪声并且凸显层位。

本发明的极地探冰雷达数据层位提取方法包括：步骤1，算法初始化，输入待处理的图像，同时给出原始图像恢复、图像层位提取、规整化系数以及层位提取规整化因子的初始值；步骤2，判断迭代终止条件是否满足，若满足，则输出原始图像恢复以及图像层位提取结果，若不满足，则进入步骤3的迭代过程；步骤3，将原始图像恢复和图像层位提取分开进行，对于原始图像恢复，采用梯度框架进行处理；对于图像层位提取，采用总变分最小化结合梯度框架进行处理；在迭代过程的最后，采用共轭梯度算法更新层位提取规整化因子的值；其中步骤3每迭代一次，都需要返回步骤2进行判断，直至满足迭代终止条件。

通过本发明的方法，利用雷达进行极地探冰时，能够得到冰川的准确层位。

附图说明

图1为本发明极地探冰雷达数据层位提取方法的流程图；

图2为本发明的层位提取规整化因子初始值取值流程图；

图3为本发明的原始图像恢复流程图；

图4为本发明的层位提取流程图；

图5为本发明的共轭梯度法更新层位提取规整化因子流程图；

图6为待处理的极地探冰雷达图像；

图7为本发明进行层位提取后的极地探冰雷达层位图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明提出的极地探冰雷达数据层位提取方法是基于鲁棒的主分量分解-总变分原理，该方法通过仿真实验进行了理论验证，理论分析和仿真实验结果证明了该方法的有效性。

图1为本发明极地探冰雷达数据层位提取方法的流程图。参照图1，该方法包括以下步骤。

在步骤1，对算法初始化，即输入待处理的图像，同时给出原始图像恢复、图像层位提取、规整化系数以及层位提取规整化因子的初始值。

在该步骤中，原始图像恢复和图像层位提取都是该算法最后得到的结果，需要将他们分别初始化为0。

在该步骤中，设λ为规整化系数，它的作用是在信号和噪声平滑性要求之间进行折衷。它是一个小的正数，取值范围在(0，1)，设μ为层位提取规整化因子，它的初始值采用对输入图像的信噪比进行估计得到。

图2为获取层位提取规整化因子μ的流程图。该因子的作用是克服问题病态，当数据带有误差/噪声时，保证解在真解附近连续地依赖于数据，不会把误差/噪声过分地放大。如图2所示，具体处理过程包括以下步骤：(a)通过 ${\mathrm{\μ}}_D=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}D(i+k,j+l)$ 计算输入图像的局部均值μ_D，其中(i，j)表示图像像素点的位置，D为输入的待处理图像，P和Q是计算窗口尺寸，一般选取P＝Q＝2(5×5窗)或者P＝Q＝1(3×3窗)；(b)通过 ${\mathrm{\σ}}_{\mathrm{DL}}^2(i,j)=\frac{1}{(2P+1)(2Q+1)}\underset{K=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}[D(i+k,j+l)-{\mathrm{\μ}}_D(i,j){]}^2$ 计算图像的局部方差σ_DL ²，其中D为输入的待处理图像，计算窗口P和Q尺寸的选取与局部均值计算中的选取保持一致；(c)用局部方差σ_DL ²的最大值和最小值之比作为输入图像信噪比的估计值；(d)将步骤(c)中得到的信噪比估计值进行5到10倍的缩小，得到的结果即为层位提取规整化因子μ的初始值。

在步骤2，判断迭代终止条件是否满足，若满足，则输出原始图像恢复以及图像层位提取结果。在满足迭代终止条件的情况下，如果还没有进入迭代过程，该输出结果是原始图像恢复和图像层位提取的初始值，该初始值由步骤1中初始化过程得到；如果已经经过若干次迭代过程，该输出结果是迭代过程的输出结果。若不满足迭代终止条件，则进入步骤3的迭代过程。

在该步骤中，迭代终止条件的选取定义如下：设A为原始图像恢复结果，E为图像层位提取结果，k为迭代次数，||·||_F为Frobenius范数，则迭代终止条件定义为：预设容忍值的选取需要综合考虑迭代次数和层位提取精度要求，一般选取范围在(10^-9，10^-6)。

在步骤3，将原始图像恢复和图像层位提取分开进行。对于原始图像恢复，采用梯度框架进行处理；对于图像层位提取，采用总变分最小化结合梯度框架进行处理；在每次迭代过程的最后，采用共轭梯度算法更新层位提取规整化因子的值。步骤3每迭代一次，都需要返回步骤2进行判断，直至满足迭代终止条件。

步骤3中所述的将原始图像恢复和图像层位提取分开进行的迭代处理过程具体包括以下步骤：

步骤301，设A为原始图像恢复结果，E为图像层位提取结果，t为迭代控制值，棱边保持规整化惩罚泛函其中，e是E的像元值；DE是E的支持域；规整化因子的取值为： δ是一个标定参数，用来调整惩罚泛函的离散程度。从而，A和E的梯度框架可以分别写成： $Y_k^A\←A_k+\frac{t_{k-1}-1}{t_k}(A_k-A_{k-1}),$ $Y_k^E\←E_k+J_{\mathrm{\μ}}^k\left(e\right)+\frac{t_{k-1}-1}{t_k}(E_k+J_{\mathrm{\μ}}^k\left(e\right)-E_{k-1}-J_{\mathrm{\μ}}^{k-1}\left(e\right)),$ 其中k表示迭代次数，μ表示层位提取规整化因子。

步骤302，对于A的恢复，采用基于梯度框架的最速梯度法进行。设D为输入的待处理图像。首先，其次，对G_k ^A做奇异值分解其中，S^A是奇异值矩阵，U^A和V^A是奇异值分解中产生的酉矩阵；最后，采用软阈值限制对S^A进行限制，即其中μ_k是每步迭代中的层位提取规整化因子，得到A的估计值其中μ_k为层位提取规整化因子，如图3所示。

步骤303，对于E的提取，采用基于梯度框架的总变分最小化方法来进行。首先，其次，根据计算待处理图像D的噪声方差σ_ξ ²，其中，m和n是D的尺寸，μ_ξ是D的噪声方差；得到噪声方差σ_ξ ²之后，计算噪声功率其中ξ是图像噪声；之后，采用总变分最小化方法对G_k ^E进行处理，得到结果为u_k；最后，采用最速梯度下降法恢复E，得到E的估计值其中λ为规整化系数，取值范围为(0，1)，μ_k为层位提取规整化因子，如图4所示。

步骤304，使用更新迭代控制值t的取值；使用共轭梯度法更新层位提取规整化因子μ_k的取值，具体步骤为：首先，将μ_k代入本问题的数学表达式 $F_k\left(X\right)=\mathrm{\μ}\left|\right|A_k{\left|\right|}_{*}+\mathrm{\μ\λ}|E_k{|}_1+\mathrm{\μ\λ}{J}_{\mathrm{\μ}}^k\left(e\right)+\frac{1}{2}||D-A_k-E_k{\left|\right|}_F^2$ 中求得F_k(X)，其中，A_k为第k次迭代原始图像恢复输出结果，E_k为第k次迭代层位提取输出结果，D为待处理的图像，为实现层位提取规整化的惩罚泛函，||·||_*为矩阵核范数，定义是矩阵所有奇异值之和，|·|₁为矩阵1-范数，定义是矩阵所有元素的绝对值之和。对F_k(X)做适当缩小ηF_k(X)，其中η＝0.9；其次，将上式写成 $\left(\right||A_k{\left|\right|}_{*}+\mathrm{\λ}|E_k{|}_1+\mathrm{\λ}{J}_{\mathrm{\μ}}^k\left(e\right))\mathrm{\μ}=\mathrm{\η}{F}_k\left(X\right)-\frac{1}{2}\left|\right|D-A_k-E_k{\left|\right|}_F^2;$ 最后，采用共轭梯度法搜索最佳μ值，作为下一次迭代的初始μ值，记作μ_k+1，如图5所示。

图6为输入的待处理的极地探冰雷达图像，层位不清晰，噪声较大。图7为使用本方法得到的层位提取后的极地探冰雷达层位图像，可以得到清晰的层位提取结果。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种极地探冰雷达数据层位提取方法，该方法包括步骤：

步骤1，输入待处理的图像，同时给出原始图像恢复、图像层位提取、规整化系数以及层位提取规整化因子的初始值；其中，所述原始图像恢复、图像层位提取的初始值为0，规整化系数为(0,1)之间的正数；设μ为层位提取规整化因子，它的初始值采用对输入图像的信噪比进行估计得到；

步骤2，判断迭代终止条件是否满足，若满足，则输出原始图像恢复以及图像层位提取结果，若不满足，则进入步骤3的迭代过程；

步骤3，将原始图像恢复和图像层位提取分开进行，对于原始图像恢复，采用梯度框架进行处理；对于图像层位提取，采用总变分最小化结合梯度框架进行处理；在每次迭代过程的最后，采用共轭梯度算法更新层位提取规整化因子的值；

其中步骤3每迭代一次，都需要返回步骤2进行判断，直至满足迭代终止条件；

其中，步骤1中通过以下步骤得到μ的初始值：

步骤101，通过 ${\mathrm{\μ}}_D=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}D(i+k,j+l)$ 计算输入图像的局部均值μ_D，其中，(i，j)表示图像像素点的位置，D为输入的待处理图像，P和Q是计算窗口尺寸，选取P＝Q＝2或者P＝Q＝1；

步骤102，通过 ${\mathrm{\σ}}_{\mathrm{DL}}^2(i,j)=\frac{1}{(2P+1)(2Q+1)}\underset{k=-P}{\overset{P}{\mathrm{\Σ}}}\underset{l=-Q}{\overset{Q}{\mathrm{\Σ}}}{[D(i+k,j+l)-{\mathrm{\μ}}_D(i,j)]}^2$ 计算图像的局部方差σ_DL ²，计算窗口P和Q尺寸的选取与局部均值计算中的选取保持一致；

步骤103，用局部方差的最大值和最小值之比作为输入图像信噪比的估计值；

步骤104，将步骤103中得到的信噪比估计值进行5到10倍的缩小，得到的结果即为层位提取规整化因子μ的初始值；

其中，步骤2中迭代终止条件的选取定义如下：设A为原始图像恢复结果，E为图像层位提取结果，k为迭代次数，||·||_F为Frobenius范数，则迭代终止条件定义为：A_k、A_k+1分别表示第k、k+1次迭代得到的原始图像恢复结果，E_k、E_k+1分别表示第k、k+1次迭代得到的图像层位提取结果；

其中，步骤3进一步包括：

步骤301，设t为迭代控制值，棱边保持规整化惩罚泛函其中，e是E的像元值；D_E是E的支持域；规整化因子的取值为： $D_{i,j}^x\left(e\right)=\frac{(e_{i,j+1}-e_{i,j})}{\mathrm{\δ}},$ δ是一个标定参数，用来调整惩罚泛函的离散程度；e_i,j、e_i,j+1、e_i+1,j分别表示E的第(i,j)、(i,j+1)、(i+1,j)象元值；从而，A和E的梯度框架可以分别写成： $Y_k^E\←E_k+J_{\mathrm{\μ}}^k\left(e\right)+\frac{t_{k-1}-1}{t_k}(E_k+J_{\mathrm{\μ}}^k\left(e\right)-E_{k-1}-J_{\mathrm{\μ}}^{k-1}\left(e\right));$ 分别表示第k次迭代得到的A和E的梯度框架；A_k-1表示第k-1次迭代得到的原始图像恢复结果；E_k-1表示第k-1次迭代得到的图像层位提取结果，t_k-1和t_k分别表示第k-1、k次迭代的迭代控制值，和分别表示第k、k-1次迭代的棱边保持规整化惩罚泛函；

步骤302，对于A的恢复，采用基于梯度框架的最速梯度法进行；设D为输入的待处理图像，首先，其次，对G_k ^A做奇异值分解 $(U^A,S^A,V^A)\←\mathrm{svd}\left(G_k^A\right);$ 最后，采用对S^A进行限制，得到A的估计值其中μ_k为第k次迭代的层位提取规整化因子；

步骤303，对于E的提取，采用基于梯度框架的总变分最小化方法来进行，首先， $G_k^E\←Y_k^E-\frac{1}{2}(Y_k^A+Y_k^E-D);$ 其次，根据 ${\mathrm{\σ}}_{\mathrm{\ξ}}^2=\frac{1}{\mathrm{mn}}\underset{k=0}{\overset{m-1}{\mathrm{\Σ}}}\underset{l=0}{\overset{n-1}{\mathrm{\Σ}}}{[\mathrm{\ξ}(k,l)-{\mathrm{\μ}}_{\mathrm{\ξ}}]}^2$ 计算待处理图像D的噪声方差σ_ξ ²，其中，m和n是D的尺寸，μ_ξ是D的噪声均值；得到噪声方差σ_ξ ²之后，计算噪声功率之后，采用总变分最小化方法对G_k ^E进行处理，得到结果为u_k；最后，采用最速梯度思想恢复E，得到E的估计值 $E_{k+1}=S_{\frac{\mathrm{\λ}{\mathrm{\μ}}_k}{2}}\left[u_k\right],$ 其中λ为规整化系数；

步骤304，使用更新迭代控制值t的取值；t_k+1表示第k+1次迭代的迭代控制值；并使用共轭梯度法更新层位提取规整化因子μ_k的取值，具体步骤为：首先，将μ_k代入 $F_k\left(X\right)={\mathrm{\μ}}_k{\left|\right|A_k\left|\right|}_{*}+{\mathrm{\μ}}_k\mathrm{\λ}{\left|E_k\right|}_1+{\mathrm{\μ}}_k\mathrm{\λ}{J}_{\mathrm{\μ}}^k\left(e\right)+\frac{1}{2}{\left|\right|D-A_k-E_k\left|\right|}_F^2$ 求得F_k(X)，将F_k(X)根据ηF_k(X)进行缩小，其中η＝0.9；其次，将上式写成 $({\left|\right|A_k\left|\right|}_{*}+\mathrm{\λ}{\left|E_k\right|}_1+\mathrm{\λ}{J}_{\mathrm{\μ}}^k\left(e\right)){\mathrm{\μ}}_k=\mathrm{\η}{F}_k\left(X\right)-\frac{1}{2}{\left|\right|D-A_k-E_k\left|\right|}_F^2;$ 最后，采用共轭梯度法搜索最佳μ值，作为下一次迭代的初始μ值，记作μ_k+1。

2.根据权利要求1所述的方法，其特征在于，其中预设容忍值选取范围是(10^-9,10^-6)。