CN105259539A - 一种针对极地探冰雷达数据的分部式噪声抑制方法 - Google Patents

Abstract

本发明涉及极地探冰雷达数据的噪声抑制方法，具体涉及一种针对极地探冰雷达数据的分部式噪声抑制方法。包括步骤为：对待处理探冰雷达数据进行预滤波处理；基于滤波后各道雷达数据的幅度和相位，对待处理冰雷达数据进行分部划分；对各分部冰雷达数据，分别进行最优低秩滤波处理；对各分部滤波处理结果，进行累加求和，得到噪声抑制后的探冰雷达数据。本发明基于极地探冰雷达数据的幅度和相位信息，对冰雷达数据的不同介质层位进行分部式滤波处理，通过对不同介质层位的分部式噪声抑制处理，可有效反映冰下微弱信号的变化，去噪同时保留冰下微弱层位特征。

Description

一种针对极地探冰雷达数据的分部式噪声抑制方法

技术领域

本发明涉及极地探冰雷达数据的噪声抑制方法，具体涉及一种针对极地探冰雷达数据的分部式噪声抑制方法。

背景技术

探冰雷达系统凭借探测深度深、分辨率高和数据信息丰富等优点，在各国的冰雪项目中发挥了重要作用。空气和冰盖交界面、冰盖内部层、冰下冰岩交界面、冰下湖是冰雷达探测的主要目标。由于探冰雷达在数据采集过程中会引入噪声干扰，通常噪声和干扰远强于深部目标的反射信号，难以提取深部层位特征。本发明主要涉及探冰雷达数据的随机噪声抑制方法。

近年来，针对随机噪声的抑制问题，人们提出了很多有效的噪声抑制方法，例如平均滤波、基于信号变换域特性的滤波(小波变换、曲波变换等)、基于图像局部平滑性的噪声抑制方法(总变分最小化、基于偏微分方程的方法等)、基于数据稀疏性的去噪方法(压缩感知、字典学习方法等)等。这些噪声抑制方法，多是对数据的全局性去噪处理，没有考虑含噪数据中不同介质的属性差异；虽然提高了图像的信噪比，但在某些情况下会丢失深部弱反射目标。

因此，如何在受噪声干扰的探冰雷达数据中，去噪同时有效保持层位边缘特征，尤其是深部微弱层位特征，对于探冰雷达成功用于极地冰川探测具有重要的意义。

发明内容

本发明的目的在于针对探冰雷达数据的不同介质属性，进行分部式滤波处理，实现噪声抑制的同时，保留层位边缘特征，即在提高信噪比的情况下，保证图像的空间分辨率。

本发明是通过以下技术方案实现的，包括以下步骤：

第一步，对待处理冰雷达数据X(r,t)进行预滤波处理，得到预处理后的数据

所述的对待处理冰雷达数据X(r,t)进行预滤波处理，具体是：对数据X(r,t)采用奇异值分解方法，进行低秩滤波处理。

第二步，对计算各道冰雷达数据的幅度和相位，基于此对待处理冰雷达数据X(r,t)进行分部划分，得到K个分部冰雷达数据X^(k)(r,t),(k＝1,2,...,K)。

所述的计算各道冰雷达数据的幅度和相位，基于此对待处理冰雷达数据X(r,t)进行分部划分，具体是：对滤波后各道冰雷达数据，计算得到幅度序列和相位序列；绘制幅度序列和相位序列的包络，设定局部极值点判定阈值，得到各包络变化局部极大值；结合幅度序列和相位序列包络变化局部极大值分布，定义不同介质的分部参考位置，对待处理冰雷达数据X(r,t)进行分部划分，得到K个分部冰雷达数据X^(k)(r,t),(k＝1,2,...,K)。

第三步，对各分部冰雷达数据X^(k)(r,t),(k＝1,2,...,K)，分别进行最优低秩滤波处理，得到各分部滤波处理结果(k＝1,2,...,K)。

所述的低秩滤波处理，具体是：对X^(k)(r,t),(k＝1,2,...,K)，分别构造其Casorati矩阵如下：

其中，代表在空间位置r_i和采样时间t_j的空间-时间采样值，N和M分别为第k个分部观测数据的空间和时间长度；对矩阵C^(k)，采用奇异值分解方法进行最优低秩估计；假设C^(k)的秩为L^(k)，则其低秩估计结果如下：

$C^{\left(k\right)*}=\underset{l=1}{\overset{L^{\left(k\right)}}{\Σ}}{\mathrm{\σ}}_l^{\left(k\right)}{u}_l^{\left(k\right)}{v}_l^{\left(k\right)H},$

其中，和分别是矩阵C^(k)的第l个奇异值、第l个左奇异向量和第l个右奇异向量，所对应的矩阵表示形式分别为Σ^(k)，U^(k)和V^(k)，H为取矩阵的共轭转置；L^(k)为C^(k)的秩，可由X^(k)(r,t)中的噪声方差||ζ^(k)||₂通过下式估计获得

${\mathrm{\λ}}_{\hat{L}+1}^{\left(k\right)}\left(C^{\left(k\right)}\right)\≤\left|\right|{\mathrm{\ζ}}^{\left(k\right)}|{|}_2\≤{\mathrm{\λ}}_{\hat{L}}^{\left(k\right)}\left(C^{\left(k\right)}\right),$

其中，为C^(k)的第个特征值；由此，各分部数据的分部滤波处理结果为 ${\hat{X}}^{\left(k\right)}(r,t)=U^{\left(k\right)}{\mathrm{\Σ}}^{\left(k\right)}{V}^{\left(k\right)}H,$ (k＝1,2,...,K)。

第四步，将各分部滤波结果累加求和，得到噪声抑制后的处理结果为：

$\hat{X}(r,t)=\underset{k=1}{\overset{K}{\Σ}}{\hat{X}}^{\left(k\right)}(r,t).$

与现有技术相比，本发明的有益效果是：

(1)本发明一种针对极地探冰雷达数据的分部式噪声抑制方法，特别是采用预滤波处理技术，可有助于实现对不同介质层位的分部划分；

(2)基于探冰雷达数据的幅度和相位信息，可有效反映出冰下不同介质属性的差异，实现对冰下介质层位的分部划分；

(3)对不同介质层位进行不同程度的分部式滤波处理，可以对不同反射强度的层位特征给予不同的处理，有助于保留深部弱反射目标，去噪同时有效保留深部微弱层位结构。

附图说明

图1为本发明极地探冰雷达数据噪声抑制方法的流程图。

图2为本发明实施例中待处理的极地探冰雷达数据成像图。

图3为发明本实施例中待处理图像第900道的幅度和相位序列信号波形图。

图4为发明本实施例中待处理图像第900道的幅度和相位序列包络图。

图5为对本发明实施例进行噪声抑制后的探冰雷达数据成像图。

图5中，(a)为采用本发明方法的噪声抑制结果，(b)为采用小波变换的噪声抑制结果。

图6为对本发明实施例进行噪声抑制后第900道的信号波形图。

图6中，(a)为采用本发明方法的噪声抑制结果，(b)为采用小波变换的噪声抑制结果。

图7为对本发明实施例进行噪声抑制后第900道的灰度梯度图。

图7中，(a)为采用本发明方法的噪声抑制结果，(b)为采用小波变换的噪声抑制结果。

具体实施方式：

下面结合附图对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和过程，但本发明的保护范围不限于下述的实施例。

本发明中所使用的数据来源于中国科学院电子学研究所研制的探冰雷达，没有具体型号，没有批产。具体参数如下：工作带宽：100MHz，中心频率：150MHz，发射信号类型：线性调频信号，天线阵列：2个。

如图1所示，本发明极地探冰雷达数据噪声抑制方法实施例的流程如下：

第一步，对图2所示待处理冰雷达数据X(r,t)(图像大小为1365×20480)进行预滤波处理，得到预处理后的数据

(1)对X(r,t)，采用奇异值分解方法，取前L＝20个奇异值，进行低秩滤波处理。

第二步，对计算各道冰雷达数据的幅度和相位，基于此对待处理冰雷达数据X(r,t)进行分部划分，得到K个分部冰雷达数据X^(k)(r,t),(k＝1,2,...,K)。

(1)对各道冰雷达数据，计算得到幅度序列和相位序列，例如第900道数据的幅度序列和相位序列，如图3所示。

(2)绘制幅度序列和相位序列的包络，设定局部极值点判定阈值，得到各包络变化局部极大值，如图4所示。

(3)结合幅度序列和相位序列包络变化局部极大值分布，定义不同介质的分部参考位置，本实施例中，将待处理冰雷达数据划分为地表、冰盖表层、冰盖内部层、冰盖底部、冰岩交界面和冰岩底部，可得6个分部冰雷达数据X^(k)(r,t),(k＝1,2,...,6)。

第三步，对各分部冰雷达数据X^(k)(r,t),(k＝1,2,...,6)，分别进行最优低秩滤波处理，得到各分部滤波处理结果(k＝1,2,...,6)。

(1)对X^(k)(r,t),(k＝1,2,...,6)，分别构造其Casorati矩阵如下：

其中，代表在空间位置r_i和采样时间t_j的空间-时间采样值，N和M分别为第k个分部观测数据的空间和时间长度。

(2)对矩阵C^(k)，采用奇异值分解方法进行最优低秩估计：

$C^{\left(k\right)*}=\underset{l=1}{\overset{L^{\left(k\right)}}{\Σ}}{\mathrm{\σ}}_l^{\left(k\right)}{u}_l^{\left(k\right)}{v}_l^{\left(k\right)H},$

其中，和分别是矩阵C^(k)的第l个奇异值、第l个左奇异向量和第l个右奇异向量，所对应的矩阵表示形式分别为Σ^(k)，U^(k)和V^(k)，H为取矩阵的共轭转置；L^(k)为C^(k)的秩，可由X^(k)(r,t)中的噪声方差||ζ^(k)||₂通过下式估计获得，本实施例中L⁽¹⁾＝2，L⁽²⁾＝10，L⁽³⁾＝10，L⁽⁴⁾＝2，L⁽⁵⁾＝10，L⁽⁶⁾＝2。

${\mathrm{\λ}}_{\hat{L}+1}^{\left(k\right)}\left(C^{\left(k\right)}\right)\≤\left|\right|{\mathrm{\ζ}}^{\left(k\right)}|{|}_2\≤{\mathrm{\λ}}_{\hat{L}}^{\left(k\right)}\left(C^{\left(k\right)}\right),$

其中，为C^(k)的第个特征值。

(3)基于奇异值分解结果，得到各分部数据的分部滤波处理结果为 ${\hat{X}}^{\left(k\right)}(r,t)=U^{\left(k\right)}{\mathrm{\Σ}}^{\left(k\right)}{V}^{\left(k\right)}H,$ (k＝1,2,...,6)。

第四步，将各分部滤波结果累加求和，得到噪声抑制后的处理结果为：

$\hat{X}(r,t)=\underset{k=1}{\overset{K}{\Σ}}{\hat{X}}^{\left(k\right)}(r,t).$

本实施例待处理冰雷达数据第900道数据幅度序列和相位序列波形如图3所示，其对应包络如图4所示。本实施例噪声抑制结果、第900道数据噪声抑制后信号波形图、第900道数据噪声抑制后灰度梯度图分别如图5、图6和图7所示。与图5(b)、图6(b)和图7(b)小波滤波处理噪声抑制结果相比可见，本发明方法用于极地探冰雷达数据噪声抑制时，可以更有效的反映出不同介质层位的属性差异，通过对不同介质层位的分部式噪声抑制处理，反映出了冰下微弱信号变化，有效保留了冰下微弱层位特征，如图6(a)和图7(a)；在提高图像信噪比的同时，有效保证了图像的空间分辨率。

上述说明已经充分介绍了本发明的具体实施方式。需要指出的是，熟悉该领域的技术人员对本发明的具体实施方式所做的任何改动均不脱离本发明的权利要求书的范围。相应的，本发明的权利要求的范围也并不仅仅局限于前述具体实施方式。

Claims (4)

1.一种针对极地探冰雷达数据的分部式噪声抑制方法，其特征在于，包括以下步骤：

(1)对待处理冰雷达数据X(r,t)进行预滤波处理，得到预处理后的数据

(2)对计算各道冰雷达数据的幅度和相位，基于此对待处理冰雷达数据X(r,t)进行分部划分，得到K个分部冰雷达数据X^(k)(r,t),(k＝1,2,…,K)；

(3)对各分部冰雷达数据X^(k)(r,t),(k＝1,2,…,K)，分别进行最优低秩滤波处理，得到各分部滤波处理结果 ${\hat{X}}^{\left(k\right)}(r,t),(k=1,2,...,K);$

(4)将各分部滤波结果累加求和，得到噪声抑制后的处理结果为：

$\hat{X}(r,t)=\underset{k=1}{\overset{K}{\Σ}}{\hat{X}}^{\left(k\right)}(r,t).$

2.根据权利要求1所述的针对极地探冰雷达数据的分部式噪声抑制方法，其特征在于，所述步骤(1)的对待处理冰雷达数据X(r,t)进行预滤波处理，具体是：对数据X(r,t)采用奇异值分解方法，进行低秩滤波处理。

3.根据权利要求1所述的针对极地探冰雷达数据的分部式噪声抑制方法，其特征在于，所述步骤(2)的计算各道雷达数据的幅度和相位，基于此对待处理冰雷达数据X(r,t)进行分部划分，具体是：对滤波后各道冰雷达数据，计算得到幅度序列和相位序列；绘制幅度序列和相位序列的包络，设定局部极值点判定阈值，得到各包络变化局部极大值；结合幅度序列和相位序列包络变化局部极大值分布，定义不同介质的分部参考位置，对待处理冰雷达数据X(r,t)进行分部划分。

4.根据权利要求1所述的针对极地探冰雷达数据的分部式噪声抑制方法，其特征在于，所述步骤(3)的对各分部冰雷达数据X^(k)(r,t),(k＝1,2,…,K)，分别进行低秩滤波处理，具体是：对X^(k)(r,t),(k＝1,2,…,K)，分别构造其Casorati矩阵如下：

其中，代表在空间位置r_i和采样时间t_j的空间-时间采样值，N和M分别为第k个分部观测数据的空间和时间长度；对矩阵C^(k)，采用奇异值分解方法进行最优低秩估计；假设C^(k)的秩为L^(k)，则其低秩估计结果如下：

$C^{\left(k\right)*}=\underset{l=1}{\overset{L^{\left(k\right)}}{\Σ}}{\mathrm{\σ}}_l^{\left(k\right)}{u}_l^{\left(k\right)}{v}_l^{\left(k\right)H},$

其中，和分别是矩阵C^(k)的第l个奇异值、第l个左奇异向量和第l个右奇异向量，所对应的矩阵表示形式分别为Σ^(k)，U^(k)和V^(k)，H为取矩阵的共轭转置；L^(k)为C^(k)的秩，可由X^(k)(r,t)中的噪声方差||ζ^(k)||₂通过下式估计获得

${\mathrm{\λ}}_{\hat{L}+1}^{\left(k\right)}\left(C^{\left(k\right)}\right)\≤\left|\right|\mathrm{\ζ}\left(k\right)|{|}_2\≤{\mathrm{\λ}}_{\hat{L}}^{\left(k\right)}\left(C^{\left(k\right)}\right),$

其中，为C^(k)的第个特征值；由此，各分部数据的分部滤波处理结果为 ${\hat{X}}^{\left(k\right)}(r,t)=U^{\left(k\right)}{\mathrm{\Σ}}^{\left(k\right)}{V}^{\left(k\right)H},(k=1,2,...,K).$