CN103149562B - 基于极化合成孔径雷达的植被下地形方位角估计方法 - Google Patents

Abstract

本发明属于雷达遥感图像处理技术领域，提出了一种基于极化合成孔径雷达的植被下地形方位角估计方法。本发明利用极化合成孔径雷达的强穿透性，对植被覆盖区域只进行一次观测，在电磁散射理论和数学建模的基础上，将雷达观测数据中植被下地形起伏信息估计出来。由于雷达可以在全天时、全天候的条件下进行地表观测，估计的植被下地形信息可以为滑坡检测、灾害应急响应、国土资源调绘等应用提供充分保障。由于只使用一张POLSAR影像，本发明可大大降低购买雷达数据的经济费用，适用范围大大增加。

Description

基于极化合成孔径雷达的植被下地形方位角估计方法

技术领域

本发明属于雷达遥感图像处理技术领域，特别设计一种利用极化合成孔径雷达（Polarimetric Synthetic Aperture Radar，POLSAR）对植被下地形进行估计的新方法。

背景技术

数字地形模型用于描述地貌形态的空间分布，其作为国家基础测绘数据，对于国家的重大基础设施建设、资源开发与利用、生态建设与环境保护以及军事中对目标的定位导航都有着十分重大的意义。1983年美国国家航天宇航局NASA统计表明：全球约33%的陆地范围为森林所覆盖，对于这些地区由于有大量的植被存在，传统测量方法获取的地形无法保证高的精度。对于我国一些大面积植被覆盖区域和西部地区南疆沙漠、青藏高原和横断山脉地区的大片国土，由于受自然地理条件和技术水平的限制，仍无1:50000和1:100000比例尺的基础地理图件，使得对这些地区的资源勘察和开发利用工作难以展开。这些地区气候条件恶劣，地形复杂，单纯依靠传统光学影像不能够完成该地区地形测量和地物覆盖分类及各项资源勘察任务，合成孔径雷达（SAR）由于全天时、全天候的特点，可以弥补光学遥感的不足，与光学遥感在西部测绘中互为补充完成各项测图任务。虽然国家在“十一五”期间启动西部侧图计划，对部分困难区域考虑合成孔径雷达干涉测量（INSAR）和光学技术联合应用来获取的1:50000和1:100000比例尺的三维地形图，已经完成了困难地区大部分区域的测图；但对于有森林植被覆盖区域，依靠INSAR技术和光学三维测量技术都不能够获取植被下DEM信息，必须采用新的手段来获取植被覆盖区域的地形。目前，利用极化干涉合成孔径雷达（Polarimetric Interferometric Synthetic Aperture Radar，POLINSAR）技术进行植被下地形提取是国内外研究的一个热点，但利用POLINSAR技术提取地形需要对同一地区进行多次重复观测，而重复观测受到经济成本制约，不能大面积生产。

发明内容

一般而言，每一次观测的POLSAR影像包含了反映地物散射机理的4个通道数据。针对现有技术的不足，本发明目的在于利用单次POLSAR雷达影像中观测的4个通道，通过数学和电磁散射模型，将植被与植被下地表散射矩阵分离，然后通过植被下散射矩阵将地形信息分离出来。

本发明的技术方案为一种基于极化合成孔径雷达的植被下地形方位角估计方法，包括以下步骤：

步骤1，雷达影像噪声去除，包括以每个待恢复点为中心开一个大小为n×n的窗口，用窗口中的所有像素进行多视处理估计极化相干矩阵；然后，在n×n的窗口内，对极化相干矩阵每一个元素进行快速傅立叶变换，在频率域构建低通滤波器，将快速傅立叶变换后的极化相干矩阵输入低通滤波器，进行逆傅立叶变化获取低通滤波以后的低噪声极化相干矩阵，作为观测值矩阵T_observer；

步骤2，针对植被下地形有关的参数θ_g，构建植被层与植被下地表模型如下，

T_observer=T_v(f_v,θ_v，,σ)+T_g(f_g,α,β，θ_g)

其中，T_v(f_v,θ_v,σ)与T_g(f_g,α,β，θ_g)分别为植被层矩阵与植被下地表矩阵，

T_v(f_v,θ_v,σ)由以下植被层模型提供，

$T_v(f_v,{\mathrm{\θ}}_v,\mathrm{\σ})=\frac{1}{16}{f}_v\left[\begin{array}{ccc}8& -4q\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)\\ -4p\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4+4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)\\ 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)& 4-4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)\end{array}\right]$

T_g(f_g,α,β，θ_g)由以下植被下地表模型提供，

$T_g(f_g,\mathrm{\α},\mathrm{\β},{\mathrm{\θ}}_g)=R\left({\mathrm{\θ}}_g\right)\·\frac{f_g}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}1& \mathrm{\β}& 0\\ {\mathrm{\β}}^{*}& {\left|\mathrm{\α}\right|}^2& 0\\ 0& 0& 0\end{array}\right]\·R{\left({\mathrm{\θ}}_g\right)}^T$

$R\left({\mathrm{\θ}}_g\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos 2\mathrm{\θ}}_g& \sin {2\mathrm{\θ}}_g\\ 0& -\sin 2{\mathrm{\θ}}_g& \cos {2\mathrm{\θ}}_g\end{array}\right]$

其中，f_v、θ_v、σ为植被层矩阵T_v(f_v,θ_v,σ)包含的三个未知参数，分别为植被层散射强度、植被层方位偏移角、植被层随机度，p(σ)与q(σ)为以σ为输入量的多项式函数；f_g、θ_g为地表散射强度、植被下地形方位偏移角，α与β为地表散射系数，β^*为β的共轭系数；

步骤3，根据步骤1所得观测值矩阵T_observer，求解步骤2中植被层与植被下地表模型的未知参数θ_g。

而且，步骤3中，采用不同的初始值，利用trust region reflective非线性方程求解法对植被层模型三个未知参数进行备选解求解，使得求解后观测矩阵与植被层矩阵之差的最小特征值为0；然后将地体强度比小于1的备选解剔除；最后，将地面散射机理最接近二面散射或单次散射的备选解视为真实解，获取真实植被下地形起伏信息。

本发明的创新之处在于，POLSAR影像通常由植被层与起伏地形下的地表层两者组成，同时观测数据又包含了大量噪声。在频率域去除噪声的基础上，本发明利用物理建模的方式，将POLSAR影像自适应的分解为植被层与起伏地表，然后在去除植被影响的条件下获取植被下地形的起伏信息。这种基于单次POLSAR雷达观测的新植被下地形估计方法，由于只使用单次观测POLSAR影像，避免了对同一地区的重复观测，可大大降低购买雷达数据的经济费用，适用范围大大增加，具有广阔的应用前景和经济价值。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的植被下地形信息估计流程图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

为了提取植被下地形，本发明实施例首先需要对原始的POLSAR雷达影像采用低通滤波，将数据中不稳定的随机噪声去除。每一次观测的POLSAR影像包含了4个通道数据，由于只有单次观测POLSAR雷达影像，建立适当的植被、植被下地形散射模型将两者分离是必要的。本实施例采用国际知名POLSAR植被层模型，同时将植被下地形起伏信息融入到地形模型中，并采用适当的求解方法，将植被下地形求解出来。

本发明技术方案可采用计算机技术实现自动运行流程。如图1所示，本发明实施例提供的流程包括以下几个步骤：

步骤1，POLSAR雷达影像噪声去除。

实施例以每个待恢复点为中心开一个大小为n×n的窗口；用窗口中的所有像素进行多视处理估计极化相干矩阵。具体实施时，本领域技术人员可以根据具体影像预设窗口尺寸n值的大小。然后，在n×n的窗口内，对极化相干矩阵每一个元素进行快速傅立叶变换（FFT）。在频率域构建低通滤波器，将FFT后的极化相干矩阵输入低通滤波器，然后进行逆傅立叶变化（IFFT），获取低通滤波以后的低噪声极化相干矩阵。

为便于实施参考起见，提供实施例的POLSAR影像噪声去除具体实现方式如下：

首先对原始极化散射矩阵S2进行Pauli基（保里基）变换，在n×n的窗口内构建极化相干矩阵T_obsever：

$S2=\left[\begin{array}{cc}{\mathrm{HH}}_i& {\mathrm{HV}}_i\\ {\mathrm{VH}}_i& {\mathrm{VV}}_i\end{array}\right]\→T_{\mathrm{oberver}}=\underset{i=1}{\overset{n\×n}{\mathrm{\Σ}}}{k}_i\·k_i^{*T}$

其中，HH_i、HV_i、VH_i、VV_i为雷达HH通道观测到的影像、雷达HV通道观测到的影像、雷达VH通道观测到的影像、雷达VV通道观测到的影像的相应n×n窗口内第i个像素，i的取值为1,2…n×n。复数Pauli矢量k_i为：

$k_i=\frac{1}{\sqrt{2}}{[{\mathrm{HH}}_i+{\mathrm{VV}}_i,{\mathrm{HH}}_i-{\mathrm{VV}}_i,{\mathrm{HV}}_i+{\mathrm{VH}}_i]}^T$

其中，T表示矩阵转置，*表示共轭操作。

然后，对估计的矩阵T_observer中每一个元素进行FFT变换，在频率域进行Hamming（汉明）低通滤波，再将低通滤波后结果通过IFFT变换回到空间域，从而完成影像去噪声操作，所得低通滤波以后的低噪声极化相干矩阵T_observer作为观测值矩阵。

步骤2，构建植被层与植被下地表模型。

实施例针对植被层采用国际知名的POLSAR植被自适应模型，其将植被层在雷达上的散射特性归纳为三个待求参数：植被散射强度、植被层方位偏移角、植被层随机度。对于植被下地形，实施例采用考虑地形起伏的改进散射模型，在雷达算法的基础上，实施例利用数学建模和电磁散射理论进行了相关的改进，所得模型包括以下几个待求参数：地表散射强度、植被下地形方位偏移角、地表散射系数等。

为便于实施参考起见，提供实施例构建植被层与植被下地表模型作为POLSAR影像自适应模型的具体实现方式如下：

T_observer=T_v(f_v,θ_v,σ)+T_g(f_g,α,β，θ_g)

T_v(f_v,θ_v,σ)与T_g(f_g,α,β，θ_g)分别为植被层矩阵与植被下地表矩阵，其中：

植被层模型如下：

$T_v(f_v,{\mathrm{\θ}}_v,\mathrm{\σ})=\frac{1}{16}{f}_v\left[\begin{array}{ccc}8& -4q\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)\\ -4p\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4+4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)\\ 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)& 4-4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)\end{array}\right]$

改进的植被下地表模型如下：

$T_g(f_g,\mathrm{\α},\mathrm{\β},{\mathrm{\θ}}_g)=R\left({\mathrm{\θ}}_g\right)\·\frac{f_g}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}1& \mathrm{\β}& 0\\ {\mathrm{\β}}^{*}& {\left|\mathrm{\α}\right|}^2& 0\\ 0& 0& 0\end{array}\right]\·R{\left({\mathrm{\θ}}_g\right)}^T$

矩阵 $R\left({\mathrm{\θ}}_g\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos 2\mathrm{\θ}}_g& \sin {2\mathrm{\θ}}_g\\ 0& -\sin 2{\mathrm{\θ}}_g& \cos {2\mathrm{\θ}}_g\end{array}\right]$

其中，f_v、θ_v、σ为植被层矩阵T_v(f_v,θ_v,σ)包含的三个未知参数，分别为植被层散射强度、植被层方位偏移角、植被层随机度，p(σ)与q(σ)为以σ为输入量的多项式函数；f_g、θ_g为地表散射强度、植被下地形方位偏移角，α与β为地表散射系数，β^*为β的共轭系数。构建的模型中θ_g是与植被下地形起伏息息相关的待求解变量，而模型中植被信息与地表信息相互混叠，需要采用特别的方法将地形信息从观测值中估计出来。

步骤3，将步骤1所得低通滤波以后的低噪声极化相干矩阵T_observer视为已知观测量，求解步骤2中的未知参数。步骤2中的模型通常是一个极其复杂的非线性方程，并且未知数个数多于观测值，所以本发明并不直接求解所有未知数，而是仅仅针对植被下地形有关的参数进行求解。本发明中，采用随机生成的不同的初始值，利用trust region reflective非线性方程求解法对植被层模型三个未知参数进行备选解求解，使得求解后观测值矩阵T_observer与植被层矩阵T_v(f_v,θ_v,σ)之差的最小特征值为0。然后在此基础上将地表-植被散射强度（即地体强度比比小于1的备选解剔除。最后，将地面散射机理最接近二面散射或单次散射的备选解视为真实解，从而获取真实植被下地形起伏信息。

如图2所示，为便于实施参考起见，提供实施例的植被下地形信息估计具体实现方式如下：

首先，以植被层模型三个未知参数f_v、θ_v、σ为待求解量，构建目标函数：

$\underset{f_v,\mathrm{\σ},{\mathrm{\θ}}_v}{\arg }\min \mathrm{imum}:\mathrm{eigenvalue}\{T_{\mathrm{observer}}-T_v(f_v,{\mathrm{\θ}}_v,\mathrm{\σ})\}=0$

其中T_obsener为观测值矩阵，eigenvalue{.}表示矩阵特征值。目标函数尝试求解三个未知数，使得观测值与植被模型差的最小特征值为0，具体求解方法可采用trust region reflective(信赖域反射算法)非线性方程求解法。满足目标函数的解通常有多个，设满足条件的M个备选解所构成解集X为：

$X=\{{\mathrm{\σ}}^j,{\mathrm{\θ}}_v^j,f_v^j|j=1,...,M\}$

其中σ^j、为满足条件的第j个解，j的取值为1,2…M。然后，对每一个解都重新构建植被层矩阵求地体强度比。定义地体强度比

$\mathrm{GrV}({\mathrm{\σ}}^j,{\mathrm{\θ}}_v^j,f_v^j)=\frac{\mathrm{Tr}\{T_{\mathrm{observer}}-T_v({\mathrm{\σ}}^j,{\mathrm{\θ}}_v^j,f_v^j)\}}{\mathrm{Tr}\left\{T_v({\mathrm{\σ}}^j,{\mathrm{\θ}}_v^j,f_v^j)\right\}}$

Tr{·}是矩阵求迹操作。在植被区由于植被层的散射强度往往大于地表层强度，所以只保留小于等于1的解。对保留的解构建并与观测矩阵做差后进行Cloude-Pottier正交分解：

$T_{\mathrm{observer}}-T_v({\mathrm{\σ}}^j,{\mathrm{\θ}}_v^j,f_v^j)=\left(u_1{u}_2{u}_3\right)\·\left(\begin{array}{ccc}0& & \\ & {\mathrm{\λ}}_2& \\ & & {\mathrm{\λ}}_3\end{array}\right)\·{\left(u_1{u}_2{u}_3\right)}^H$

其中u₁、u₂、u₃分别为矩阵三个特征值对应的正交特征矢量，λ₂、λ₃为第二、第三特征值。

H表示共轭矩阵。

获取极化散射角α_j：

α_j=acosu₂(1)

u₂(1)为第二特征矢量的第一个元素。由于α_j表征了当前散射机理类型，纯净的地表散射机理通常为0度或90度。可以通过寻找α_j最接近0度或90度的解作为最优解α：

α=argmax{|α_j-45°|}

获取最优解α对应的进行极化方位角旋转，从而得到真实植被下地表起伏信息θ_g。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (1)

1.一种基于极化合成孔径雷达的植被下地形方位角估计方法，其特征在于，包括以下步骤：

步骤1，雷达影像噪声去除，包括以每个待恢复点为中心开一个大小为n×n的窗口，用窗口中的所有像素进行多视处理估计极化相干矩阵；然后，在n×n的窗口内，对极化相干矩阵每一个元素进行快速傅立叶变换，在频率域构建低通滤波器，将快速傅立叶变换后的极化相干矩阵输入低通滤波器，进行逆傅立叶变化获取低通滤波以后的低噪声极化相干矩阵，作为观测值矩阵T_observer；

步骤2，针对植被下地形有关的参数θ_g，构建植被层与植被下地表模型如下，

T_observer＝T_v(f_v,θ_v,σ)+T_g(f_g,α,β,θ_g)

其中，T_v(f_v,θ_v,σ)与T_g(f_g,α,β,θ_g)分别为植被层矩阵与植被下地表矩阵，

T_v(f_v,θ_v,σ)由以下植被层模型提供，

$T_v(f_v,{\mathrm{\θ}}_v,\mathrm{\σ})=\frac{1}{16}{f}_v\left[\begin{array}{ccc}8& -4p\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)\\ -4p\left(\mathrm{\σ}\right)\cos \left({2\mathrm{\θ}}_v\right)& 4+4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)\\ 4p\left(\mathrm{\σ}\right)\sin \left({2\mathrm{\θ}}_v\right)& -4q\left(\mathrm{\σ}\right)\sin \left({4\mathrm{\θ}}_v\right)& 4-4q\left(\mathrm{\σ}\right)\cos \left({4\mathrm{\θ}}_v\right)\end{array}\right]$

T_g(f_g,α,β,θ_g)由以下植被下地表模型提供，

$T_g(f_g,\mathrm{\α},\mathrm{\β},{\mathrm{\θ}}_g)=R\left({\mathrm{\θ}}_g\right)\·\frac{f_g}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}1& \mathrm{\β}& 0\\ {\mathrm{\β}}^{*}& {\left|\mathrm{\α}\right|}^2& 0\\ 0& 0& 0\end{array}\right]\·R{\left({\mathrm{\θ}}_g\right)}^{}$

$R\left({\mathrm{\θ}}_g\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos 2\mathrm{\θ}}_g& \sin {2\mathrm{\θ}}_g\\ 0& {-\sin 2\mathrm{\θ}}_g& {\cos 2\mathrm{\θ}}_g\end{array}\right]$

其中，f_v、θ_v、σ为植被层矩阵T_v(f_v,θ_v,σ)包含的三个未知参数，分别为植被层散射强度、植被层方位偏移角、植被层随机度，p(σ)与q(σ)为以σ为输入量的多项式函数；f_g、θ_g为地表散射强度、植被下地形方位偏移角，α与β为地表散射系数，β^*为β的共轭系数；步骤3，根据步骤1所得观测值矩阵T_observer，求解步骤2中植被层与植被下地表模型的未知参数θ_g，实现方式如下，

采用不同的初始值，利用trust region reflective非线性方程求解法对植被层模型三个未知参数进行备选解求解，使得求解后观测矩阵与植被层矩阵之差的最小特征值为0；然后将地体强度比小于1的备选解剔除；最后，将地面散射机理最接近二面散射或单次散射的备选解视为真实解，获取真实植被下地形起伏信息。