干雪深度反演的同极化相位差模型
宋依娜1
, 肖鹏峰1,2,3
, 张学良1,卓越1,马威1
1. 南京大学地理与海洋科学学院, 江苏 南京 210023;
2. 自然资源部国土卫星遥感应用重点实验室, 江苏 南京 210023;
3. 江苏省地理信息技术重点实验室, 江苏 南京 210023
基金项目:国家自然科学基金(41671344);国家科技基础资源调查专项(2017FY100502)
关键词:克兰河 积雪深度 SAR 积雪微结构 同极化相位差
引文格式:宋依娜, 肖鹏峰, 张学良, 等. 干雪深度反演的同极化相位差模型[J]. 测绘学 ,2021,50(7):905-915. DOI:
10.11947/j.AGCS.2021.20200125
SONG Yina, XIAO Pengfeng, ZHANG Xueliang, et al. The co-polarized phase difference model for dry snow depth inversion[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(7): 905-915. DOI: 10.11947/j.AGCS.2021.20200125
阅读全文:
http://xb.sinomaps.com/article/2021/1001-1595/2021-7-905.htm
引 言
积雪是冰冻圈和水圈重要的组成部分,在全球变化研究中具有重要地位[1]。一方面,积雪作为一种特殊的地表覆盖物,通过反射大部分太阳辐射以及在消融过程中吸收大量能量,对全球气候系统产生了显著的降温作用[2-3]。同时,积雪融水作为全球水资源的重要组成部分,为全球1/6人口提供了水源[4],并进一步产生了水力发电、农业灌溉等其他用途[5]。另一方面,积雪具有时空不稳定性,过量或非适时的积雪也会引发自然灾害,给人类生命和财产安全带来巨大的威胁和损失[6]。而积雪深度作为积雪的重要因子,表征局部气候环境特征与水资源条件,是水文预测、气候模拟、雪灾监测与评价等模型的重要输入参数[7]。获取高精度雪深空间分布信息可为水资源管理、气候变化研究和防灾减灾工程等提供科学支撑。
合成孔径雷达(synthetic aperture radar, SAR)因其有较高的空间分辨率[8]且对积雪特性有着较高的敏感度等优势[9-10],在流域尺度的雪深反演中具有较强的适用性[11-12]。极化SAR通过回波获取关于积雪的散射特征[13],被广泛应用于积雪相关研究,如干湿雪的分类、雪湿度和雪密度等参数的反演[14-16],但是目前的研究对极化SAR的极化相位信息的分析较少[17]。极化SAR计算同一天线、不同极化方式的极化相干,可获取极化相干系数和极化相位信息,包括同极化相位差(co-polarized phase difference, CPD)。积雪为各向异性介质,导致HH和VV极化信号在雪层中具有不同的传播速度,因此产生CPD。文献[18]根据此原理建立积雪CPD模型,定量分析了CPD与积雪各向异性结构、雪深的关系。但是,根据CPD模型反演干雪深度需要冰粒的各向异性结构(以冰粒轴间比表示)和积雪密度两个参数,其中冰粒轴间比参数比较抽象,难以通过野外观测获得。而且,目前国际上基于CPD模型进行雪深反演的研究主要使用X波段,暂未见使用C波段极化SAR数据的研究 道。文献[19]假设研究区积雪各向异性介电常数不变,使用一个实测站点处的积雪密度0.07 g/cm3和假定轴间比1.5作为雪深反演算法的输入,利用全极化TerraSAR-X数据反演整个影像的雪深,并以该实测站点的雪深进行验证,获得较高的反演精度。但是,由于积雪的各向异性、密度空间分布不均,仅使用一个实测点的积雪密度和假定轴间比反演整个影像的雪深,将给反演结果带来较大误差。文献[20]根据实测点处CPD的正负,结合该点处实测雪深、雪密度等参数,代入CPD正演模型,求得冰粒轴间比,再借助积雪密度反演算法[21]获取全局雪深分布,但此反演方法较为复杂,且积雪密度反演算法可能会带来新的误差。
1 研究区与数据1.1 研究区概况
克兰河流域位于新疆北部,阿尔泰山南麓,在北温带大陆性气候条件影响下,所形成的积雪为典型的“干寒型”积雪[22],年平均雪深在40 cm以上。克兰河属融雪径流补给为主的河流,年平均径流量6.0×108m3[23]。积雪融水对当地农业灌溉及畜牧业发展具有重要影响,研究该流域的积雪深度具有重要意义。
研究区范围如图 1所示,位于克兰河流域上游东支小东沟源头区域。研究区西南部为戈壁,地形起伏较为平缓;东北部是阿尔泰山南麓侵蚀中山山地,海拔较高,地形起伏较大,是典型的山区环境[24]。
图 1 研究区位置
Fig. 1 Location of the study area
1.2 研究数据
表 1 所用的GF-3影像的成像参数
Tab. 1 Parameters of the used GF-3 image
成像参数 | 参数值 |
成像日期 | 2018-01-17 |
成像时间 | 北京时间12:22:04 |
产品形式 | 复数据产品(SLC) |
极化方式 | 全极化 |
轨道方向 | 降轨 |
成像波段 | C波段 |
电磁波长/cm | 5.6 |
距离向分辨率/m | 2.25 |
方位向分辨率/m | 4.87 |
近端入射角/(°) | 43.39 |
(2) 野外实测数据。CPD反演雪深的可行性分析、模型参数的确定与反演结果的精度验证均需地面实测数据的支持,因此于2018年1月17日对研究区进行星地同步积雪观测,设计了积雪样带观测和积雪属性剖面观测相结合的观测方案。此时为研究区的积雪稳定期,观测期间最高气温在0℃以下,未发生降水。积雪样带观测共计44个样点,主要测量雪深、积雪密度、雪层温度等参数。积雪剖面观测共计8个剖面点,主要测量雪深、雪层密度、雪层温度、雪层介电常数、液态水含量等参数。雪深直接由钢制直尺测量,取观测点周围3 m的范围内3次测量的平均值,测量精度为0.1 cm;积雪密度由钢制直尺、托板和自制的量雪筒测量,测量精度为1 g/cm3;自积雪表层向下每10 cm划分雪层,测得积雪参数代表该雪层积雪参数,其中雪层温度利用经过标定的针式温度计进行测量,测量精度为0.1℃,雪层介电常数、密度和液态水含量通过雪特性分析仪SnowFork测量3次取平均值得到[27]。
雪深测量结果如图 2所示,可见研究区雪深空间分布不均,西南部戈壁区域积雪较浅,东北部山地区域积雪较深。这是因为来自西南方向的暖湿气团在东移过程与阿尔泰山成近正交角,地形对气流的抬升作用显著,有利于山区降水,使得山麓的降水多于河谷平原地区[28]。为了探究模型在不同积雪条件下的适用性,将研究区分为深雪区与浅雪区两个区域,分别进行反演与分析。浅雪区的雪深观测点16个,平均雪深27.5 cm;深雪区的雪深观测点36个,平均雪深72.6 cm。
图 2 观测雪深点的空间分布
Fig. 2 Position of the measured snow depth points
深雪区观测剖面5个,平均积雪密度0.23 g/cm3,平均体积含水量0.64%;浅雪区观测剖面3个,平均积雪密度0.17 g/cm3,平均体积含水量0.27%。雪层含水量都接近于0,可认为积雪为干雪,满足CPD模型条件。分层积雪密度与液态水含量的观测结果如图 3所示。可以看出浅雪区密度与含水量分层垂直差异小,深雪区分层积雪密度存在明显的垂直差异,雪层密度呈中间大、底部和顶部小的特征。深雪区底层积雪的含水量高于表层积雪,垂直分布表现出较大的差异。图 4为研究区雪粒图像,可以看出深雪区较浅雪区雪粒垂直分异更明显,积雪重结晶变质作用更强。
图 3 积雪剖面观测结果
Fig. 3 Distribution of snow depth, snow density, and snow wetness of the measured snow samples
图 4 各雪层的雪粒观测结果
Fig. 4 Observed results of snow particle of each layer
(3) 其他数据。根据阿勒泰国家基准气候站提供的气温和雪深数据,在星地同步观测前,1月13日前后发生了一次较大的降雪,自1月13日至2月底无较大降雪天气过程,平均雪深约为14.5 cm。此外,采用国家基础地理信息中心发布的全球30 m地表覆盖数据GlobeLand 30掩膜森林等复杂散射机理的区域,使其不参与雪深反演;采用30 m分辨率的SRTM DEM数据作为数字高程数据。
2.1 积雪各向异性结构
积雪的物理特性与其微观几何结构有关。在理想状态下,干雪可以看作是冰粒和空气的二相混合物[29]。为了模拟雪的各向异性结构,假定干雪冰粒为规则的椭球体,以椭球体中心为原点,建立三维笛卡儿坐标系(x,y,z),z轴为冰粒晶体对称轴,平行于重力方向,x和y平面平行于平坦地表。ax、ay和az是椭球体在x、y、z方向上的半轴长,假设ax=ay,冰粒轴间比r=ax/az,表示冰粒的形状。r=1时,为各向同性随机结构的新降雪;r> 1时,积雪在自身重力压实作用下迅速密实化[30],逐渐转变为各向异性、圆盘状水平排列结构的雪粒[31],r越大,各向异性越强;r< 1时,积雪受温度梯度作用发生重结晶和变质作用,形成各向异性、针状垂直排列的变质雪粒[32-33],r越趋近于0,各向异性越强。由新降雪到发生重结晶与变质作用变成陈雪的过程,伴随冰粒晶体形状及轴间比r的演变如图 5所示。
图 5 积雪冰粒演化过程
Fig. 5 Evolution of snow particle shape
2.2 同极化相位差模型
如图 6所示,当电磁波以θ入射角穿透一个深度为SD的干雪层,冰粒在电磁波极化电场作用下极化,依据Laplace方程及其边界条件,计算出3个方向的退极化因子Ni,i∈{x,y,z}[38]
(1)
(2)
图 6 极化信号穿透雪层产生相位差示意图(改编自文献[37]) Fig. 6 Electromagnetic wave penetrates the snow layer and results in the phase difference (adapted from reference[37])
根据Maxwell-Garnett等效介质理论计算冰粒各向异性介电常数εeff,i,i∈{x,y,z}
(3)
式中,Ni为退极化因子;εice为冰的等效介电常数;εair为空气的等效介电常数;fvol为积雪与冰的密度比,故各向异性介电常数可由积雪密度ρsnow、冰粒轴间比r计算得到。
因为折射率n与介电常数εe存在εe=n2的关系,故HH极化和VV极化的折射率nH和nV与3个轴向上的等效介电常数εeff,i、入射角θ满足以下关系
HH和VV极化信号在雪层中的折射率不同,导致后向散射信号所传播的路径和距离不同,从而产生了相位差异。利用波长与相位的关系,以及极化信号传播路径和雪层的几何关系,得到不同雪深引起的同极化相位差
(6)
2.3 雪深反演方法
根据式(6)可计算厚度为SD的积雪导致的相位差,由此可以推导雪深反演模型为
(7)
根据nV、nH与εeff,i的关系、εeff,i与积雪密度ρsnow、Ni的关系以及Ni与冰粒轴间比r的关系可知,在确定SAR数据计算CPD之后,输入积雪密度和积雪冰粒轴间比参数即可得到雪深。
(8)
图 7 积雪各向异性介电常恒定下CPD与雪深的关系
Fig. 7 Relationship between CPD and snow depth based on constant snow anisotropic relative permittivity
利用CPD模型反演雪深,首先需要根据SAR数据计算得到CPD。CPD是同一天线、同一次成像的VV与HH极化数据极化相干得到的复相干系数的相位角。复相干系数
(9)
。
3 反演结果分析
3.1 反演精度分析
图 8 反演结果R和RMSE与滤波窗口大小的关系
Fig. 8 Relationship between R, RMSE and the window size of filter
结果表明,浅雪区的最优滤波窗口为59×59像元,相关系数R为0.83,均方根误差RMSE为2.72 cm;深雪区的最优滤波窗口为33×33像元,相关系数R为0.54,均方根误差RMSE为11.69 cm。随着滤波窗口大小的增加,浅雪区和深雪区的R均呈现出先上升再下降的趋势,RMSE呈现先大幅上涨波动、然后减小、最后稳定的趋势。在滤波窗口25×25像元以下的结果,R与RMSE出现较大波动,这是因为滤波窗口较小则不能有效抑制相干斑噪声,相位的波动使相干性降低;当滤波窗口过大,则会平滑CPD的变化,削弱图像的细节,也导致相干性降低。
不同区域最优滤波窗口大小下雪深的空间分布如图 9所示。图 9(a)显示深雪区反演得到的雪深自西向东逐渐减少,南部和北部的雪深小于中部区域。图 9(c)为浅雪区雪深反演的结果,可见除西北部反演雪深明显大于其他区域以外,其他区域反演雪深呈相间分布。掩膜森林、水体等具有复杂的散射机制及受人类活动影响严重的人造地表区域的雪深结果如图 9(b)、(d)所示。
图 9 最优滤波窗口下的雪深反演结果
Fig. 9 Spatial distribution ofretrieved snow depth with optimal window size of filter
3.2 反演不确定性分析
地形通过影响水热条件从而影响积雪空间分布,地形起伏带来雷达构像几何关系的变化从而影响雷达散射机制,对雪深的反演结果产生影响[39]。为分析地形对反演精度的影响,选取坡度作为变量与反演误差进行相关性分析。如图 10所示,雪深反演误差与坡度显著相关,随着坡度的增加,雪深的反演误差呈现增加趋势,且浅雪区的相关性强于深雪区。由浅雪区实测点坡度分布范围可知,浅雪区的地形起伏较深雪区更大,由于构建CPD模型时假定冰粒的z轴平行于重力方向,复杂的地形会导致极化信号穿透雪层的几何关系发生改变,影响模型的适用性,同时电磁波回波由于受到地形起伏干扰作用会引入过多的噪声,从而增大反演误差。
图 10 反演误差与坡度的关系
Fig. 10 Correlation between slope and inversion errors
图 11 C波段极化信号积雪穿透深度模拟
Fig. 11 Simulation the penetration depth of C band polarimetric signal
图 12 模拟CPD与不同入射角之间的关系
Fig. 12 Simulation the relationship between CPD and incidence angle
3.3 反演方法比较
反演方法 | 子区域 | 相关系数R | 均方根误差RMSE/cm |
深雪区 | 0.54 | 11.69 | |
浅雪区 | 0.83 | 2.72 | |
文献[20]的方法 | 深雪区 | 0.25 | 34.43 |
浅雪区 | 0.46 | 7.81 | |
文献[19]的方法 | 深雪区 | 0.52 | 73.66 |
浅雪区 | 0.49 | 23.92 |
图 13 已有雪深反演方法的反演结果
Fig. 13 Spatial distribution of retrieved snow depth based on the Patil’s and Majumdar’s methods
4 结论
(1) 在积雪各向异性介电常数视为恒定的理想情况下,CPD仅是雪深的函数,可根据全部实测点的雪深数据与从SAR数据获取的CPD数据,拟合得到最适合研究区积雪特性的模型经验系数,反演研究区雪深。反演精度的高低与计算CPD过程中使用的滤波器的窗口大小有关,浅雪区的最优滤波窗口为59×59像元,此时R为0.83, RMSE为2.72 cm,深雪区的最优滤波窗口为33×33像元,此时R为0.54, RMSE为11.69 cm。
(2) 雪深反演误差与坡度显著相关,随着坡度的增加,雪深的反演误差呈现增加的趋势,雪深反演不确定性受雪层变质程度、含水量及卫星入射角观测几何条件影响,反演方法对于干燥、雪层变质结晶程度低、均质的积雪,以及具有大入射角的SAR卫星有更好的适用性。
初审:张艳玲
复审:宋启凡
终审:金 君
往期推荐
资讯
高光谱观测卫星成功发射
关于开展2021年度注册测绘师资格考试工作的通知
9个省份公布2021注册测绘师 名时间
中科院院士增选初步候选人公示材料情况汇总
会议
第十届全国地图学与地理信息系统学术大会(第一轮)
第三届中国空间数据智能学术会议SpatialDI 2022 (一号通知)
关于召开“2021空间规划与实景三维高峰论坛”的通知
会议通知 | 2021中国地理信息科学理论与方法学术年会(第二号通知)
《测绘学 》
智能化测绘专刊 | 张继贤:人机协同的自然资源要素智能提取方法
《测绘学 》智能化测绘专刊 | 陈军院士:智能化测绘的基本问题与发展方向
《测绘学 》2021年第8期:智能化测绘专刊,客座主编陈军院士
《测绘学 》2021年第7期 | 何秀凤:顾及边界信号及垂直约束的GNSS水汽层析方法
《测绘通 》
利用遥感生态指数新方案评价矿区生态环境
基于星载SAR的滑坡灾害形变特征分析
图像全站仪及图像测量发展与展望
《测绘通 》2021年第7期目录
《北京测绘》
《北京测绘》2021年第6期摘要推荐(下)
《北京测绘》2021年第6期摘要推荐(下)
《北京测绘》2021年第5期摘要推荐(下)
《北京测绘》2021年第5期摘要推荐(上)
《测绘科学技术学 》
摘要 |《测绘科学技术学 》2021年第3期摘要推荐
摘要 |《测绘科学技术学 》2021年第2期摘要推荐
摘要 |《测绘科学技术学 》2021年第1期摘要推荐
《测绘科学技术学 》2021年第1期重点论文推荐
《地球信息科学学 》
《地球信息科学学 》2021年第8期佳文推介
《测绘工程》
摘要 |《测绘工程》2021年第4期摘要推荐
摘要 |《测绘工程》2021年第3期摘要推荐
摘要 |《测绘工程》2021年第2期摘要推荐
摘要 |《测绘工程》2021年第1期摘要推荐
《中国空间科学技术》
《中国空间科学技术》“空间可展开结构” 专栏征文通知
《中国空间科学技术》入驻“智绘科服”融媒体平台!
《卫星应用》
敬请期待
《Journal of Geodesy and Geoinformation Science》
精选论文 |《测绘学 (英文版)》“GNSS和LBS”专刊
精选论文 |《测绘学 (英文版)》“GNSS和LBS”专刊
Yaoming MA et al. |《测绘学 (英文版)》(JGGS)精选论文
Bing HAN et al. |《测绘学 (英文版)》(JGGS)精选论文
《Satellite Navigation》
多源导航(一) | 《Satellite Navigation》重点论文推荐
导航定位应用| Satellite Navigation文章精选
《Satellite Navigation 》| “智绘科服” 迎来新伙伴加入!
声明:本站部分文章内容及图片转载于互联 、内容不代表本站观点,如有内容涉及侵权,请您立即联系本站处理,非常感谢!