CN110133657A - 基于PolInSAR地形效应补偿的植被参数反演方法 - Google Patents

基于PolInSAR地形效应补偿的植被参数反演方法 Download PDF

Info

Publication number
CN110133657A
CN110133657A CN201910590564.8A CN201910590564A CN110133657A CN 110133657 A CN110133657 A CN 110133657A CN 201910590564 A CN201910590564 A CN 201910590564A CN 110133657 A CN110133657 A CN 110133657A
Authority
CN
China
Prior art keywords
polinsar
matrix
lut
vegetation
coherence
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.)
Pending
Application number
CN201910590564.8A
Other languages
English (en)
Inventor
索志勇
薛超
李真芳
赵秉吉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Publication of CN110133657A publication Critical patent/CN110133657A/zh
Pending legal-status Critical Current

Links

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/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
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • 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
    • G01S13/904—SAR modes
    • G01S13/9076—Polarimetric features in SAR
    • 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
    • G01S13/9094—Theoretical aspects

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)
  • Image Processing (AREA)

Abstract

本发明提出了基于PolInSAR地形效应补偿的植被参数反演方法，通过对PolInSAR极化方位角和地形效应补偿，旨在提高植被参数反演的精度。实现步骤为：对PolInSAR干涉相位图进行去平地；获取PolInSAR主图像的Pauli基矢量矩阵和辅图像的Pauli基矢量矩阵；对Pauli基矢量矩阵进行POA校正；获取PolInSAR的极化自相干矩阵集合和极化互相干矩阵集合；对自相干矩阵集合以及互相干矩阵集合进行地形干涉相位坡度补偿；获取PolInSAR的地表相干系数矩阵和体散射相干系数矩阵；对植被场景的参数进行反演。本发明降低了极化相干矩阵之间的地形干涉相位坡度的差异性，提高了植被参数反演的精度。

Description

基于PolInSAR地形效应补偿的植被参数反演方法
技术领域
本发明属于信号处理技术领域，涉及一种植被参数反演方法，具体涉及一种基于PolInSAR地形效应补偿的植被参数反演方法，可用于森林的分类和制图、大面积植被检测以及复杂地形高精度三维地形测绘等。
背景技术
合成孔径雷达SAR是一种全天时、全天候的空间微波遥感成像雷达，能同时实现对地面目标的距离向和方位向的高分辨率成像。干涉合成孔径雷达InSAR利用不同天线的回波数据进行干涉处理，可以对地面的高程进行估计，对海流进行测高和测速等。极化合成孔径雷达PolSAR通过对目标进行全极化测量，实现目标的地物分类，荒漠化评估等。极化干涉合成孔径雷达PolInSAR集PolSAR和InSAR测量技术于一体，可以把目标的精细物理特征与空间分布特性结合起来，能进行植被参数反演，岩层分析，目标的检测和识别等。
对于植被场景，PolInSAR图像分辨单元内存在来自地面、植被茎干、植被冠层等多种极化散射回波的组合，从而使分辨单元内存在地面、植被茎干、植被冠层等多种极化相干系数，这些极化相干系数包含植被的参数信息，通过极化相干系数与植被参数之间关系，能反演植被的参数。因此，准确地获取这些极化相干系数是PolInSAR精确反演植被参数的关键。通过PolInSAR回波的散射矩阵来构建相干矩阵，并利用相干矩阵来估计植被的体散射相干系数和地表相干系数，进而实现对植被参数的反演。
当植被场景存在地形坡度时，PolInSAR回波的极化状态与发射波相比会发生变化，导致散射矩阵的幅度和相位发生变化，即散射矩阵发生旋转，无法准确地描述目标的散射特性。此外，由于地形坡度引起的相邻位置地表存在水平高度差异，导致相邻的处理像素之间的干涉相位存在差异，造成多视处理过程中极化相干矩阵样本之间的差异性，从而降低了用极化相干矩阵来估计植被体散射相干系数和地表相干系数的准确度，进而降低植被参数反演的精度。
目前，利用PolInSAR来对带地形坡度的植被场景进行植被参数反演的方法中，应用最广泛的是基于PolInSAR极化方位角(POA)补偿的植被参数反演方法，例如：文献“龙江平,丁晓利,汪长城.极化方位角补偿信息支持下的植被参数反演[J].测绘学报,2014,43(10):1051-1060.”公开的一种极化方位角补偿信息支持下的植被参数反演方法，该方法对PolInSAR干涉相位图进行去平地，获取PolInSAR的Pauli基矢量矩阵，对Pauli基矢量矩阵进行POA校正，获取PolInSAR的极化相干矩阵集合，获取PolInSAR的地表相干系数矩阵和体散射相干系数矩阵，对植被场景的参数进行反演。该方法对Pauli基矢量矩阵进行POA校正，消除了Pauli矢量矩阵里的每个Pauli矢量的旋转效应，提高了Pauli矢量矩阵描述植被散射特性的准确度，从而提高了由Pauli矢量矩阵得到的极化相干矩阵集合描述植被相干特性的准确度，但忽略了由PolInSAR地形效应引起的相邻极化相干矩阵之间干涉相位坡度的差异性，降低了用极化相干矩阵来估计植被体散射相干系数和地表相干系数的准确度，使得植被参数反演的精度较低。
发明内容
本发明的目的在于克服现有技术的不足，提出一种基于PolInSAR地形效应补偿的植被参数反演方法，通过对PolInSAR极化方位角和地形效应补偿，旨在提高植被参数反演的精度。
为实现上述目的，本发明采取的技术方案包括如下步骤：
(1)对PolInSAR的干涉相位图进行去平地：
(1a)对PolInSAR包含M×N个像元的干涉相位图的平地干涉相位进行估计，得到与干涉相位图大小相同的PolInSAR平地干涉相位图，其中M≥2，N≥2；
(1b)对PolInSAR平地干涉相位图进行复数域转换，得到M×N个相位修正因子，并通过每个相位修正因子对PolInSAR辅图像对应位置的像元进行修正，得到修正后的PolInSAR辅图像；
(1c)将PolInSAR主图像与修正后的PolInSAR辅图像进行共轭内积，得到去平地后的PolInSAR干涉相位图；
(2)获取PolInSAR主图像的Pauli基矢量矩阵和辅图像的Pauli基矢量矩阵：
计算PolInSAR主图像上每个像元的Pauli基矢量k1，得到包含M×N个k1的Pauli基矢量矩阵{k1}，同时计算修正后的PolInSAR辅图像上每个像元的Pauli基矢量k2，得到包含M×N个k2的Pauli基矢量矩阵{k2}；
(3)对Pauli基矢量矩阵{k1}和{k2}进行POA校正：
(3a)计算PolInSAR主图像上每个像元对应的主天线POA值βm，得到包含M×N个βm的主天线POA值矩阵{βm}，同时计算PolInSAR辅图像上每个像元对应的辅天线POA值βs，得到包含M×N个βs的辅天线POA值矩阵{βs}；
(3b)对{βm}中每一个βm进行极化基变换，得到包含M×N个旋转矩阵Gm的矩阵{Gm}，并通过每个Gm对{k1}中相同位置的k1进行校正，得到校正后的PolInSAR主图像Pauli基矢量矩阵{k′1}，同时对{βs}中每一个βs进行极化基变换，得到包含M×N个旋转矩阵Gs的矩阵{Gs}，并通过每个Gs对{k2}中相同位置的k2进行校正，得到校正后的PolInSAR辅图像Pauli基矢量矩阵{k′2}；
(4)获取PolInSAR的极化自相干矩阵集合和极化互相干矩阵集合：
(4a)获取PolInSAR的极化自相干矩阵集合：
对{k′1}中的每个k′1与k′1的共轭转置矢量(k′1)*T进行外积，得到包含M×N个PolInSAR主图像极化自相干矩阵T11的集合{T11}，同时对{k′2}中的每个k'2与k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR辅图像极化自相干矩阵T22的集合{T22}；
(4b)获取PolInSAR的极化互相干矩阵集合：
对{k′1}中的每个k′1与{k'2}中与每个k′1相同位置的k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR极化互相干矩阵Ω12的集合{Ω12}；
(5)对自相干矩阵集合{T11}和{T22}，以及互相干矩阵集合{Ω12}进行地形干涉相位坡度补偿：
(5a)计算去平地后的PolInSAR干涉相位图包含的每一像元的方位向干涉相位坡度ωa和距离向干涉相位坡度ωr，得到方位向干涉相位坡度矩阵{ωa}和距离向干涉相位坡度矩阵{ωr}；
(5b)通过{ωa}和{ωr}对{T11}、{T22}和{Ω12}中的每一个矩阵进行地形干涉相位坡度补偿，得到补偿后的PolInSAR极化自相干矩阵集合{T′11}和{T′22}，以及互相干矩阵集合{Ω'12}；
(6)获取PolInSAR的地表相干系数矩阵和体散射相干系数矩阵：
(6a)设置含有K个不同随机散射机理矢量W的集合{W}，其中，K为正整数，且500≤K≤10000；
(6b)对{T′11}、{T′22}和{Ω'12}进行组合，得到极化相干矩阵组合的集合{(T′11,T′22,Ω'12)}，并通过每个(T′11,T′22,Ω'12)与{W}计算植被的随机相干系数，得到包含M×N个植被的随机相干系数分布区域Aγ的集合{Aγ}，然后对每个Aγ进行直线拟合，得到包含M×N条直线Lγ的集合{Lγ}，其中，每个Aγ中均含有K个植被的随机相干系数；
(6c)对{Lγ}中的每个Lγ与单位圆进行相交，得到包含M×N个交点对(γuc1,γuc2)的集合{(γuc1,γuc2)}，同时对{Lγ}中的每个Lγ与{Aγ}中对应的Aγ的进行相交，得到包含M×N个交点对(γcre1,γcre2)的集合{(γcre1,γcre2)}；
(6d)对{(γuc1,γuc2)}中的每个(γuc1,γuc2)进行地表散射判别，得到包含M×N个PolIInSAR地表相干系数γg的矩阵{γg}，并对{(γcre1,γcre2)}中的每个(γcre1,γcre2)进行体散射判别，得到包含M×N个PolInSAR体散射相干系数γv的矩阵{γv}；
(7)对植被场景的参数进行反演：
(7a)对{γv}与{γg}进行共轭内积，得到去地表后的PolInSAR体散射相干系数矩阵{γ′v}；
(7b)设置包含U×V个植被高度hlut的查找矩阵{hlut}和包含U×V个消光系数σlut的查找矩阵{σlut}，并对{hlut}和{σlut}进行组合，得到植被参数组合的矩阵{(hlut,σlut)}，其中，400≤U≤2000，400≤V≤2000；
(7c)获取PolInSAR视角θ，并计算PolInSAR有效垂直波数kz，然后通过kz、θ和{(hlut,σlut)}中的每个(hlut,σlut)计算植被的体散射相干系数γlut，得到包含U×V个植被体散射相干系数γlut的查找矩阵{γlut}；
(7d)在{γlut}里获取与{γv'}中的每个γv'最为相似的相干系数的二维矩阵坐标，得到包含M×N个二维矩阵坐标pos的矩阵{pos}，并在{hlut}里提取坐标与{pos}中的每个pos相同的植被高度，得到包含M×N个植被高度的矩阵Htree。
本发明与现有技术相比，具有如下优点：
本发明在对植被场景的参数进行反演的过程中，首先利用极化方位角对PolInSAR的Pauli矢量矩阵进行补偿，实现PolInSAR极化方位角补偿，并利用地形干涉相位坡度矩阵对PolInSAR的极化相干矩阵集合进行补偿，实现PolInSAR的地形效应补偿，降低了现有技术中仅用极化方位角对PolInSAR的Pauli矢量矩阵进行补偿具有的极化相干矩阵之间地形干涉相位坡度的差异性，地表相干系数和体散射相干系数能更准确描述植被的参数特征，提高了植被参数反演的精度。
附图说明
图1为本发明的实现流程图；
图2为本发明对植被场景进行仿真获取的干涉相位图；
图3为本发明和现有技术的植被参数反演结果的植被高度分布对比图。
具体实施方式
下面结合附图和具体实施例，对本发明作进一步详细描述：
参照图1，基于PolInSAR地形效应补偿的植被参数反演方法，包括如下步骤：
步骤1)对PolInSAR的干涉相位图进行去平地：
步骤1a)对PolInSAR包含M×N个像元的干涉相位图的平地干涉相位进行估计，得到与干涉相位图大小相同的PolInSAR平地干涉相位图，且每个像元的平地干涉相位的估计公式为：
其中，φflat为像元的平地干涉相位，λ为PolInSAR的工作波长，Rflat1为像元对应的平地与PolInSAR主天线的斜距，Rflat2为像元对应的平地与PolInSAR辅天线的斜距，M≥2，N≥2；限定M≥2，N≥2主要是考虑后续的对相干矩阵的多视处理，要使多视处理后的相干矩阵的秩为3，则至少需要3个相干矩阵样本，而PolInSAR图像是矩阵形式，所以图像须至少有2×2个像元；
步骤1b)对PolInSAR平地干涉相位图进行复数域转换，得到M×N个相位修正因子，并通过每个相位修正因子对PolInSAR辅图像对应位置的像元进行修正，得到修正后的PolInSAR辅图像；
步骤1c)将PolInSAR主图像与修正后的PolInSAR辅图像进行共轭内积，得到去平地后的PolInSAR干涉相位图；
步骤2)获取PolInSAR主图像的Pauli基矢量矩阵和辅图像的Pauli基矢量矩阵：
计算PolInSAR主图像上每个像元的Pauli基矢量k1，得到包含M×N个k1的Pauli基矢量矩阵{k1}，同时计算修正后的PolInSAR辅图像上每个像元的Pauli基矢量k2，得到包含M×N个k2的Pauli基矢量矩阵{k2}；
步骤3)对Pauli基矢量矩阵{k1}和{k2}进行POA校正：
步骤3a)计算PolInSAR主图像上每个像元对应的主天线POA值βm，得到包含M×N个βm的主天线POA值矩阵{βm}，同时计算PolInSAR辅图像上每个像元对应的辅天线POA值βs，得到包含M×N个βs的辅天线POA值矩阵{βs}，对于每个βm和βs，计算公式分别为：
其中，ω为像元的方位向地形坡度角，γ为像元的距离向地形坡度角，θm为PolInSAR主天线照射到像元的视角，θs为PolInSAR雷达辅天线照射到像元的视角。
步骤3b)对{βm}中每一个βm进行极化基变换，得到包含M×N个旋转矩阵Gm的矩阵{Gm}，并通过每个Gm对{k1}中相同位置的k1进行校正，得到校正后的PolInSAR主图像Pauli基矢量矩阵{k′1}，同时对{βs}中每一个βs进行极化基变换，得到包含M×N个旋转矩阵Gs的矩阵{Gs}，并通过每个Gs对{k2}中相同位置的k2进行校正，得到校正后的PolInSAR辅图像Pauli基矢量矩阵{k'2}；
步骤4)获取PolInSAR的极化自相干矩阵集合和极化互相干矩阵集合：
步骤4a)获取PolInSAR的极化自相干矩阵集合：
对{k′1}中的每个k′1与k′1的共轭转置矢量(k′1)*T进行外积，得到包含M×N个PolInSAR主图像极化自相干矩阵T11的集合{T11}，同时对{k'2}中的每个k'2与k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR辅图像极化自相干矩阵T22的集合{T22}；
步骤4b)获取PolInSAR的极化互相干矩阵集合：
对{k′1}中的每个k′1与{k'2}中与每个k′1相同位置的k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR极化互相干矩阵Ω12的集合{Ω12}；
步骤5)对自相干矩阵集合{T11}和{T22}，以及互相干矩阵集合{Ω12}进行地形干涉相位坡度补偿：
步骤5a)计算去平地后的PolInSAR干涉相位图包含的每一像元的方位向干涉相位坡度ωa和距离向干涉相位坡度ωr，得到方位向干涉相位坡度矩阵{ωa}和距离向干涉相位坡度矩阵{ωr}；
步骤5b)通过{ωa}和{ωr}对{T11}、{T22}和{Ω12}中的每一个矩阵进行地形干涉相位坡度补偿，得到补偿后的PolInSAR极化自相干矩阵集合{T′11}和{T′22}，以及互相干矩阵集合{Ω'12}，对于每次的地形干涉相位坡度补偿操作，补偿公式为：
其中，m和n均为非负整数，且0≤m＜M，0≤n＜N，{T11}为PolInSAR主图像极化自相干矩阵的集合，{T22}为PolInSAR辅图像极化自相干矩阵的集合，{Ω12}为PolInSAR极化互相干矩阵的集合，T11(x+k,y+l)为{T11}中第x+k行y+l列的元素，T22(x+k,y+l)为{T22}中第x+k行y+l列的元素，Ω12(x+k,y+l)为{Ω12}中第x+k行y+l列的元素，{ωa}为方位向干涉相位坡度矩阵，{ωr}为距离向干涉相位坡度矩阵，ωa(x,y)为{ωa}中第x行y列的元素，ωr(x,y)为{ωr}中第x行y列的元素，且m+1≤x≤M-m，n+1≤x≤N-n。
步骤6)获取PolInSAR的地表相干系数矩阵和体散射相干系数矩阵：
步骤6a)设置含有K个不同随机散射机理矢量W的集合{W}，其中，K为正整数，且500≤K≤10000；限定500≤K≤10000是根据多次实验的结果，K小于500时，导致后续得到随机散射相干系数的个数较少，无法较好地反映散射分布，而K大于10000时，会影响后续处理的效率；
步骤6b)对{T′11}、{T′22}和{Ω'12}进行组合，得到极化相干矩阵组合的集合{(T′11,T′22,Ω'12)}，并通过每个(T′11,T′22,Ω'12)与{W}计算植被的随机相干系数，得到包含M×N个植被的随机相干系数分布区域Aγ的集合{Aγ}，然后对每个Aγ进行直线拟合，得到包含M×N条直线Lγ的集合{Lγ}，其中，每个Aγ中均含有K个植被的随机相干系数；
步骤6c)对{Lγ}中的每个Lγ与复数域单位圆进行相交，得到包含M×N个交点对(γuc1,γuc2)的集合{(γuc1,γuc2)}，同时对{Lγ}中的每个Lγ与{Aγ}中对应的Aγ的进行相交，得到包含M×N个交点对(γcre1,γcre2)的集合{(γcre1,γcre2)}；
步骤6d)对{(γuc1,γuc2)}中的每个(γuc1,γuc2)进行地表散射判别，得到包含M×N个PolIInSAR地表相干系数γg的矩阵{γg}，并对{(γcre1,γcre2)}中的每个(γcre1,γcre2)进行体散射判别，得到包含M×N个PolInSAR体散射相干系数γv的矩阵{γv}，对于每次的地表散射判别和体散射判别，采用差异最大化方法，实现步骤为：
步骤6d1)计算HH+VV极化相干系数γHH+VV、HH-VV极化相干系数γHH-VV和HV极化相干系数γHV；
步骤6d2)判断γuc1的相位是否小于γHV的相位，若是，则将γuc1作为PolInSAR地表相干系数γg，否则，将γuc2作为PolInSAR地表相干系数γg；
步骤6d3)计算γcre1与γg的距离Rcre1，γcre2与γg的距离Rcre2，γHH-VV与γg的距离RHH+VV，γHH-VV与γg的距离RHH-VV；
步骤6d4)判断Rcre1是否大于RHH+VV，且Rcre1是否大于RHH-VV，若是，则将γcre1作为PolInSAR体散射相干系数γv，否则，则将γcre2作为PolInSAR体散射相干系数γv。
步骤7)对植被场景的参数进行反演：
步骤7a)对{γv}与{γg}进行共轭内积，得到去地表后的PolInSAR体散射相干系数矩阵{γ′v}；
步骤7b)设置包含U×V个植被高度hlut的查找矩阵{hlut}和包含U×V个消光系数σlut的查找矩阵{σlut}，并对{hlut}和{σlut}进行组合，得到植被参数组合的矩阵{(hlut,σlut)}，其中，400≤U≤2000，400≤V≤2000；限定400≤U≤2000，400≤V≤2000主要是根据多次实验的结果，当U和V小于400时，会降低植被参数的精度，当U和V小于400时，会降低植被参数的查找效率；
步骤7c)获取PolInSAR视角θ，并计算PolInSAR有效垂直波数kz，然后通过kz、θ和{(hlut,σlut)}中的每个(hlut,σlut)计算植被的体散射相干系数γlut，得到包含U×V个植被体散射相干系数γlut的查找矩阵{γlut}，且对于每个植被的体散射相干系数的计算操作，计算公式为：
其中，kz为有效垂直波数，θ为PolInSAR视角，α为植被场景距离向地形坡度，hlut为植被高度，σlut为消光系数。
步骤7d)在{γlut}里获取与{γ′v}中的每个γ′v最为相似的相干系数的二维矩阵坐标，得到包含M×N个二维矩阵坐标pos的矩阵{pos}，并在{hlut}里提取坐标与{pos}中的每个pos相同的植被高度，得到包含M×N个植被高度的矩阵Htree。
以下通过仿真实验，对本发明的技术效果作进一步说明：
1、仿真条件和内容：
本实验所用的仿真数据为ESA PolSARPro仿真软件生成的带地形坡度的植被场景PolInSAR全极化数据，所设置的仿真参数如表1所示：
表1PolSARPro仿真参数
仿真完成后，在MATLAB2015b软件上读取仿真全极化数据，得到HH通道受干涉相位图如图2所示。
对现有技术和本发明得到的植被场景高度图进行统计，得到植被高度统计对比直方图如图3所示，分别求取现有技术和本发明植被场景高度图的高度均值和高度标准差，如表2所示：
表2植被参数反演结果
2、仿真结果分析
从表1可知，所设置的仿真参数中植被的高度为11m。从表2可知，现有技术的植被参数反演结果中，植被高度为9.9499m，标准差为1.0907m，而本发明的植被参数反演结果中，植被高度为9.9742m，标准差为1.0647m；本发明反演的植被高度更接近真实植被高度。并且，由图3可知，相对于现有技术，本发明得到的靠近真实树高的植被高度点数更多。因此，本发明对带地形坡度的植被场景的参数反演精度比现有技术高。

Claims (6)

1.一种基于PolInSAR地形效应补偿的植被参数反演方法，其特征在于，包括如下步骤：
(1)对PolInSAR的干涉相位图进行去平地：
(1a)对PolInSAR包含M×N个像元的干涉相位图的平地干涉相位进行估计，得到与干涉相位图大小相同的PolInSAR平地干涉相位图，其中M≥2，N≥2；
(1b)对PolInSAR平地干涉相位图进行复数域转换，得到M×N个相位修正因子，并通过每个相位修正因子对PolInSAR辅图像对应位置的像元进行修正，得到修正后的PolInSAR辅图像；
(1c)将PolInSAR主图像与修正后的PolInSAR辅图像进行共轭内积，得到去平地后的PolInSAR干涉相位图；
(2)获取PolInSAR主图像的Pauli基矢量矩阵和辅图像的Pauli基矢量矩阵：
计算PolInSAR主图像上每个像元的Pauli基矢量k1，得到包含M×N个k1的Pauli基矢量矩阵{k1}，同时计算修正后的PolInSAR辅图像上每个像元的Pauli基矢量k2，得到包含M×N个k2的Pauli基矢量矩阵{k2}；
(3)对Pauli基矢量矩阵{k1}和{k2}进行POA校正：
(3a)计算PolInSAR主图像上每个像元对应的主天线POA值βm，得到包含M×N个βm的主天线POA值矩阵{βm}，同时计算PolInSAR辅图像上每个像元对应的辅天线POA值βs，得到包含M×N个βs的辅天线POA值矩阵{βs}；
(3b)对{βm}中每一个βm进行极化基变换，得到包含M×N个旋转矩阵Gm的矩阵{Gm}，并通过每个Gm对{k1}中相同位置的k1进行校正，得到校正后的PolInSAR主图像Pauli基矢量矩阵{k’1}，同时对{βs}中每一个βs进行极化基变换，得到包含M×N个旋转矩阵Gs的矩阵{Gs}，并通过每个Gs对{k2}中相同位置的k2进行校正，得到校正后的PolInSAR辅图像Pauli基矢量矩阵{k'2}；
(4)获取PolInSAR的极化自相干矩阵集合和极化互相干矩阵集合：
(4a)获取PolInSAR的极化自相干矩阵集合：
对{k’1}中的每个k’1与k’1的共轭转置矢量(k’1)*T进行外积，得到包含M×N个PolInSAR主图像极化自相干矩阵T11的集合{T11}，同时对{k'2}中的每个k'2与k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR辅图像极化自相干矩阵T22的集合{T22}；
(4b)获取PolInSAR的极化互相干矩阵集合：
对{k’1}中的每个k’1与{k'2}中与每个k’1相同位置的k'2的共轭转置矢量(k'2)*T进行外积，得到包含M×N个PolInSAR极化互相干矩阵Ω12的集合{Ω12}；
(5)对自相干矩阵集合{T11}和{T22}，以及互相干矩阵集合{Ω12}进行地形干涉相位坡度补偿：
(5a)计算去平地后的PolInSAR干涉相位图包含的每一像元的方位向干涉相位坡度ωa和距离向干涉相位坡度ωr，得到方位向干涉相位坡度矩阵{ωa}和距离向干涉相位坡度矩阵{ωr}；
(5b)通过{ωa}和{ωr}对{T11}、{T22}和{Ω12}中的每一个矩阵进行地形干涉相位坡度补偿，得到补偿后的PolInSAR极化自相干矩阵集合{T’11}和{T’22}，以及互相干矩阵集合{Ω'12}；
(6)获取PolInSAR的地表相干系数矩阵和体散射相干系数矩阵：
(6a)设置含有K个不同随机散射机理矢量W的集合{W}，其中，K为正整数，且500≤K≤10000；
(6b)对{T’11}、{T’22}和{Ω'12}进行组合，得到极化相干矩阵组合的集合{(T’11,T’22,Ω'12)}，并通过每个(T’11,T’22,Ω'12)与{W}计算植被的随机相干系数，得到包含M×N个植被的随机相干系数分布区域Aγ的集合{Aγ}，然后对每个Aγ进行直线拟合，得到包含M×N条直线Lγ的集合{Lγ}，其中，每个Aγ中均含有K个植被的随机相干系数；
(6c)对{Lγ}中的每个Lγ与单位圆进行相交，得到包含M×N个交点对(γuc1,γuc2)的集合{(γuc1,γuc2)}，同时对{Lγ}中的每个Lγ与{Aγ}中对应的Aγ的进行相交，得到包含M×N个交点对(γcre1,γcre2)的集合{(γcre1,γcre2)}；
(6d)对{(γuc1,γuc2)}中的每个(γuc1,γuc2)进行地表散射判别，得到包含M×N个PolIInSAR地表相干系数γg的矩阵{γg}，并对{(γcre1,γcre2)}中的每个(γcre1,γcre2)进行体散射判别，得到包含M×N个PolInSAR体散射相干系数γv的矩阵{γv}；
(7)对植被场景的参数进行反演：
(7a)对{γv}与{γg}进行共轭内积，得到去地表后的PolInSAR体散射相干系数矩阵{γ’v}；
(7b)设置包含U×V个植被高度hlut的查找矩阵{hlut}和包含U×V个消光系数σlut的查找矩阵{σlut}，并对{hlut}和{σlut}进行组合，得到植被参数组合的矩阵{(hlut,σlut)}，其中，400≤U≤2000，400≤V≤2000；
(7c)获取PolInSAR视角θ，并计算PolInSAR有效垂直波数kz，然后通过kz、θ和{(hlut,σlut)}中的每个(hlut,σlut)计算植被的体散射相干系数γlut，得到包含U×V个植被体散射相干系数γlut的查找矩阵{γlut}；
(7d)在{γlut}里获取与{γ’v}中的每个γ’v最为相似的相干系数的二维矩阵坐标，得到包含M×N个二维矩阵坐标pos的矩阵{pos}，并在{hlut}里提取坐标与{pos}中的每个pos相同的植被高度，得到包含M×N个植被高度的矩阵Htree。
2.根据权利要求1所述的基于PolInSAR地形效应补偿的植被参数反演方法，其特征在于，步骤(1a)中所述的对PolInSAR包含M×N个像元的干涉相位图的平地干涉相位进行估计，是指对PolInSAR干涉相位图的每个像元的平地干涉相位进行估计，估计公式为：
其中，φflat为像元的平地干涉相位，λ为PolInSAR的工作波长，Rflat1为像元对应的平地与PolInSAR主天线的斜距，Rflat2为像元对应的平地与PolInSAR辅天线的斜距。
3.根据权利要求1所述的基于PolInSAR地形效应补偿的植被参数反演方法，其特征在于，步骤(3a)中所述的计算PolInSAR主图像上每个像元对应的主天线POA值βm，以及计算PolInSAR辅图像上每个像元对应的辅天线POA值βs，计算公式分别为：
其中，ω为像元的方位向地形坡度角，γ为像元的距离向地形坡度角，θm为PolInSAR主天线照射到该像元的视角，θs为PolInSAR雷达辅天线照射到该像元的视角。
4.根据权利要求1所述的基于PolInSAR地形效应补偿的植被参数反演方法，其特征在于，步骤(5b)中所述的通过{ωa}和{ωr}对{T11}、{T22}和{Ω12}中的每一个矩阵进行地形干涉相位坡度补偿，补偿公式为：
其中，m和n均为非负整数，且0≤m＜M，0≤n＜N，{T11}为PolInSAR主图像极化自相干矩阵的集合，{T22}为PolInSAR辅图像极化自相干矩阵的集合，{Ω12}为PolInSAR极化互相干矩阵的集合，T11(x+k,y+l)为{T11}中第x+k行y+l列的元素，T22(x+k,y+l)为{T22}中第x+k行y+l列的元素，Ω12(x+k,y+l)为{Ω12}中第x+k行y+l列的元素，{ωa}为方位向干涉相位坡度矩阵，{ωr}为距离向干涉相位坡度矩阵，ωa(x,y)为{ωa}中第x行y列的元素，ωr(x,y)为{ωr}中第x行y列的元素，且m+1≤x≤M-m，n+1≤x≤N-n。
5.根据权利要求1所述的基于PolInSAR地形效应补偿的植被参数反演方发法，其特征在于，步骤(6d)中所述的对{(γuc1,γuc2)}中的每个(γuc1,γuc2)进行地表散射判别，并对{(γcre1,γcre2)}中的每个(γcre1,γcre2)进行体散射判别，采用差异最大化方法，实现步骤为：
(6d1)计算HH+VV极化相干系数γHH+VV、HH-VV极化相干系数γHH-VV和HV极化相干系数γHV；
(6d2)判断γuc1的相位是否小于γHV的相位，若是，则将γuc1作为PolInSAR地表相干系数γg，否则，将γuc2作为PolInSAR地表相干系数γg；
(6d3)计算γcre1与γg的距离Rcre1，γcre2与γg的距离Rcre2，γHH-VV与γg的距离RHH+VV，γHH-VV与γg的距离RHH-VV；
(6d4)判断Rcre1是否大于RHH+VV，且Rcre1是否大于RHH-VV，若是，则将γcre1作为PolInSAR体散射相干系数γv，否则，则将γcre2作为PolInSAR体散射相干系数γv。
6.根据权利要求1所述的基于PolInSAR地形效应补偿的植被参数反演方法，其特征在于，步骤(7c)中所述的通过kz、θ和{(hlut,σlut)}中的每个(hlut,σlut)计算植被的体散射相干系数γlut，计算公式为：
其中，kz为有效垂直波数，θ为PolInSAR视角，α为植被场景距离向地形坡度，hlut为植被高度，σlut为消光系数。
CN201910590564.8A 2019-01-22 2019-07-02 基于PolInSAR地形效应补偿的植被参数反演方法 Pending CN110133657A (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910058439 2019-01-22
CN2019100584392 2019-01-22

Publications (1)

Publication Number Publication Date
CN110133657A true CN110133657A (zh) 2019-08-16

Family

ID=67566493

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910590564.8A Pending CN110133657A (zh) 2019-01-22 2019-07-02 基于PolInSAR地形效应补偿的植被参数反演方法

Country Status (1)

Country Link
CN (1) CN110133657A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111398957A (zh) * 2020-04-01 2020-07-10 中国林业科学研究院资源信息研究所 改进相干性计算方法的短波长双天线InSAR森林高度反演方法
CN112882026A (zh) * 2021-01-11 2021-06-01 清华大学 一种树高反演的方法、装置、计算机存储介质及终端
CN113466857A (zh) * 2021-05-11 2021-10-01 中国地质大学（武汉） 基于非局部平均的TomoSAR林下地形反演方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5552787A (en) * 1995-10-10 1996-09-03 The United States Of America As Represented By The Secretary Of The Navy Measurement of topography using polarimetric synthetic aperture radar (SAR)
JP2005127874A (ja) * 2003-10-24 2005-05-19 National Institute Of Information & Communication Technology 偏波合成開口レーダ較正方法及び装置
CN107144842A (zh) * 2017-06-27 2017-09-08 哈尔滨工业大学 一种改进的极化干涉sar植被高度联合反演方法
CN109188391A (zh) * 2018-08-29 2019-01-11 哈尔滨工业大学 一种倾斜森林地区三层S-RVoG散射模型的确定及分析方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5552787A (en) * 1995-10-10 1996-09-03 The United States Of America As Represented By The Secretary Of The Navy Measurement of topography using polarimetric synthetic aperture radar (SAR)
JP2005127874A (ja) * 2003-10-24 2005-05-19 National Institute Of Information & Communication Technology 偏波合成開口レーダ較正方法及び装置
CN107144842A (zh) * 2017-06-27 2017-09-08 哈尔滨工业大学 一种改进的极化干涉sar植被高度联合反演方法
CN109188391A (zh) * 2018-08-29 2019-01-11 哈尔滨工业大学 一种倾斜森林地区三层S-RVoG散射模型的确定及分析方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
JONG-SEN LEE ET AL.: "Polarimetric SAR Data Compensation for Terrain Azimuth Slope Variation", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
刘琦等: "极化干涉SAR森林冠层高反演的地形坡度改正", 《东北林业大学学报》 *
卢红喜: "极化干涉合成孔径雷达与层析成像技术研究", 《中国优秀博硕士学位论文全文数据库(博士) 信息科技辑》 *
卢红喜等: "双基线极化干涉SAR植被参数反演", 《西安电子科技大学学报》 *
龙江平等: "极化方位角补偿信息支持下的植被参数反演", 《测绘学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111398957A (zh) * 2020-04-01 2020-07-10 中国林业科学研究院资源信息研究所 改进相干性计算方法的短波长双天线InSAR森林高度反演方法
CN112882026A (zh) * 2021-01-11 2021-06-01 清华大学 一种树高反演的方法、装置、计算机存储介质及终端
CN112882026B (zh) * 2021-01-11 2023-04-11 清华大学 一种树高反演的方法、装置、计算机存储介质及终端
CN113466857A (zh) * 2021-05-11 2021-10-01 中国地质大学（武汉） 基于非局部平均的TomoSAR林下地形反演方法及系统
CN113466857B (zh) * 2021-05-11 2022-11-04 中国地质大学（武汉） 基于非局部平均的TomoSAR林下地形反演方法及系统

Similar Documents

Publication Publication Date Title
Schneider et al. Polarimetric and interferometric characterization of coherent scatterers in urban areas
CN110133657A (zh) 基于PolInSAR地形效应补偿的植被参数反演方法
CN106405547B (zh) 一种多极化sar简缩极化模式定标方法
JP6349937B2 (ja) 変動検出装置、変動検出方法および変動検出用プログラム
Touzi Target scattering decomposition in terms of roll-invariant target parameters
Wenxue et al. Extended three-stage polarimetric SAR interferometry algorithm by dual-polarization data
CN110109111B (zh) 极化干涉sar稀疏植被高度反演方法
CN100545676C (zh) 基于相关加权的干涉合成孔径雷达干涉相位估计方法
CN110412573A (zh) 极化sar图像分解方法及存储介质
Liao et al. The impacts of spatial baseline on forest canopy height model and digital terrain model retrieval using P-band PolInSAR data
Shen et al. A novel polarimetric PSI method using trace moment-based statistical properties and total power interferogram construction
Yu et al. Full-polarization target classification using single-polarization ground penetrating radars
CN103940834B (zh) 采用合成孔径雷达技术测量土壤湿度的方法
Touzi Target scattering decomposition of one-look and multi-look SAR data using a new coherent scattering model: The TSVM
Shi et al. An improved dual-baseline PolInSAR method for forest height inversion
CN112630741A (zh) 一种全极化合成孔径雷达图像目标补偿peoc方法
Liu et al. Building damage mapping based on Touzi decomposition using quad-polarimetric ALOS PALSAR data
Titin-Schnaider Characterization and recognition of bistatic polarimetric mechanisms
CN110058230B (zh) 基于全极化sar的陡峭地形取向角的估计方法及系统
Zhang et al. Research on inversion models for forest height estimation using polarimetric SAR interferometry
Hu et al. Model-based polarimetric target decomposition with power redistribution for urban areas
Nghia et al. Forest height estimation from PolInSAR image using adaptive decomposition method
Trisakti et al. DEM generation from stereo ALOS PRISM and its quality improvement
Lavalle et al. PolInSAR for forest biomass retrieval: PALSAR observations and model analysis
Esmaeili et al. PSInSAR Improvement using Amplitude Dispersion Index Optimization of Dual Polarimetry Data

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
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20190816

WD01 Invention patent application deemed withdrawn after publication