一种基于起伏地形的SAR森林二次散射有效路径计算方法

摘要

本发明公开了一种基于起伏地形的SAR森林二次散射有效路径计算方法，属于合成孔径雷达森林遥感技术领域。步骤如下：一、将地表沿方位向和地距向均匀划分为一定尺度的大尺度面元；二、对于大尺度面元Ai,j，计算Ai,j的法向量根据和电磁波出射矢量计算Ai,j满足镜面反射的入射矢量令点P(xP,yP,zP)为任一散射介电粒子的中心坐标，过P作方向为的直线与Ai,j所在平面相交于点Ri,j；判断Ri,j若在Ai,j内部，则P与Ai,j存在二次散射有效路径，若存在，则经P作方向为的直线与Ai,j所在的平面相交于点Ei,j；以步骤二遍历或者以一定的顺序遍历步骤一拟合得到的大尺度面元。本方法适用于地形起伏情况下的SAR森林场景。

说明书

技术领域

本发明涉及属于合成孔径雷达森林遥感技术领域。

背景技术

近年来，合成孔径雷达SAR在林业遥感中得到了广泛的应用，成为定量提 取生物量、平均树高、垂直结构等森林参数的有效方式。

对于三维森林场景，对其地表进行模拟或者依据实测数据可以获得地面数 字高程DEM图，由合成孔径雷达SAR照射三维森林场景时，由于树干表面通 常情况下较为粗糙，因此呈现出较强的漫散射特征，镜面方向反射不明显；而 地表在远场效应、大尺度条件下呈现较强的镜面反射特性，非镜面反射方向强 度明显较弱。上述树干双基地散射、地表镜面反射的散射路径为二次散射的主 要成分，此即为SAR森林二次散射有效路径。

由此可以看出，地形起伏是影响森林参数定量反演精度的重要因素之一。 地形起伏使雷达局部入射角发生改变，造成森林雷达散射机理发生变化，大大 增加了森林参数反演的难度。只有深入了解地形对森林雷达散射的影响，并能 够在地形起伏情况下对SAR森林二次散射进行准确计算，才能够正确地反演森 林结构参数。

对于山区及丘陵地带，真实的地表大都存在坡度的变化，此时，随着地表 坡度的变化，SAR森林二次散射有效路径可能不存在也可能存在多条。而目前 已有的SAR森林二次散射计算方法均基于水平地表或小倾斜的粗糙度表，仅针 对地表的斜率固定、坡度通常不发生变化的情况适用，在这种情况下，SAR森 林二次散射有效路径最多仅有一条。因此采用目前已有的计算方法，对于山区 或丘陵地带进行二次散射有效路径的计算时也仅仅考虑最多仅有一条有效路径 的情况则会引入较大的误差，导致计算结果与真实自然环境的特点不一致，严 重影响了计算结果在SAR森林遥感信号模拟和结构参数反演方面的应用。

发明内容

有鉴于此，本发明提供了一种基于起伏地形的SAR森林二次散射有效路径 计算方法，打破了现有方法仅考虑最多仅有一条有效路径的情况仅针对水平地 表或小倾斜的粗糙度表的局限性，降低了对山区或丘陵地带进行二次散射有效 路径的计算误差。

为达到上述目的，本发明的技术方案为：一种基于起伏地形的SAR森林二 次散射有效路径计算方法，包括如下步骤：

步骤一、针对任一三维场景，获得其地表的地面数字高程DEM图，所述 DEM图由一定数量的小面元组成，对于每个小面元进行均匀采样获取采样点， 选取一定范围块对所有采样点分别沿方位向和地距向进行分块，将分于同一块 的采样点进行平面拟合，得到的拟合平面记为大面元，则地表被划分出大面元 集合；

选取范围块的依据是：对范围块内的采样点进行拟合得到的大面元为矩形 并且边长大于入射波波长；

对于拟合得到的大面元A_i,j，其方位向坡度角为α_{a_i,j}和地距向坡度角为 α_{gr_i,j}，其中i为方位向大面元序号，j为地距向大面元序号；

步骤二、由合成孔径雷达SAR向所述三维场景发射电磁波，令为电磁波 的入射矢量，为示电磁波的出射矢量，与平行，选择三维场景中一散射 点作为散射介电粒子P，本步骤分为以下步骤：

步骤201、根据α_{a_i,j}和α_{gr_i,j}计算A_i,j的法向量

步骤202、根据和计算A_i,j满足镜面反射的入射矢量

步骤203、经过点P作方向为的直线与面元A_i,j所在的平面相交于点R_i,j， 则R_i,j即为满足P与面元A_i,j有效二次散射的镜面反射点；

步骤204、判断若R_i,j在A_i,j内部，则P与A_i,_j存在二次散射有效路径，否则 P与A_i,j不存在二次散射有效路径；

步骤205、若P与面元A_i,j存在二次散射有效路径，则经过点P作方向为的直线与A_i,j所在的平面相交于点E_i,j，E_i,j即为P与A_i,j有效二次散射的等效散 射相位中心；

步骤206、将步骤一得到的大面元均针对P进行步骤201～步骤205的处理， 即可得到P的所有二次散射有效路径及相应的等效相位中心;

步骤三、三维场景中所有散射点均进行步骤二的处理，即可得到三维场景 中任一散射点的所有二次散射有效路径及相应的等效相位中心。

进一步地，本方案中步骤206为：

步骤2061、在步骤一划分的大面元中选择散射介电粒子在地表投影点所对 应的地距向面元；

步骤2062、在步骤2061获得的地距向面元中选择比散射介电粒子在地表投 影点更接近雷达一侧的面元；

步骤2063、在步骤2062获得的面元中选择散射介电粒子在地表投影点设定 范围内的面元；

步骤2064、将步骤2063获得的面元进行步骤201～步骤205的处理。

进一步地，本方案中步骤三为：

步骤31、将三维森林场景均匀划分为三维网格，所有的散射介电粒子根据 空间位置分布在对应的网格内，选取网格的最大边长小于SAR二维分辨率；

步骤32、按照步骤二计算每个网格中心处散射介电粒子的二次散射有效路 径，建立网格中心与二次散射有效路径的关系；

步骤33、针对任一散射介电粒子Q，获取Q所在网格中心的二次散射有效 路径作为Q的二次散射有效路径。

有益效果：

1、本发明深入了解地形对森林雷达散射的影响，充分考虑了地表具有大尺 度起伏特征时森林雷达二次散射可能不存在或存在多条二次散射有效路径的情 况，能够有效地表征地形起伏情况下森林的雷达二次散射特征，从而打破了现 有方法仅针对水平地表或小倾斜的粗糙度表的局限性，使得可为山区、丘陵地 带森林的SAR遥感数据模拟以及森林参数反演算法研究提供支撑。

2、本发明考虑到算法中计算量较大的问题，结合实际情况，提出了有效的 快速算法，在保证算法精度的基础上大大减少了计算量。

附图说明

图1为本发明实施例中机载激光雷达实测的徂徕山DEM图；

图2为本发明实施例中步骤1）拟合得到的大尺度面元DEM图；

图3为本发明实施例中二次散射有效路径计算示意图；

图4为本发明实施例中二次散射有效路径快速计算示意；

图5为本发明实施例中计算结果示意图，（a）为阔叶树的几何示意图，（b） 为采用基于单一坡度地表模型的计算结果，（c）为采用本发明的计算方法的计 算结果。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本实施例采用机载激光雷达获取的山东省徂徕山地区的地面数字高程DEM 图作为起伏地形实例，如图1所示。

一种基于起伏地形的SAR森林二次散射有效路径计算方法，步骤如下：

步骤一、将地表地形沿方位向和地距向均匀划分为大面元，大面元的要求 为：大面元的边长大于入射波波长。

其中所获得的地表的地面数字高程DEM图由一定数量的小面元组成，其每 个面元可以看成具有一定坡度（二维）的无限大介电平面。

由于真实地表通常叠加了小尺度的粗糙起伏，为了避免小尺度粗糙的影响， 所以本发明使用平面拟合的方式将地表沿方位向和地距向均匀划分出大面元， 具体为：在小面元的DEM图中，对于每个小面元进行均匀采样获取采样点，选 取一定范围块，对所有采样点沿方位向和地距向以该范围块进行分块，将划分 于同一块的采样点进行拟合，得到的拟合平面记为大面元，则地表被划分出大 面元集合。

选取范围块的依据是：对范围块内的采样点进行拟合得到的大面元为矩形 并且边长大于入射波波长。

对于拟合得到的大面元A_i,j，其方位向坡度角为α_{a_i,j}和地距向坡度角为 α_{gr_i,j}，其中i为方位向大面元序号，j为地距向大面元序号。

其中方位向为雷达运动方向，地距向为电磁波传播方向的水平分量。

本步骤中拟合得到的大面元DEM如图2所示。

步骤二、使用合成孔径雷达SAR发射电磁波对森林场景进行散射模拟，在 步骤1拟合得到的大平面下，二次散射有效路径由入射角度、大尺度面元的坡 度坡向和中心位置，以及散射介电粒子空间位置的共同决定。

和分别表示SAR的电磁波入射矢量和出射矢量，P点为森林场景中任 一散射介电粒子，则求解P与A_i,j之间二次散射有效路径的具体步骤为：

步骤201、根据α_{a_i,j}和α_{gr_i,j}计算A_i,j的法向量假设建立如图3所示的 空间坐标系，以竖直向上的方向为z轴，以水平面为xoy面，其中方位向为y 轴，地距向为x轴，则对于A_i,j来说，α_{a_i,j}为A_i,j面与y轴夹角，α_{gr_i,j}为A_i,j面 与x轴夹角，则根据A_i,j边长与角度的关系，可建立A_i,j边长的矢量，两相邻边 长矢量作叉乘即得A_i,j的法向量

步骤202、由法向量和电磁波出射矢量计算A_i,j满足镜面反射的入射矢 量和为镜面反射，关于对称且在和组成的面上。

步骤203、由于地表在远场效应、大尺度条件下呈现较强的镜面反射特性， 非镜面反射方向强度明显较弱；因此对于点P，其与A_i,j的有效二次散射的反射 点即为镜面反射点，求该镜面反射点的具体方法为：

经过点P作方向为的直线与面元A_i,j所在的平面相交于点R_i,j，则R_i,j即 为满足散射介电粒子P与面元A_i,j有效二次散射的镜面反射点；

步骤204、由于散射介电粒子P与面元A_i,j可能并不存在二次散射有效路径， 因此需要作出如下判断：判断R_i,j是否在面元A_i,j内部，若在内部，则散射介电 粒子P与面元A_i,j存在二次散射有效路径，否则不存在；

步骤205、若存在二次散射有效路径，则经过点P作方向为的直线与面 元A_i,j所在的平面相交于点E_i,j，则E_i,j即为P与A_i,_j有效二次散射的等效散射相 位中心。

步骤206、针对步骤1拟合得到的大尺度平面的每个面元进行上述步骤的处 理，即可得到P的在起伏地表下的所有二次散射有效路径及相应的等效相位中 心。

步骤三、所述三维场景中所有散射点均进行步骤二的处理，即可得到三维 场景中任一散射点的所有二次散射有效路径及相应的等效相位中心。

针对三维森林场景中所有的散射介电粒子进行上述步骤的处理和判断，则 可以得到整个森林场景的二次散射有效路径和对应的等效散射中心；但是，当 场景较大时，对所有面元进行遍历判断的计算量较大，因此对于步骤206本发 明采用如下方法限定所处理面元的范围：

步骤2061、由于散射介电粒子与地表发生镜面反射主要集中在沿着电磁波 入射的方向上，即地距向，因此仅选择散射介电粒子在地表投影点所对应的地 距向面元；在本实施例中，假设点P在地表投影点位于大面元A_2,3，则选取A_2,3所对应的地距向面元，即A_2,1、A_2,2、A_2,3。

步骤2062、由于散射介电粒子不可能远距方向的地表面元发生二次散射， 因此在步骤2061选择的地距向面元中选择比散射介电粒子在地表投影点更接近 雷达一侧的面元；在本实施例中，由SAR的电磁波入射矢量可以看出，SAR 应位于点P的左侧，则在步骤2061所选择的A_2,1、A_2,2、A_2,3中A_2,1、A_2,2为 更接近雷达一侧的面元。

步骤2063、进行小范围搜索，由于树冠对电磁波产生一定衰减作用，因此 即使距散射介电粒子较远处的面元存在二次散射有效路径，二次散射波返回雷 达的能量也相对较小，因此在步骤2062选择的更接近雷达一侧的面元中选择散 射介电粒子在地表投影点小范围内的面元；在本实施例中，可对该小范围进行 设定，并进行实验，选择实验结果最佳的设定值作为该小范围的设定值。在本 实施例中假设该小范围的设定值恰好仅包含与其最接近的A_2,2面，则在本次二 次散射有效路径的计算中，仅计算P与面元A_2,2的二次散射有效路径即可，由 此大大减少了计算量。

在进行二次散射有效路径的计算时，选择森林场景中的散射点作为散射介 电粒子，则三维森林场景中每株树木均由成千上万个的散射介电粒子组成，整 个森林场景的散射粒子数量可能达到百万甚至千万的量级。因此在步骤三中， 对于数量巨大的散射介电粒子，每个散射介电粒子均使用本实施例中提供的步 骤二进行计算，为减少计算量，本实施例提供如下方法，采用划分三维网格的 方式，实现二次散射有效路径的快速计算，具体为：

步骤31、将三维森林场景均匀划分为三维正方体网格，所有的散射介电粒 子根据空间位置分布在对应的网格内，选取网格的尺寸小于SAR二维分辨率， 如图4所示。

步骤32、计算每个网格中心处散射介电粒子的二次散射有效路径，建立网 格中心与二次散射有效路径的关系。

步骤33、针对任一散射介电粒子Q，获取Q所在网格中心的二次散射有效 路径作为Q的二次散射有效路径。

针对三维森林场景中所有散射介电粒子均以其所在网格的网格中心的二次 散射有效路径作为其二次散射有效路径，由此获取整个三维森林场景的二次散 射有效路径。

本发明针对在步骤1拟合得到的徂徕山区大尺度面元DEM图上生长一株阔 叶树，树高20m，阔叶树被离散化为1000个散射介电粒子，分别采用基于单一 坡度地表模型的计算方法和本发明的计算方法计算阔叶树的有效二次散射情 况，计算结果如图5（a）（b）（c）所示，其中，基于单一坡度地表模型的计算 方法计算得到的二次散射有效路径数量为1000条，即每个散射介电粒子存在且 仅存在1条，而且其等效散射相位中心分布在单一坡度的地表上；采用本发明 的方法，二次散射有效路径共计算得到2171条，平均每个散射介电粒子存在2 条，而且其等效散射中心分布范围较大，与地表面的最大高度差接近5m，与真 实地形特性相符，真实的山地、丘陵等地形中大尺度起伏特征是普遍存在的， 采用本发明的方法计算这些地区森林的SAR二次散射有效路径更加逼近实际情 况，这说明，按照本发明提供的技术方案进行起伏地形下的SAR森林二次散射 有效路径计算可以达到预期目的，可为山区、丘陵地带森林的SAR遥感数据模 拟以及森林参数反演算法研究提供支撑。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保 护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等， 均应包含在本发明的保护范围之内。