本发明涉及利用遥感SAR影像的林下地形反演领域，提供一种基于非局部平均的TomoSAR林下地形反演方法及系统，包括步骤：S1：获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；S2：计算所述观测值数据集中，各像素的最佳协方差矩阵；S3：通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。本发明既最大限度地利用了邻域信息，保证了联合估计出的中心像素的真实性，又通过加权的方式剔除了邻域内的噪声信息或者异质信息等干扰信息，提升了中心像素估计值的准确性，进而有效提高了TomoSAR反演林下地形的精度和普适性。
1.一种基于非局部平均的TomoSAR林下地形反演方法，其特征在于，包括步骤：S1：获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
S2：计算所述观测值数据集中，各像素的最佳协方差矩阵；
S21：选取所述观测值数据集中的某一像素作为目标像素x0，所述目标像素x0位于搜索窗口W的中心坐标(m，n)处，计算所述目标像素x0的样本协方差矩阵 和所有邻域像素xi的样本协方差矩阵 其中i表示邻域像素的编号，0＜i＜T，T为搜索窗口W内的像素总数；
S22：计算所述目标像素x0的匹配窗口P与邻域像素xi的匹配窗口的空间相似度fs(x0，xi)和辐射相似度fr(x0，xi)；
S23：通过所述空间相似度fs(x0，xi)和所述辐射相似度fr(x0，xi)，计算获得所述邻域像素xi的权重 计算公式为：S24：重复步骤S22到S23，获得所有邻域像素的权重，联合所有所述邻域像素的样本协方差矩阵及其权重进行加权计算，获得所述目标像素x0的最佳协方差矩阵，计算公式为S25：重复步骤S21到S24，获得各像素的最佳协方差矩阵；
S3：通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。
2.根据权利要求1所述的基于非局部平均的TomoSAR林下地形反演方法，其特征在于，步骤S1中，所述对各所述SLC影像数据进行预处理，具体为：将各所述SLC影像数据，进行单视复数影像序列配准、去平地效应和相位补偿操作，所述相位补偿包括：去斜、去大气扰动和去轨道误差。
3.根据权利要求1所述的基于非局部平均的TomoSAR林下地形反演方法，其特征在于，步骤S21中，像素的样本协方差矩阵R的计算公式为：其中，L表示多视数，g(l)表示l处SLC影像HH极化通道观测值，H为共轭转置操作符。
4.根据权利要求1所述的基于非局部平均的TomoSAR林下地形反演方法，其特征在于，所述搜索窗口W的大小为w×w，所述目标像素x0的匹配窗口P的大小为p×p，其中，w的取值范围为11‑21，p的取值范围为3‑5；
步骤S22中，所述空间相似度fs(x0，xi)的计算公式为：所述辐射相似度fr(x0，xi)的计算公式为：
其中，x0表示目标像素，xi表示邻域像素，P表示目标像素x0的匹配窗口P中像素的数量；
其中，γ为自定义滤波尺度因子，γs为控制着滤波的空间范围，γr为根据两个像素之间的辐射相似度来控制过滤量；
函数 为两个厄米特矩阵之间的仿射不变距离，表示为：其中，||·||F表示Frobenius范数，log表示求矩阵对数。
5.根据权利要求1所述的基于非局部平均的TomoSAR林下地形反演方法，其特征在于，步骤S3中，所述谱估计方法包括：波束形成方法、自适应波束形成法和多重信号分类法；
其中，a(zd)为映射矩阵A(z)＝[a(z1)，a(z2)，…，a(zD)]的第d个映射矢量，N为SLC影像数据的景数，为像素的最佳协方差矩阵；
6.一种基于非局部平均的TomoSAR林下地形反演系统，其特征在于，包括以下模块：观测值数据集获取模块，用于获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
最佳协方差矩阵计算模块，用于计算所述观测值数据集中，各像素的最佳协方差矩阵；
S21：选取所述观测值数据集中的某一像素作为目标像素x0，所述目标像素x0位于搜索窗口W的中心坐标(m，n)处，计算所述目标像素x0的样本协方差矩阵 和所有邻域像素xi的样本协方差矩阵 其中i表示邻域像素的编号，0＜i＜T，T为搜索窗口W内的像素总数；
S22：计算所述目标像素x0的匹配窗口P与邻域像素xi的匹配窗口的空间相似度fs(x0，xi)和辐射相似度fr(x0，xi)；
S23：通过所述空间相似度fs(x0，xi)和所述辐射相似度fr(x0，xi)，计算获得所述邻域像素xi的权重 计算公式为：S24：重复步骤S22到S23，获得所有邻域像素的权重，联合所有所述邻域像素的样本协方差矩阵及其权重进行加权计算，获得所述目标像素x0的最佳协方差矩阵，计算公式为S25：重复步骤S21到S24，获得各像素的最佳协方差矩阵；
林下地形生成模块，用于通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。
基于非局部平均的TomoSAR林下地形反演方法及系统
技术领域
[0001] 本发明涉及利用遥感SAR影像的林下地形反演领域，尤其涉及一种基于非局部平均的TomoSAR林下地形反演方法及系统。
背景技术
[0002] 林下地形作为重要的森林资源调查参数，不仅影响森林资源的空间分布，还与森林生态系统的稳定性密切相关。然而，在森林覆盖区，传统航测或光学遥感手段只能获取森林冠层顶部的高度信息，无法获得真正的林下地形。长波合成孔径雷达层析(Tomographic Synthetic Aperture Radar，TomoSAR)技术，尤其是长波长SAR系统，能够穿透森林冠层到达地面并记录森林垂直结构信息，为林下地形测绘提供了可能。SAR层析技术通过引入谱估计理论，对每个分辨单元内的散射回波在高度向上进行断层扫描，获得每个高度位置的散射回波，将传统二维成像扩展为三维成像，能够有效分离同一分辨单元内不同高度位置的散射体，已被大量用于林下地形反演。
[0003] 对于森林这类分布式散射体，其垂直向的后向散射功率包含在协方差矩阵的幅度和相位中，因而森林SAR层析三维成像均是对协方差矩阵进行处理。但是真实的协方差矩阵无法获取，通常是通过样本协方差矩阵来直接代替协方差矩阵，以此来进行SAR层析三维聚焦。目前，常用的样本协方差矩阵估计是局部平均(local means,LM)法，即将滑动窗口内的所有像素进行统计平均。该方法虽然简单易行，但是仅当窗口内的像素统计信息一致时才能获得正确的估计值。如果不加区分，会造成协方差矩阵估计不准确，易引起不同散射机制的混合叠加，层析谱不精细，损失细节信息，降低地形信息估计的精度，如坡面经过LM法处理后会变成平坦的地面。
[0004] 上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。
发明内容
[0005] 本发明的主要目的在于，解决现有技术中，协方差矩阵估计不准确，易引起不同散射机制的混合叠加，层析谱不精细，损失细节信息，降低地形信息估计的精度的技术问题。
[0006] 为实现上述目的，本发明提供一种基于非局部平均的TomoSAR林下地形反演方法，包括步骤：
[0007] S1：获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
[0008] S2：计算所述观测值数据集中，各像素的最佳协方差矩阵；
[0009] S3：通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。
[0010] 优选地，步骤S1中，所述对各所述SLC影像数据进行预处理，具体为：
[0011] 将各所述SLC影像数据，进行单视复数影像序列配准、去平地效应和相位补偿操作，所述相位补偿包括：去斜、去大气扰动和去轨道误差。
[0012] 优选地，步骤S2具体为：
[0013] S21：选取所述观测值数据集中的某一像素作为目标像素x0，所述目标像素x0位于搜索窗口W的中心坐标(m,n)处，计算所述目标像素x0的样本协方差矩阵 和所有邻域像素xi的样本协方差矩阵 其中i表示邻域像素的编号，0＜i＜T，T为搜索窗口W内的像素总数；
[0014] S22：计算所述目标像素x0的匹配窗口P与邻域像素xi的匹配窗口的空间相似度fs(x0,xi)和辐射相似度fr(x0,xi)；
[0015] S23：通过所述空间相似度fs(x0,xi)和所述辐射相似度fr(x0,xi)，计算获得所述邻域像素xi的权重 计算公式为：
[0016]
[0017] S24：重复步骤S22到S23，获得所有邻域像素的权重，联合所有所述邻域像素的样本协方差矩阵及其权重进行加权计算，获得所述目标像素x0的最佳协方差矩阵，计算公式为
[0018]
[0019] S25：重复步骤S21到S24，获得各像素的最佳协方差矩阵。
[0020] 优选地，步骤S21中，像素的样本协方差矩阵R的计算公式为：
[0021]
[0022] 其中，L表示多视数，g(l)表示l处SLC影像HH极化通道观测值，H为共轭转置操作符。
[0023] 优选地，所述搜索窗口W的大小为w×w，所述目标像素x0的匹配窗口P的大小为p×p，其中，w的取值范围为11‑21，p的取值范围为3‑5；
[0024] 步骤S22中，所述空间相似度fs(x0,xi)的计算公式为：
[0025]
[0026] 所述辐射相似度fr(x0,xi)的计算公式为：
[0027]
[0028] 其中，x0表示目标像素，xi表示邻域像素，P2表示目标像素x0的匹配窗口P中像素的数量；
[0029] 和 的定义式为；
[0030] gγ(x)＝exp[‑(γ‑1x)2]
[0031] 其中，γ为自定义滤波尺度因子，γs为控制着滤波的空间范围，γr为根据两个像素之间的辐射相似度来控制过滤量；
[0032] 函数 为两个厄米特矩阵之间的仿射不变距离，表示为：
[0033]
[0034] 其中，‖·‖F表示Frobenius范数，log表示求矩阵对数。
[0035] 优选地，步骤S3中，所述谱估计方法包括：波束形成方法、自适应波束形成法和多重信号分类法；
[0036] 所述波束形成方法的计算公式为：
[0037]
[0038] 其中，a(zd)为映射矩阵A(z)＝[a(z1),a(z2),…,a(zD)]的第d个映射矢量，N为SLC影像数据的景数，为像素的最佳协方差矩阵；
[0039] 所述自适应波束形成法的计算公式为：
[0040]
[0041] 所述多重信号分类法的计算公式为：
[0042]
[0043] 其中，U表示观测值的特征值向量，T为转置操作符。
[0044] 一种基于非局部平均的TomoSAR林下地形反演系统，包括以下模块：
[0045] 观测值数据集获取模块，用于获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
[0046] 最佳协方差矩阵计算模块，用于计算所述观测值数据集中，各像素的最佳协方差矩阵；
[0047] 林下地形生成模块，用于通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。
[0048] 本发明具有以下有益效果：本发明提出的方法既最大限度地利用了邻域信息，保证了联合估计出的中心像素的真实性，又通过加权的方式剔除了邻域内的噪声信息或者异质信息等干扰信息，提升了中心像素估计值的准确性，进而有效提高了TomoSAR反演林下地形的精度和普适性。
附图说明
[0049] 图1为本发明实施例的方法流程图；
[0050] 图2为本发明实施例的谱估计剖面线位置图；
[0051] 图3为本发明实施例的层析谱图；
[0052] 图4为本发明实施例的林下地形剖面线高度图；
[0053] 图5为本发明实施例的Lidar观测得到的数字地形模型图；
[0054] 图6为本发明实施例的林下地形图；
[0055] 图7为本发明实施例的系统结构图；
[0056] 本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0057] 应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0058] 本发明提供一种基于非局部平均的TomoSAR林下地形反演方法，包括步骤：
[0059] S1：获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
[0060] S2：计算所述观测值数据集中，各像素的最佳协方差矩阵；
[0061] S3：通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，最后寻找地面散射相位中心，获得高精度林下地形。
[0062] 本实施例中，步骤S1中所述对各所述SLC影像数据进行预处理，具体为：
[0063] 将各所述SLC影像数据，进行单视复数影像序列配准、去平地效应和相位补偿操作，所述相位补偿包括：去斜、去大气扰动和去轨道误差；
[0064] 具体实现中，单视复数影像序列配准包括粗配准和精配准，配准方法可选用基于区域的或者基于特征的，主影像的选取一般选择时间基线、空间基线居中的影像，以尽量减小时间去相关与空间去相关效应，保证干涉图的质量；去平地效应主要是去除因平坦地面引起距离向和方位向呈现周期性变化的平地相位，可选用基于观测数据的频率估计方法或基于成像几何系统参数的频率估计方法；相位补偿可选用基于永久散射体干涉(Persistent Scatterer InSAR,PSInSAR)的相位补偿或者基于小基线集干涉(Small Baseline Subsets,SBAS)的相位补偿方法。
[0065] 本实施例中，步骤S2具体为：
[0066] S21：选取所述观测值数据集中的某一像素作为目标像素x0，所述目标像素x0位于搜索窗口W的中心坐标(m,n)处，计算所述目标像素x0的样本协方差矩阵 和所有邻域像素xi的样本协方差矩阵 其中i表示邻域像素的编号，0＜i＜T，T为搜索窗口W内的像素总数；
[0067] S22：计算所述目标像素x0的匹配窗口P与邻域像素xi的匹配窗口的空间相似度fs(x0,xi)和辐射相似度fr(x0,xi)；
[0068] S23：通过所述空间相似度fs(x0,xi)和所述辐射相似度fr(x0,xi)，计算获得所述邻域像素xi的权重 计算公式为：
[0069]
[0070] S24：重复步骤S22到S23，获得所有邻域像素的权重，联合所有所述邻域像素的样本协方差矩阵及其权重进行加权计算，获得所述目标像素x0的最佳协方差矩阵，计算公式为
[0071]
[0072] S25：重复步骤S21到S24，获得各像素的最佳协方差矩阵。
[0073] 本实施例中，步骤S21中，像素的样本协方差矩阵R的计算公式为：
[0074]
[0075] 其中，L表示多视数，g(l)表示l处SLC影像HH极化通道观测值，H为共轭转置操作符。
[0076] 本实施例中，所述搜索窗口W的大小为w×w，所述目标像素x0的匹配窗口P的大小为p×p，其中，w的取值范围为11‑21，p的取值范围为3‑5；
[0077] 步骤S22中，所述空间相似度fs(x0,xi)的计算公式为：
[0078]
[0079] 所述辐射相似度fr(x0,xi)的计算公式为：
[0080]
[0081] 其中，x0表示目标像素，xi表示邻域像素，P2表示目标像素x0的匹配窗口P中像素的数量；
[0082] 和 的定义式为；
[0083] gγ(x)＝exp[‑(γ‑1x)2] (6)
[0084] 其中，γ为自定义滤波尺度因子，γs为控制着滤波的空间范围，γr为根据两个像素之间的辐射相似度来控制过滤量；
[0085] 函数 为两个厄米特矩阵之间的仿射不变距离，表示为：
[0086]
[0087] 其中，‖·‖F表示Frobenius范数，log表示求矩阵对数。
[0088] 本实施例中，步骤S3中所述谱估计方法包括：波束形成方法、自适应波束形成法和多重信号分类法；
[0089] 所述波束形成方法的计算公式为：
[0090]
[0091] 其中，a(zd)为映射矩阵A(z)＝[a(z1),a(z2),…,a(zD)]的第d个映射矢量，N为SLC影像数据的景数， 为像素的最佳协方差矩阵；
[0092] 所述自适应波束形成法的计算公式为：
[0093]
[0094] 所述多重信号分类法的计算公式为：
[0095]
[0096] 其中，U表示观测值的特征值向量，T为转置操作符。
[0097] 参考图2‑6，通过实验对本发明的基于非局部平均的TomoSAR林下地形反演方法进行验证，其中，图2为本实验选取的谱估计剖面线位置图，底图为实验区域的Pauli基影像；图3(a)、(c)、(e)分别为利用基于LM的Beamforming、Capon、MUSIC方法估计出的层析谱，图3(b)、(d)、(f)分别为利用基于NLM(Nonlocal Means)的Beamforming、Capon、MUSIC方法估计出的层析谱；图4(a)为利用基于LM的Beamforming、Capon、MUSIC方法得到的林下地形剖面线高度，图4(b)为利用基于NLM的Beamforming、Capon、MUSIC方法得到的林下地形剖面线高度；图5为Lidar观测得到的数字地形模型(Digital Terrain Model,DTM)，即为本实验验证所参考的真值；图6(a)、(c)、(e)分别为利用基于LM的Beamforming、Capon、MUSIC方法反演出的林下地形，图6(b)、(d)、(f)分别为利用基于NLM的Beamforming、Capon、MUSIC方法反演出的林下地形。
[0098] 实验区域及数据集描述：研究区域是瑞典北部Kryclan流域的寒带森林，包括针叶林，例如苏格兰松树、挪威云杉以及桦树。这里的年平均气温约为1℃，年平均降水量为600mm。而且，该区域地形地势为丘陵，高程从190m到290m不等。平均树高约为18m，最大树高为30m。表1和表2分别列出了E‑SAR机载系统参数的详细信息以及干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)对的基线信息。
[0099] 表1 E‑SAR机载系统参数
[0100]
[0101] 表2 InSAR对基线信息
[0102]
[0103]
[0104] 从图3中可以明显地看出，基于LM的谱估计结果在很多地方未能成功分离出地面和冠层散射贡献，而是将冠层的散射相位中心视为地面的散射相位中心。反观基于NLM的谱估计结果，其基本上能将地面和冠层的散射相位中心区分开。除此之外，后者估计得到的层析谱较前者更为清晰和连续，这主要得益于利用NLM估计出的更为真实准确的协方差矩阵，考虑更多的非局部邻域信息，通过寻找与自身空间相似度和辐射相似度较高的像素来加权反映自身信息，从而避免了直接草率地采用局部均值所带来的错误信息，如噪声亦或是与其特征相差较大的像素。为了定量研究新方法的有效性，本文计算了不同方法估计结果的均方根误差值(Root mean square error,RMSE)，从表3中可以直观地看出，引入NLM方法后，三种谱估计方法的精度有了较为明显的改善，分别提升了34.87％、38.28％、31.61％。
[0105] 表3剖面分析不同方法RMSE对比
[0106]
[0107] 比较图5和图6中可以得到，基于LM的谱估计方法和基于NLM的谱估计方法反演出的林下地形和Lidar的DTM测量值都较为接近，但仔细对比可以看出，基于LM方法的结果中存在较多的林下地形高估区域，而基于NLM方法的结果中较少出现这种情况。分析原因可知，主要是基于LM方法在部分区域未能区分开地面和冠层的散射相位中心，从而导致将冠层的高度误当成地面高度。这和前面的剖面层析谱分析比较中的结论是一样的，除此之外，NLM方法在地形高度变化较大的地方也表现出了不错的反演性能。因此，NLM方法在地势起伏较大的大范围区域内也是适用的，且反演精度较高。表4为两种方法反演林下地形的RMSE结果，该结果从定量的方面证明了NLM方法较LM方法在林下地形反演精度上提升了30％以上。
[0108] 表4林下地形反演不同方法RMSE对比
[0109]
[0110] 综上所述，本发明提及的NLM方法相比于传统的LM方法可得到更为精确的林下地形。
[0111] 参考图7，本发明提供一种基于非局部平均的TomoSAR林下地形反演系统，包括以下模块：
[0112] 观测值数据集获取模块10，用于获取N景SLC影像数据，对各所述SLC影像数据进行预处理，获得观测值数据集；
[0113] 最佳协方差矩阵计算模块20，用于计算所述观测值数据集中，各像素的最佳协方差矩阵；
[0114] 林下地形生成模块30，用于通过谱估计公式对各所述像素的最佳协方差矩阵进行谱估计，获得高精度林下地形。
[0115] 需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。
[0116] 上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。
[0117] 以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。
