CN112363161B - 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 - Google Patents
基于散射机制分解的植被垂直结构及林下地形反演方法及装置 Download PDFInfo
- Publication number
- CN112363161B CN112363161B CN202011164833.3A CN202011164833A CN112363161B CN 112363161 B CN112363161 B CN 112363161B CN 202011164833 A CN202011164833 A CN 202011164833A CN 112363161 B CN112363161 B CN 112363161B
- Authority
- CN
- China
- Prior art keywords
- scattering
- canopy
- covariance matrix
- polarization
- forest
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于散射机制分解的植被垂直结构及林下地形反演方法及装置,其方法为:获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地处理,得到极化层析SAR数据;对极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,得到对应的稀疏贝叶斯学习模型;采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;根据地表散射和冠层散射的层析谱,分别提取森林反演区域的林下地形和植被高度。本发明可以利用少量、非均匀分布的SAR数据高精度地反演植被高和林下地形。
Description
技术领域
本发明涉及雷达遥感对地观测技术,具体涉及一种在少量非均匀采样数据条件下基于散射机制分解的超分辨率极化SAR层析植被垂直结构及林下地形反演方法及装置。
背景技术
生态文明建设是我国发展的重要布局之一,森林资源的调查与保护在其中扮演着重要的角色。森林中包含着种类繁多的植被,是地球上最宝贵的自然资源之一。森林在维持全球的碳氧平衡、维护生态***的多样性、调节大气湿度等方面都有着巨大的作用。森林中植被的垂直结构是量化分析森林的重要参数,对全球生物量的估计、林下地形测量、森林动态监测等有着重要的意义。
长波段合成孔径雷达技术(Synthetic Aperture Radar,SAR)具有较强的穿透力,在森林的透视测绘中有着重要的应用。合成孔径雷达层析技术(Synthetic ApertureRadar Tomography,TomoSAR)通过传感器多次近似平行的飞行,在高度向和方位向上分别形成合成孔径,具备三维成像能力(图1)。在雷达信号中,森林植被与林下地形有着不同的散射机制,而SAR影像不同极化通道对不同散射机制的敏感性不同。因此极化层析SAR技术(Polarimetric Synthetic Aperture Radar Tomography,Pol-TomoSAR)正成为森林植被参数反演的重要技术手段。
在植被的SAR层析领域,周期图法、自适应波束形成法、多重信号分类法等已广泛用于植被的三维成像。但这类常规谱分析方法的分辨率和重构精度往往受到采样数量和采样分布的限制。观测数量少、基线分布不均匀往往会使得层析结果出现严重的旁瓣效应、分辨率低、受噪声影响大。而增加观测数据则意味着大量的观测成本,同时带来时间失相干的影响。压缩感知方法可以很好地克服数据采样的限制,但传统基于l1范数凸优化的压缩感知方法(Compressive Sensing,CS)需要设置用户参数,很难选择一个普遍适用于所有模型和场景的参数。此外,凸优化运算量巨大,限制了该方法的应用。
发明内容
基于上述现有技术存在的技术问题,本发明提供一种基于散射机制分解的植被垂直结构及林下地形反演方法及装置,将森林反演区域冠层和地表两种不同散射机制的信号进行稀疏表达,进而将其纳入协方差矢量的稀疏贝叶斯学习模型中以求解对应的层析谱,从而实现少量的观测数据重构高精度的植被垂直结构和林下地形。本发明具有更接近于l0范数的全局最优解,无需设置用户参数,具有更广泛的实用价值。
为实现上述技术目的,本发明采用如下技术方案:
一种基于散射机制分解的植被垂直结构及林下地形反演方法,包括以下步骤:
步骤1,获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地处理,得到极化层析SAR数据;
步骤2,对步骤1得到的极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;
步骤3,针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型;
步骤4,采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;
步骤5,根据地表散射的层析谱提取森林反演区域的林下地形,根据冠层散射的层析谱提取森林反演区域的植被高度。
在更优的技术方案中,步骤2中多基线多极化协方差矩阵的计算方法为:
式中,y为步骤1得到的极化层析SAR数据,gHV=gVH,gHH、gHV、gVH、gVV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测矩阵,N为基线的数量;(*)H为共轭转置操作,E[·]为求期望运算;W为计算得到的多基线多极化协方差矩阵;Cg和Rg分别为地表散射的极化协方差矩阵和干涉协方差矩阵,Cv和Rv分别为冠层散射的极化协方差矩阵和干涉协方差矩阵,为克罗内克积;
对多基线多极化协方差矩阵W进行SVD分解,使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg。
在更优的技术方案中,使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg的方法为:
首先,取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R1和R2,使用R1和R2表示冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg为:
式中,a、b为一对实数;
然后,基于冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg的半正定约束,即Rv和Rg的特征值非负,求解实数a、b的取值区间;
再后,根据以下散射机制最大化分离的准则,求解实数a、b:
其中,trace()为求矩阵的迹,||·||F为矩阵的Frobenius范数。最后,根据求解得到的实数a、b和矩阵R1和R2,即可得到冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg。
在更优的技术方案中,由步骤3得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为:
式中,Rg、Rv分别为地表散射和冠层散射的干涉协方差矩阵,vec(·)为矢量化操作,yg、yv分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示;
B为中间参数,且B=A⊙conj(A),⊙为Khatri-Rao积,A为中间参数,A=exp(j2πξnsd),sd为散射体的空间高度,ξn为空间频率;Ψg表示地表散射的小波稀疏基,采用单位正交基,Ψv表示冠层散射的小波稀疏基;
pg、pv分别为地表散射和冠层散射的后向散射功率矢量;
σg、σv分别为地表散射和冠层散射对应的噪声功率,tg=vec(σgI)、tv=vec(σvI)分别为对应噪声功率矩阵的矢量化操作,I为单位矩阵;δg和δv分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差。
在更优的技术方案中,步骤4中采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数并获得冠层散射和地表散射的层析谱,具体为:
步骤a,定义先验分布:
为了便于阐述,将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为:Y=Gx+δ;其中,G=[BΨH,I]表示在稀疏基下新的测量矩阵,x=[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量,Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示;
根据稀疏操作后的散射信号是稀疏这一先验信息,定义其服从高斯分布,即x~N(0,Λ);其中Λ=diag(Γ),为稀疏贝叶斯学习模型的超参数矢量,表征了x的空间分布;前D个超参数τ1,...τD控制高度向上D个散射体后向散射功率的稀疏性,后N2个超参数为***噪声;则Y关于Γ的概率密度函数为:
p(Y|Γ)=∫p(Y|x)p(x|Γ)dx;
步骤b,参数初始值定义:
设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件;
步骤c,计算后验概率分布:
基于贝叶斯准则和步骤a定义的先验分布,计算x的均值矢量和协方差矩阵:
ΣY=V+GΛ(i)GH
式中,上标i表示当前迭代次数,μ为x的均值矢量,Σx和ΣY分别为x和Y的协方差矩阵,V为样本协方差矩阵估计误差δ的二阶统计。
步骤d,基于步骤c得到的均值矢量和协方差矩阵,采用EM算法学习更新超参数向量:
步骤e,重复步骤c、d,直至满足迭代终止条件,使用超参数向量Γ中的前D个超参数Γ(1…D)计算冠层散射和地表散射的层析谱:
在更优的技术方案中,迭代终止条件为||Γi+1-Γi||2/||Γi||2≤tol,和/或达到最大迭代次数。
在更优的技术方案中,根据地表散射的层析谱提取森林反演区域的林下地形的方法为:提取地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。
在更优的技术方案中,根据冠层散射的层析谱提取森林反演区域的植被高度的方法为:通过功率衰减法确定冠层顶部的位置,即:以冠层相位中心处为起始点,以预设的步长向上衰减,当功率衰减处提取的高度值和Lidar数据提取的冠层高度(DSM)的均方根误差达到最小时,将此时的高度值作为冠层高度的估计值;通过冠层高度减去地面高度即为植被高度。
本发明还提供一种基于散射机制分解的植被垂直结构及林下地形反演装置,包括:
数据预处理模块,用于:获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地处理,得到极化层析SAR数据;
协方差计算分解模块,用于:对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;
模型构建模块,用于:针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型;
层析谱计算模块,用于:采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;
植被垂直结构及林下地形反演模块,用于:根据地表散射的层析谱提取森林反演区域的林下地形,根据冠层散射的层析谱提取森林反演区域的植被高度。
有益效果
本发明使用多基线全极化的SAR数据,对冠层和地表两种不同散射机制进行了分离与层析。一方面,可以获得更加准确的相位中心估计值,提高林下地形与植被高度的反演精度;另一方面,使用不同稀疏基对不同散射机制进行稀疏表达,使得本发明方法能在森林植被区域合理地运用。
而且,本发明方法是一种近似于l0范数的稀疏重构方法,相较于传统基于l1范数凸优化的压缩感知方法,可以获得更稀疏的全局最优解,能够以较少的测量数据满足RIP条件(即感知矩阵G列不相关);无需艰难地设置噪声上限参数,可以自适应噪声的影响,具备更高的植被高和林下地形重构精度;同时,本发明的方法具备更高的运算效率,适用于大范围、大尺度的林下地形测绘以及植被高度反演,具有更高的生产实用价值。
另外,本发明顾及了层析SAR数据样本协方差矩阵的估计误差,该方法更加稳健。
附图说明
图1是层析SAR森林三维成像示意图;
图2是本发明实施例所述方法的流程图;
图3是本发明实施例所述方法各步骤子流程图。
具体实施方式
为了更加清楚地说明本发明的目的、实施流程和技术优势,本实施例选取Tropi-SAR 2009项目ONERA SETHI飞行***获得的P波段全极化的SAR数据,对本发明所提方法进行进一步的详细说明。应当理解,此处所描述的具体实施仅用于解释本发明,并不用于限定本发明。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例,都属于本发明保护的范围。
本发明基于散射机制分解的植被垂直结构及林下地形反演方法,为超分辨率的稀疏重构方法。为了说明方法的优越性能,具体实施中分析了不同层析SAR方法在不同数目的全极化SAR数据下的性能,然后仅采用3条基线的全极化SAR数据反演林下地形和植被高度。
本实施例的具体实施步骤参照图2、图3所示,包括:
步骤1,获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地处理,得到极化层析SAR数据。具体为:
(1)选取TropiSAR2009项目全极化SAR影像作为实验数据。
(2)配准:对于指定的极化通道(例如:HH极化通道),在获取的N景SAR影像中,选择其中一景影像为参考影像,其余N-1景影像应用最大化相干系数配准到参考影像的坐标框架中。然后将剩余极化通道(HV、VH、VV)的SAR影像采用与HH极化通道相同的多项式进行配准。
(3)去平地:即去除由平地效应引起的干涉相位。
(4)数据离散化:经上述处理后,对于选定极化方式的极化SAR层析数据,距离-方位向平面像素(r,x)聚焦后的复数值为:
式中bn为垂直基线,Δs为高度向的采样范围,γ(r,x,s)为高度向上的反射率函数,ξn=2πbn/λr为空间频率。
HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值分别为gHH、gHV、gVH、gVV;N为基线的数量;在单站后向散射体制下gHV=gVH,因此可将步骤1最终得到的极化层析SAR数据表示为对于多基线的层析SAR数据集,对高度向上连续的信号进行D次离散的采样后,可将(1)表达为以下矩阵形式:
g=Aγ+e (2)
其中g是N维观测值。A=exp(j2πξnsd),sd为散射体的空间高度,e是N维的噪声向量。
步骤2,对极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵。
对于森林植被的SAR层析而言,往往利用不同极化通道对不同散射机制的敏感性不同来提取林下地形和植被高度,例如在HH极化通道的层析谱中提取林下地形,在HV极化通道的层析谱中提取冠层高度。但由于森林场景较为复杂,各极化后向散射能量并不完全具备如此规律的分布,前述提取不同散射机制的策略可能会造成一定的偏差。此外,本发明基于散射机制分解的极化SAR层析技术要求所测量的信号是稀疏的(信号矢量中只有少量元素不为零)。在森林区域,植被信号的散射是一个连续的过程,而地面信号本身具备稀疏性,因此本发明对冠层和地表两种不同散射机制的信号采用不同的稀疏基进行稀疏表达。
为了更加准确地提取出不同散射机制相位中心的位置,以及对不同的散射机制进行稀疏表达,本发明对步骤1得到的极化层析SAR数据计算其协方差矩阵,即得到多基线多极化协方差矩阵,进而使用代数合成的方法对两种不同的散射机制(地表散射、冠层散射)进行分离:
(*)H为共轭转置操作,E[·]为求期望运算;W为计算得到的多基线多极化协方差矩阵;Cg和Rg分别为地表散射的极化协方差矩阵和干涉协方差矩阵,Cv和Rv分别为冠层散射的极化协方差矩阵和干涉协方差矩阵,为克罗内克积。
得到多基线多极化协方差矩阵W后,对其进行SVD分解,使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg:
首先,取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R1和R2,使用矩阵R1和R2表示冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg为:
式中,a、b为一对实数;
然后,基于冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg的半正定约束,即Rv和Rg的特征值非负,求解实数a、b的取值区间;
再后,根据以下散射机制最大化分离的准则,求解实数a、b:
其中,trace()为求矩阵的迹,||·||F为矩阵的Frobenius范数。最后,根据求解得到的实数a、b和矩阵R1和R2,即可得到冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg。
步骤3,针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型。
在步骤2分解得到冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg后,分别进行矢量化和稀疏表达,得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为:
式中,Rg、Rv分别为地表散射和冠层散射的干涉协方差矩阵,vec(·)为矢量化操作,yg、yv分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示;
B为中间参数,且B=A⊙conj(A),⊙为Khatri-Rao积,A为中间参数,A=exp(j2πξnsd),sd为散射体的空间高度,ξn为空间频率;Ψg表示地表散射的小波稀疏基,采用单位正交基,Ψv表示冠层散射的小波稀疏基;
pg、pv分别为地表散射和冠层散射的后向散射功率矢量;
σg、σv分别为地表散射和冠层散射对应的噪声功率,tg=vec(σgI)、tv=vec(σvI)分别为对应噪声功率矩阵的矢量化操作,I为单位矩阵;δg和δv分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差,其二阶统计量分别为Vg,Vv。
步骤4,采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱。具体为:
步骤a,定义先验分布:
为了便于阐述,将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为:Y=Gx+δ;其中,G=[BΨH,I]表示在稀疏基下新的测量矩阵,x=[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量,Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示;
根据稀疏操作后的散射信号(即x)是稀疏这一先验信息,定义其服从高斯分布,即x~N(0,Λ);其中Λ=diag(Γ),为稀疏贝叶斯学习模型的超参数矢量,表征了x的空间分布;前D个超参数τ1,...τD控制高度向上D个散射体后向散射功率的稀疏性,后N2个超参数为***噪声;则Y关于Γ的概率密度函数为:
p(Y|Γ)=∫p(Y|x)p(x|Γ)dx;
步骤b,参数初始值定义:
设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件,如:(Γ)0是一个D+N2维的随机非负矢量,Λ(0)=diag(Γ(0));迭代终止误差||Γi+1-Γi||2/||Γi||2≤tol=10-4,最大迭代次数为10次;
步骤c,计算后验概率分布:
基于贝叶斯准则和步骤a定义的先验分布,计算x的均值矢量和相关协方差矩阵:
ΣY=V+GΛ(i)GH
式中,上标i表示当前迭代次数,μ为x的均值矢量,Σx和ΣY分别为x和Y的协方差矩阵,V为样本协方差矩阵估计误差δ的二阶统计;
步骤d,基于步骤c得到的均值矢量和协方差矩阵,采用EM算法(即期望最大化算法)学习更新超参数向量:
超参数向量的更新同对以下目标函数进行最小化进行:
步骤e,重复步骤c、d,直至满足任一个迭代终止条件,使用超参数向量Γ中的前D个超参数Γ(1…D)计算冠层散射和地表散射的层析谱(即空间功率谱):
由于已对不同散射机制的信号进行了稀疏化操作,故用于表达后向散射功率的小波系数在稀疏域中大部分位置的值都是0值,故本发明中假设稀疏化后的森林观测信号x服从高斯分布符合先验信息,具有合理性。同时,高斯先验分布十分有利于稀疏贝叶斯学习的开展。
步骤5,根据地表散射的层析谱计算森林反演区域的林下地形,根据冠层散射的层析谱计算森林反演区域的植被高度。具体为:
林下地形的计算方法为:将地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。
植被高度的计算方法为:通过功率衰减法确定冠层顶部的位置,即:以冠层相位中心处为起始点,以预设的步长-0.1db向上衰减,当功率衰减处提取的高度值和Lidar数据DSM的均方根误差达到最小时,将此时的高度值作为冠层高度的估计值;通过冠层高度减去地面高度即为植被高度。
多基线多极化的SAR数据中的每个像元,均按上述方法遍历,即可获得整个森林反演区域的林下地形和植被高度。
通过对实例数据反演的不同散射层析谱、林下地形和植被高度进行分析,本发明方法能够在少量观测数据下获得清晰的层析谱,几乎不受噪声影像,分辨率高。而传统谱分析方法Beamforming和Capon即使使用了全极化的数据,仍然不能完全地区分不同散射机制的散射体。随着基线数量的下降,两种方法的性能受到了很大的影响,以致不能区分不同高度的散射体。而本发明方法受影像数量的影响较小,即使在仅使用3条基线SAR数据的情况下,仍然保持着较高的重构性能。此外,本发明方法仅在仅使用3条基线全极化SAR数据的情况下,也获得了可靠的林下地形和植被高反演结果。
本发明还提供一种基于散射机制分解的植被垂直结构及林下地形反演装置,与上述反演方法对应,包括:
数据预处理模块,用于:获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地处理,得到极化层析SAR数据;
协方差计算分解模块,用于:对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;
模型构建模块,用于:针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型;
层析谱计算模块,用于:采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;
植被垂直结构及林下地形反演模块,用于:根据地表散射的层析谱提取森林反演区域的林下地形,根据冠层散射的层析谱提取森林反演区域的植被高度。
以上实施例为本申请的优选实施例,本领域的普通技术人员还可以在此基础上进行各种变换或改进,在不脱离本申请总的构思的前提下,这些变换或改进都应当属于本申请要求保护的范围之内。
Claims (7)
1.一种基于散射机制分解的植被垂直结构及林下地形反演方法,其特征在于,包括以下步骤:
步骤1,获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地,得到极化层析SAR数据;
步骤2,对步骤1得到的极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;
其中,多基线多极化协方差矩阵的计算方法与分解方法为:
式中,y为步骤1得到的极化层析SAR数据,gHV=gVH,gHH、gHV、gVH、gVV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值;(*)H为共轭转置操作,E[·]为求期望运算;W为计算得到的多基线多极化协方差矩阵;Cg和Rg分别为地表散射的极化协方差矩阵和干涉协方差矩阵,Cv和Rv分别为冠层散射的极化协方差矩阵和干涉协方差矩阵,为克罗内克积;
对多基线多极化协方差矩阵W进行SVD分解,使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg:
首先,取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R1和R2,使用R1和R2表示冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg为:
式中,a、b为一对实数;
然后,基于冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg的半正定约束,即Rv和Rg的特征值非负,求解实数a、b的取值区间;
再后,根据以下散射机制最大化分离的准则,求解实数a、b:
其中,trace()为求矩阵的迹,||·||F为矩阵的Frobenius范数;
最后,根据求解得到的实数a、b和矩阵R1和R2,即可得到冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg;
步骤3,针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型;
步骤4,采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;
步骤5,根据地表散射的层析谱提取森林反演区域的林下地形,根据冠层散射的层析谱提取森林反演区域的植被高度。
2.根据权利要求1所述的方法,其特征在于,由步骤3得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为:
式中,Rg、Rv分别为地表散射和冠层散射的干涉协方差矩阵,vec(·)为矢量化操作,yg、yv分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示;
B为中间参数,且B=A⊙conj(A),⊙为Khatri-Rao积,A为中间参数,A=exp(j2πξnsd),sd为散射体的空间高度,ξn为空间频率;Ψg表示地表散射的小波稀疏基,采用单位正交基,Ψv表示冠层散射的小波稀疏基;
pg、pv分别为地表散射和冠层散射的后向散射功率矢量;
σg、σv分别为地表散射和冠层散射对应的噪声功率,tg=vec(σgI)、tv=vec(σvI)分别为对应噪声功率矩阵的矢量化操作,I为单位矩阵;δg和δv分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差。
3.根据权利要求2所述的方法,其特征在于,步骤4中采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数并获得冠层散射和地表散射的层析谱,具体为:
步骤a,定义先验分布:
为了便于阐述,将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为:Y=Gx+δ;其中,G=[BΨH,I]表示在稀疏基下新的测量矩阵,x=[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量,Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示;
根据待求稀疏域下的后向散射功率和噪声矢量x是稀疏这一先验信息,定义其服从高斯分布,即x~N(0,Λ);其中Λ=diag(Γ),为稀疏贝叶斯学习模型的超参数矢量,表征了x的空间分布;前D个超参数τ1,...τD控制高度向上D个散射体后向散射功率的稀疏性,后N2个超参数为***噪声;则Y关于Γ的概率密度函数为:
p(Y|Γ)=∫p(Y|x)p(x|Γ)dx;
步骤b,参数初始值定义:
设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件;
步骤c,计算后验概率分布:
基于贝叶斯准则和步骤a定义的先验分布,计算x的均值矢量和协方差矩阵:
ΣY=V+GΛ(i)GH
式中,上标i表示当前迭代次数,μ为x的均值矢量,Σx和ΣY分别为x和Y的协方差矩阵,V为样本协方差矩阵估计误差δ的二阶统计;
步骤d,基于步骤c得到的均值矢量和协方差矩阵,采用EM算法学习更新超参数向量:
步骤e,重复步骤c、d,直至满足迭代终止条件,使用超参数向量Γ中的前D个超参数Γ(1…D)计算冠层散射和地表散射的层析谱:
4.根据权利要求3所述的方法,其特征在于,迭代终止条件为||Γi+1-Γi||2/||Γi||2≤tol,和/或达到最大迭代次数。
5.根据权利要求1所述的方法,其特征在于,根据地表散射的层析谱提取森林反演区域的林下地形的方法为:提取地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。
6.根据权利要求3所述的方法,其特征在于,根据冠层散射的层析谱提取森林反演区域的植被高度的方法为:通过功率衰减法确定冠层顶部的位置,即:以冠层相位中心处为起始点,以预设的步长向上衰减,当功率衰减处提取的高度值和Lidar数据提取的冠层高度DSM的均方根误差达到最小时,将此时的高度值作为冠层高度的估计值;通过冠层高度减去地面高度即为植被高度。
7.一种基于散射机制分解的植被垂直结构及林下地形反演装置,其特征在于,包括:
数据预处理模块,用于:获取森林反演区域多基线多极化的SAR数据集,并进行配准、去平地,得到极化层析SAR数据;
协方差计算分解模块,用于:对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵,并分解计算冠层散射和地表散射的干涉协方差矩阵;
其中,多基线多极化协方差矩阵的计算方法与分解方法为:
式中,y为步骤1得到的极化层析SAR数据,gHV=gVH,gHH、gHV、gVH、gVV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值;(*)H为共轭转置操作,E[·]为求期望运算;W为计算得到的多基线多极化协方差矩阵;Cg和Rg分别为地表散射的极化协方差矩阵和干涉协方差矩阵,Cv和Rv分别为冠层散射的极化协方差矩阵和干涉协方差矩阵,为克罗内克积;
对多基线多极化协方差矩阵W进行SVD分解,使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg:
首先,取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R1和R2,使用R1和R2表示冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg为:
式中,a、b为一对实数;
然后,基于冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg的半正定约束,即Rv和Rg的特征值非负,求解实数a、b的取值区间;
再后,根据以下散射机制最大化分离的准则,求解实数a、b:
其中,trace()为求矩阵的迹,||·||F为矩阵的Frobenius范数;
最后,根据求解得到的实数a、b和矩阵R1和R2,即可得到冠层散射的干涉协方差矩阵Rv和地表散射的干涉协方差矩阵Rg;
模型构建模块,用于:针对冠层散射和地表散射的干涉协方差矩阵,均进行矢量化和稀疏表达,分别得到冠层散射和地表散射的稀疏贝叶斯学习模型;
层析谱计算模块,用于:采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数,进而获得冠层散射和地表散射的层析谱;
植被垂直结构及林下地形反演模块,用于:根据地表散射的层析谱提取森林反演区域的林下地形,根据冠层散射的层析谱提取森林反演区域的植被高度。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011164833.3A CN112363161B (zh) | 2020-10-27 | 2020-10-27 | 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011164833.3A CN112363161B (zh) | 2020-10-27 | 2020-10-27 | 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112363161A CN112363161A (zh) | 2021-02-12 |
CN112363161B true CN112363161B (zh) | 2022-12-20 |
Family
ID=74510786
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011164833.3A Active CN112363161B (zh) | 2020-10-27 | 2020-10-27 | 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112363161B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113466857B (zh) * | 2021-05-11 | 2022-11-04 | 中国地质大学(武汉) | 基于非局部平均的TomoSAR林下地形反演方法及*** |
CN114545407B (zh) * | 2021-09-24 | 2023-03-24 | 中国科学院精密测量科学与技术创新研究院 | 一种基于分布式压缩感知的星载差分层析sar成像方法 |
CN114252480A (zh) * | 2021-12-21 | 2022-03-29 | 中南大学 | 基于单极化sar数据邻域像素反演裸土区土壤湿度方法 |
CN115166741B (zh) * | 2022-09-08 | 2022-11-29 | 中国科学院空天信息创新研究院 | 一种基于简化模型的双相位中心极化层析分解方法 |
CN115166739B (zh) * | 2022-09-08 | 2022-11-29 | 中国科学院空天信息创新研究院 | 一种基于多基线层析极化目标分解的目标高度估计方法 |
CN115343712B (zh) * | 2022-10-18 | 2022-12-20 | 中国电子科技集团公司第十四研究所 | 一种反演植被高程的高低频极化干涉试验*** |
CN117452432B (zh) * | 2023-12-21 | 2024-03-15 | 西南林业大学 | 一种基于森林穿透补偿的森林冠层高度估测方法 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2784537A1 (en) * | 2013-05-15 | 2014-10-01 | Institute of Electronics, Chinese Academy of Sciences | Inversion method and apparatus based on polarimetric interferometric synthetic aperture radar |
CN105005047A (zh) * | 2015-07-17 | 2015-10-28 | 武汉大学 | 后向散射优化的森林复杂地形校正及树高反演方法、*** |
KR20160050322A (ko) * | 2014-10-29 | 2016-05-11 | 한국해양과학기술원 | 타깃 탐지의 정확성 개선 방법 |
CN105974412A (zh) * | 2016-06-07 | 2016-09-28 | 电子科技大学 | 一种用于合成孔径雷达的目标特征提取方法 |
CN110109111A (zh) * | 2019-04-28 | 2019-08-09 | 西安电子科技大学 | 极化干涉sar稀疏植被高度反演方法 |
CN110161499A (zh) * | 2019-05-09 | 2019-08-23 | 东南大学 | 改进的稀疏贝叶斯学习isar成像散射系数估计方法 |
CN110569624A (zh) * | 2019-09-20 | 2019-12-13 | 哈尔滨工业大学 | 适用于PolInSAR反演的森林三层散射模型的确定及分析方法 |
CN110703220A (zh) * | 2019-10-12 | 2020-01-17 | 中南大学 | 一种顾及时间去相干因子的多基线PolInSAR植被参数反演方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US9417323B2 (en) * | 2012-11-07 | 2016-08-16 | Neva Ridge Technologies | SAR point cloud generation system |
CN108983229B (zh) * | 2018-05-03 | 2022-04-19 | 电子科技大学 | 基于sar层析技术的高压输电铁塔高度及形变提取方法 |
-
2020
- 2020-10-27 CN CN202011164833.3A patent/CN112363161B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2784537A1 (en) * | 2013-05-15 | 2014-10-01 | Institute of Electronics, Chinese Academy of Sciences | Inversion method and apparatus based on polarimetric interferometric synthetic aperture radar |
KR20160050322A (ko) * | 2014-10-29 | 2016-05-11 | 한국해양과학기술원 | 타깃 탐지의 정확성 개선 방법 |
CN105005047A (zh) * | 2015-07-17 | 2015-10-28 | 武汉大学 | 后向散射优化的森林复杂地形校正及树高反演方法、*** |
CN105974412A (zh) * | 2016-06-07 | 2016-09-28 | 电子科技大学 | 一种用于合成孔径雷达的目标特征提取方法 |
CN110109111A (zh) * | 2019-04-28 | 2019-08-09 | 西安电子科技大学 | 极化干涉sar稀疏植被高度反演方法 |
CN110161499A (zh) * | 2019-05-09 | 2019-08-23 | 东南大学 | 改进的稀疏贝叶斯学习isar成像散射系数估计方法 |
CN110569624A (zh) * | 2019-09-20 | 2019-12-13 | 哈尔滨工业大学 | 适用于PolInSAR反演的森林三层散射模型的确定及分析方法 |
CN110703220A (zh) * | 2019-10-12 | 2020-01-17 | 中南大学 | 一种顾及时间去相干因子的多基线PolInSAR植被参数反演方法 |
Non-Patent Citations (5)
Title |
---|
The Impact of Forest Density on Forest Height Inversion Modeling from Polarimetric InSAR Data;Wang, Changcheng et al.;《Remote Sensing》;20161231;全文 * |
基于Freeman-Durden分解理论的植被高度反演;宋桂萍等;《赤峰学院学报(自然科学版)》;20190525(第05期);全文 * |
基于散射机制分解的ESPRIT植被高度反演新方法;宋桂萍等;《现代测绘》;20130925(第05期);全文 * |
干涉、极化干涉SAR技术森林高度估测算法研究进展;张王菲等;《遥感技术与应用》;20171215(第06期);全文 * |
极化干涉SAR相干最优理论及其验证分析;尚玉双等;《测绘与空间地理信息》;20140125(第01期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN112363161A (zh) | 2021-02-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112363161B (zh) | 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 | |
CN108983229B (zh) | 基于sar层析技术的高压输电铁塔高度及形变提取方法 | |
Cazcarra-Bes et al. | Comparison of tomographic SAR reflectivity reconstruction algorithms for forest applications at L-band | |
Pannekoucke et al. | Background‐error correlation length‐scale estimates and their sampling statistics | |
CN103954950A (zh) | 一种基于样本协方差矩阵稀疏性的波达方向估计方法 | |
CN110703249B (zh) | 稳健高效合成孔径雷达多元特征增强成像方法 | |
CN107085206B (zh) | 一种基于自适应稀疏保持投影的一维距离像识别方法 | |
CN111766577B (zh) | 一种基于三阶段算法p波段的输电线路通道树木高度反演方法 | |
CN110113085A (zh) | 一种基于协方差矩阵重构的波束形成方法及*** | |
Wipf et al. | Beamforming using the relevance vector machine | |
Wu et al. | Superresolution radar imaging via peak search and compressed sensing | |
Berenger et al. | A deep learning approach for sar tomographic imaging of forested areas | |
CN113421198B (zh) | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 | |
CN110940638A (zh) | 一种高光谱图像亚像元级水体边界探测方法及探测*** | |
CN108828482B (zh) | 结合稀疏和低秩特性的欠采样磁共振扩散谱的重建方法 | |
CN110231625B (zh) | 一种基于多尺度融合的综合孔径成像方法 | |
CN113406560B (zh) | 一种非相干分布宽带源的角度和频率参数估计方法 | |
CN113640891B (zh) | 一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法 | |
Gómez-Dans et al. | Indoor C-band polarimetric interferometry observations of a mature wheat canopy | |
Dou et al. | Deep learning imaging for 1-D aperture synthesis radiometers | |
CN109633635B (zh) | 基于结构化递归最小二乘的米波雷达测高方法 | |
CN108830799B (zh) | 基于相对极化全变差的极化sar影像相干斑抑制方法 | |
CN113435366A (zh) | 一种小波域的多时高光谱图像贝叶斯解混方法 | |
CN111899226A (zh) | 一种基于多任务稀疏学习的高光谱影像目标先验优化方法 | |
CN113466857B (zh) | 基于非局部平均的TomoSAR林下地形反演方法及*** |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |