吴志鹏, 张 平, 李 震, 黄 磊, 刘 畅, 高 硕
(1. 中国科学院空天信息创新研究院, 北京 100094; 2. 三亚中科遥感研究所海南省地球观测重点实验室, 海南 三亚 572029; 3. 中国科学院大学, 北京 100049)
植被高度是森林资源质量评价体系的重要指标参数,表征了森林的生长状况及其健康程度,为科学合理地制定森林资源经营管理措施提供了有效依据。作为主要的森林结构参数之一,植被高度在林业遥感领域成为森林蓄积量、地上生物量(aboveground biomass, AGB)和碳储量估测研究的主要因子[1]。传统的基于地面实测数据的森林资源调查方法对植被高度数据的获取需要耗费大量人力物力,时效性差,无法实现森林高度的动态监测,且难以获取大范围森林高度观测信息,不能满足森林资源管理、监督和可持续利用的需求。卫星遥感技术能够快速、准确、宏观地以不同时空尺度获取地面森林资源分布、结构、动态变化及过程信息,已经广泛地应用于大尺度范围的植被高度信息获取[2-4]。而对于较小尺度森林生态系统的监测,卫星遥感的分辨率和精度往往很难满足应用需求。
近年来,轻小型无人机与遥感传感器的快速发展,为森林资源的遥感监测提供了新的手段[5]。与卫星遥感相比,无人机平台的飞行高度通常不超过500 m,属于近地观测范畴,能够获得更高分辨率和精度的目标区域观测数据[6];与航空遥感相比,无人机能够由内置的飞行控制系统实现航线规划和自动执飞功能,具有低成本、易操作、便于携带等特点,并且不受起飞场地的限制,机动性强。其搭载的不同传感器能够获得多层面的森林资源调查数据,丰富了数据的多样性,大大拓展了森林资源调查的深度和广度[7-8]。
目前,国内外无人机遥感的植被高度反演研究大都基于激光雷达(light detection and ranging, LIDAR)[9-11]或合成孔径雷达(synthetic aperture radar, SAR)[12]载荷。LIDAR利用激光扫描测量距离和角度[13],惯性测量单元(inertial measurement unit, IMU)测量姿态以及通过差分全球定位系统(differential global positioning system, DGPS)达到厘米级定位,能够获取高分辨率、高精度的数字地表模型进而提取植被高度等信息[14]。已有学者尝试使用激光扫描仪集成的无人机LIDAR系统进行植被高度反演研究,例如Zhang等[15]基于RIEGL VUX-1无人机载LIDAR系统获取了呼伦贝尔地区0.1 m分辨率的草地植被高度分布信息,并进一步估测了该区域的AGB。Guo等[16]利用配备了VLP-16激光扫描仪的无人机LIDAR系统,绘制了中国地区3种生态系统的高分辨率冠层高度模型。Sankey等[17]将无人机LIDAR数据与高光谱数据融合,提升了高生物量地区植被冠层高度和数字高程模型(digital elevation model, DEM)反演的准确率。但高昂的数据获取成本限制了LIDAR的大规模应用。SAR作为一种主动式的微波对地观测系统,其发射的信号波长较长,能够穿透部分树枝、树叶到达冠层下部,具有获得指示植被垂直方向参数的能力,已经被广泛地应用于森林高度信息提取[18]。Hensley等[19]首次将多基线极化干涉无人机SAR数据应用到森林结构反演中,先利用Rvog模型分离地体散射,然后对来自地表和冠层信号采用频谱分析的方法分别进行层析成像,进而提取植被高度。El Moussawi等[20]和Asopa等[21]利用Capon谱估计的方法对L波段无人机SAR数据进行层析成像,得到了热带雨林的三维结构,并进一步估算出了森林冠层高度。可以看出,目前主要通过SAR层析数据得到植被高度维剖面进而反演植被高度[22-23],而层析处理算法较为复杂,且要求高质量的多航过数据,因此雷达系统较为复杂,通常使用大型无人机,成本较高。本文研制了一种高度集成、轻量化、高可靠性的无人机雷达系统,能够满足轻小型无人机载荷的需要,并提出一种利用该系统获取的数据反演植被高度信息的方法,该方法模型简单、处理便捷,且只需一次航过便可反演得到植被高度值。
无人机雷达系统采用模块化设计,具备高度集成、轻量化、高可靠性等特点,平台的搭建以模块化的信号收发器为基础,利用单板计算机(single board computer,SBC)完成雷达信号发射/接收通信管理,利用超宽带信号处理器产生脉冲信号,通过对数周期天线发射、接收信号,利用数据记录板采集数据信息,整个系统集成后挂载在无人机平台下。无人机雷达协同的硬件组成如图1所示。
图1 无人机雷达协同的硬件组成
雷达系统的主控模块选用通用单板计算机,分别通过以太网、通用串行总线(universal serial bus, USB)接口完成与超宽带信号处理器、存储模块的通信连接。雷达模块选用PulsOn 440超宽带信号处理器(简称P440)产生脉冲信号,其工作频率在3.1~4.8 GHz之间,中心频率为3.95 GHz。P440能够对每一个发射脉冲的能量进行相干叠加,在不改变发射强度的基础上,累积更多的脉冲能量以增加接收信号的信噪比,有效改善了系统的接收性能。雷达系统的天线采用发射、接收双天线设计,选择对数周期天线,天线波束宽度为30°,增益设计为6 dB,频率范围可达2~11 GHz,以下视配置安装在雷达上。
为了对无人机雷达系统的有效性进行评估,于2021年4月21日在海南省文昌市八门湾红树林保护区内进行测试实验,实验区域如图2所示。共飞行2个架次,一号航线飞行距离约315 m,飞行时间约5 min,二号航线飞行距离约584 m,飞行时间约10 min。无人机飞行高度均设置为44 m,飞行速度为1 m/s,航向采样率约为5 Hz。雷达观测的快时间窗设置为166.1~332.6 ns,在600 s的慢时间范围内采集雷达原始数据,每条扫描数据在高度向上有3 264个采样点。
图2 实验区域
轻小型无人机在飞行过程中容易受到不平稳气流的影响,倘若飞行平台在气流影响下造成严重颠簸,会使获得的雷达图像失真,从而给目标检测和定位带来困难。飞控系统中记录的两次航过的飞行姿态数据显示,无人机平台飞行高度在43.9~44.0 m的范围内波动,俯仰角和偏航角的摆动幅度均不超过5°,在高度测量误差允许的范围内,可以近似认为航线姿态是稳定的。
被观测目标在无人机雷达图像上的特征取决于目标的位置、形状、散射特性、电磁波传播速度以及飞行平台的速度等[24]。由于实际应用场景中雷达波的介质环境非常复杂,接收信号的成分既包含目标回波,也包含噪声和干扰等,通常信噪比较低。为了提高无人机雷达图像质量,实现对目标信息的精确提取,需要对获得的原始数据进行处理。首先,通过回波信号变换将原始实数数据转换为雷达复数数据;然后,再降低直耦波和噪声对目标后向散射的影响;最后,利用基于互相关信息的后向投影方法进一步增强目标信号,以提取到更可靠的植被高度。数据处理方法的具体流程如图3所示。
图3 数据处理方法流程图
雷达发射高斯信号,信号形式如下:
(1)
式中:A表示信号幅度;t表示时间变量;tR表示目标位置对应的采样时刻;f0表示信号的中心频率。根据信号形式进行模拟分析,模拟结果如图4所示。
图4 雷达信号模拟分析
可以看出,将信号直接进行希尔伯特变换,即可得到分辨率较好的尖脉冲,同时获得雷达复数数据。图5(a)和图5(b)分别为1号航线和2号航线获取的原始实数数据经回波信号变换后得到的雷达复数图像数据。其中,横坐标为水平测量点的位置,每个测量点间隔0.236 m;纵坐标为高度方向采样点位置,零点表示信号截断位置,每个采样点间隔为0.01 m。
图5 雷达复数图像数据
图6 二维滤波法抑制直耦波
经直耦波抑制之后的雷达信号仍然掺杂着噪声的干扰,这些噪声会对雷达图像的分辨率造成影响,掩盖图像的细节信息,从而降低植被高度反演精度。主元分析法是一种线性的子空间投影算法,该方法基于最小均方误差准则,能够有效地压制噪声能量,提高信噪比[26-27]。目前,主元的选择方法主要采用经验值法[28]和特征能量百分比法[29],这两种方法受主观影响较大,稳定性差,难以实现主元的自主选择。本文针对于上述应用场景,提出一种联合剩余图像熵与奇异值分解的自适应主元分析法抑制噪声,以去除前k个主元之后剩余图像熵值达到局部最小值为准则,实现主元的自主选择,以达到最佳的去噪效果。
假设SR∈CM×N代表去除直耦波之后雷达数据,SR的奇异值分解表示为
(2)
(3)
(4)
ω(k)=Ω(W(k)),k=1,2,…,Q
(5)
以1号航线获得的无人机雷达数据为例,通过式(5)得到熵值ω(k)随k值变化,如图7所示。可以看出,一开始图像能量主要包含目标能量和噪声能量,且大部分能量主要集中在目标能量部分,整幅图像的能量较分散,熵值较大,随着k值增大,目标部分的能量被不断去除,整幅图像的能量趋于集中,熵值逐渐减小。当k=94 时,ω(k)达到局部最小值,此时可以认为目标能量全部被去除,整幅图像的能量都集中在噪声部分,即第94个主元之后所对应的信号全为噪声信号,那么选择第1~93个主元重构该图像,能够保留目标信息且达到较好的去噪效果。
图7 熵值随k值变化
图8为1号航线无人机雷达数据去噪结果,对比图6(b)可以看出,大量的背景噪声得到了有效去除,同时保留了目标的边缘和细节特征,目标信号更加清晰。该方法解决了传统主元分析法在噪声抑制中主元依赖人工选择、阈值不稳定的问题,避免了针对不同数据需人机交互判定主元数目带来的麻烦,具有一定的鲁棒性,能够适应不同的应用场景。
图8 基于剩余图像熵的自适应主元分析法去噪结果
后向投影算法是一种时域重建目标影像的方法,其核心是补偿每一成像点相对于每一时刻雷达天线位置的相位,并将该点对应的合成孔径内所有回波数据相干叠加,实现雷达方位向的高分辨率成像。由于本系统雷达天线具有宽视角的特点,所以其在每个时延位置处收集到的散射响应包含了天线照射范围内所有点在该时延位置处的响应,容易出现杂波干扰能量较大的问题。根据各接收通道之间目标信号和杂波相关性不同这一特性,利用后向投影算法处理时对各接收通道之间回波信号进行互相关运算可以消除杂波的干扰,增强目标响应[32]。本文利用基于互相关信息的后向投影算法对噪声抑制后的数据进行处理,如图9所示,与图8相比,杂波的影响得到了有效地抑制,目标信号与背景的对比度获得了较大提升,提高了目标信息的识别能力。
图9 不同航线基于互相关信息的后向投影算法处理结果图
当电磁波照射到随机分布,截面较大的体散射单元上时,虽然每个散射单元都有各自方向的散射响应,但总体来看其取向是随机的,后向散射回波在各个方向均匀分布。但当电磁波照射到光滑表面上时则会发生镜面反射,散射体的回波信号只在一个方向叠加,其他方向相互抵消,雷达几乎无法接收到后向散射回波。所以,雷达入射到较为平静的水体发生的主要为镜面反射,水体雷达后向散射系数几乎为0,入射到水面上的电磁波后向散射回波也近乎为0,但是水面上的植被冠层将发生体散射,加之水面与树干之间的二次散射,将导致植被区域的散射较水体散射强很多,如图10所示。因此,根据回波信号在高度向上迅速衰减这一特性,可以识别出林下地表所对应的高程值,从而直接从无人机雷达图像中反演出植被高度,避免了应用物理模型的先验信息。
图10 植被散射机制示意图
根据上述分析,雷达接收到的回波信号中,来自植被体散射的信号能量较强。基于回波能量在高度向上的分布形状,可以估测该区域植被的冠层高度。将雷达图像上回波能量分布分为相位中心区、功率损失区、噪声区3个区域[33],来自植被的回波能量绝大多数集中在相位中心区,损失区的能量由于植被密度的降低而逐渐减弱,继续沿高度向向上的能量基本来自噪声。所以,可以通过识别损失区和噪声区的分界位置提取冠层高度。从图9可以看出,经过基于通道互相关的后向投影算法处理之后,噪声区的能量已得到较好的压制。在雷达图像上应用Sobel边缘检测算子提取回波信号轮廓线所对应的高度值,即可获取到被测区域的植被冠层高度,如图11所示,图中也可以发现回波信号在高度向上约1.5 m处迅速衰减,表明了水体位置。
图11 不同航线雷达数据植被冠层高度反演结果
开展地面植被高度测量时,由于红树林大都生长在潮间带的浅滩和淤泥等测量人员难以进入的区域中,获得大面积的地面实测数据较为困难,所以在进行雷达测量的同时,仅对部分树高采用地面激光测高仪采集了树高实测值。另外一部分植被根系位于海面以下的区域采用数字地表模型作为植被高度验证数据,对基于无人机雷达数据提取到的树高进行验证,并以相关系数R、均方根误差(root mean squared error,RMSE)作为定量评价指标,具体公式如下:
(6)
(7)
图12 无人机雷达数据反演植被高度的精度
考虑到现有的植被高度遥感反演方法存在模型复杂、数据获取难度大等问题,本文设计了一种轻小型无人机雷达系统,并提出利用该系统反演植被高度信息的方法,利用飞行实验对该方法的有效性进行了评估。无人机雷达数据获取系统采用了模块化设计,具备高度集成、轻量化、高可靠性等特点;在雷达数据处理方面,采用二维滤波和主元分析法对杂波及噪声进行抑制,提出基于剩余图像熵的主元自动选择方法,以解决传统经验法和能量百分比法主元选择不稳定的问题,并利用基于互相关信息的后向投影算法进一步减弱杂波的影响,得到了更高精度的植被高度反演结果。验证结果表明,该方法提取得到的植被高度能够具有较高的精度。系统与其他方法相比,数据表现形式更加直观、模型更加简单,只需一次航过便可反演得到植被高度值,降低了数据获取成本,为反演高时空分辨率、高精度的植被高度信息提供了新的思路。