基于SAR图像纹理的北极海冰厚度的反演研究

于淼 卢鹏 李志军 石立坚



基于SAR图像纹理的北极海冰厚度的反演研究

于淼1卢鹏1李志军1石立坚2,3

(1大连理工大学海岸和近海工程国家重点实验室, 辽宁 大连 116024;2国家卫星海洋应用中心, 北京 100081;3国家海洋局空间海洋遥感与应用研究重点实验室, 北京 100081)

基于七景北极Radarsat-2 SAR图像以及中国第六次北极科学考察走航期间利用船侧录像观测获得的平整冰厚度数据, 通过灰度共生矩阵计算纹理, 确定了最适合反演海冰厚度的纹理参数。并分析了海冰厚度与纹理之间的相关关系, 探讨了纹理反演海冰厚度的可能性。选取了最合适的纹理特征进行拟合, 并利用所得经验方程进行反演验证, 结果与实测数据吻合较好, 平均相对误差13.7%。与传统的仅依靠后向散射系数反演海冰厚度进行对比, 新方法的误差更小, 证明了纹理特征反演冰厚的优势。

北极 海冰厚度 灰度共生矩阵 纹理

0 引言

海冰的存在影响了大气与海洋之间的直接能量交换, 在全球大气和海洋之间起到重要的调节作用。近年来, 北极海冰面积快速减小, 其变化影响着北极气候乃至全球的生态环境[1]。同时, 北极具有丰富的自然资源, 以及不菲的商业利益, 北极航道的商业利用符合我国的北极战略[2], 而海冰正是航道利用中需要考虑的主要因素。

由于北极地区复杂恶劣的气候和海洋环境, 北极海冰物理参数实测数据的获取很困难, 其中, 厚度是最重要的参数之一, 也是最难获取的参数。目前海冰厚度的获取最准确的方法是钻孔测量, 同样也是最消耗资源的方法。于是各种非接触测量技术相继出现, 如仰视声纳、电磁感应、走航摄影观测、卫星遥感等[3]。近年来遥感技术的发展充分发挥了其大尺度、可重复和成本低的优势, 为长期动态获取北极海冰信息提供了新思路。其中合成孔径雷达(Synthetic Aperture Radar, SAR)由于具有对光照条件不敏感、不受云层遮挡的特点, 特别适合复杂的天气条件, 能够实现全天候、全天时的冰情监测。目前已经广泛用于海冰范围、分类、密集度、速度以及厚度等基础参数的观测研究[4-10]]。

影响海冰SAR图像散射特征主要有两个因素, 海冰表面粗糙度和复介电常数[11-12], 两者均与海冰生长过程密切相关, 可与海冰厚度联系起来用于冰厚的反演。表面粗糙程度主要是通过纹理特征以及后向散射系数(0)体现: 纹理目前在冰型分类的方面应用比较多[13], 在厚度方面应用相对较少, 而且精度不够理想[14]。０用于海冰厚度反演有一定的积累, 但是存在着低估海冰厚度的情况以及适用地域限制[15], 一般需要结合其他的产品共同使用。复介电常数的应用主要是与之相关的极化参数, Alpha角()是从全极化SAR数据中提取出来关于海冰介电常数和入射角的函数, 已经被证明可以用来反演海冰厚度[16]。

本文尝试利用纹理特征反演海冰厚度主要基于以下经验理论: 海冰越厚, 生长时间越长, 所受到的外力作用也会越多越复杂, 相应的变形及表面的粗糙程度也会发生变化[15]。本文的纹理特征分析方法采用灰度共生矩阵(GLCM), 它能够反映出图像要素的自身特征以及像素之间的空间关系, 在遥感领域纹理特征分析方面的应用最为广泛[17-19]。利用∑AP图像以及匹配的平整冰厚度数据, 分析纹理特征和冰厚之间的相关性, 选出最合适的纹理特征拟合经验方程, 然后进行冰厚的反演及误差分析, 并与传统的厚度反演方法进行比较分析。

1 研究区域与数据源

本文使用的SAR图像数据是七景北冰洋Radar­sat-2双极化(HH+HV)SGF产品, 是中国第六次北极科学考察的途经区域, 图像区域大部分被一年冰覆盖, 除8月11日图像有部分开阔水面, 其他图像几乎没有明显冰裂缝或开阔水, 适合提取冰表面纹理进行研究。图像详细参数见表1。研究区域如图1所示, 图中颜色代表2014年8月平均海冰密集度, 数据来自德国不莱梅大学发布的AMSR2/ASI产品, 黑色矩形代表SAR图像位置, 黑色虚线代表第六次北极科考船航线｡

表1 SAR图像参数

图1 研究区域位置示意图

Fig.1. Map of study sites

现场冰厚数据来自第六次北极科考, 是利用船侧观测录像提取的。与卫星等大尺度遥感方式相比, 数字摄影产品可以有针对性的对一个区域进行观测, 并且具有更高的分辨率, 能够弥补大尺度遥感的分辨率不足和观测位置的局限。首次北极科学考察就已经对密集冰区进行了航空数字相机遥感观测, 得到了一批不同类型的高分辨率海冰数字影像[20]。目前船侧录像已经广泛用在极地考察中进行冰情的观测, 它的主要优点就是自动化程度高, 可以自动记录走航路线中的冰情, 而且观测结果可以达到很高的精度[21]。其基本原理如图2所示,当船在冰区前进, 冰层被压碎后随着船体向后滑动, 由于船体的阻碍而发生侧翻露出断面, 根据断面的宽度及参照物就可以计算出相应的海冰厚度[22]。即从观测录像中筛选出带有侧翻海冰的信息, 根据图上翻冰厚度和参照球直径的比例, 乘以参照球的直径就可以得到冰层的实际厚度, 此方法在平整冰的应用上能达到很高的精确度。

图2 冰厚获取示意图

Fig.2. Method of getting ice thickness

在本文的研究中使用的海冰厚度数据是经过预处理, 提取完成的厚度数据, 包括拍摄时间, 经纬度以及对应的海冰厚度, 共计981个测量点。为保证精确度, 数据中冰厚最大和最小的两个区间部分数据过于零散, 很难得到有意义的统计数据, 予以剔除, 在余下970个测量点的基础上对平整海冰进行厚度反演研究。如图3所示为厚度数据的地理位置及厚度大小, 图中红线为科考船航线。并将数据随机分为两组, 其中一组数据用于相关分析, 确定待定系数拟合反演方程,另一组数据用于反演厚度并计算误差。

图3 厚度数据位置及大小示意图

Fig.3. Location and value of ice thickness

2 研究方法

2.1 图像预处理

SAR数据预处理主要包括辐射校正、几何校正、滤波降噪处理和入射角校正[23]。其中辐射校正是将图像的像元亮度值消除误差后转化为真实的后向散射系数, 并以dB为单位保存。再通过几何校正修正图像的几何畸变, 并使用Polar Stereo­graphic投影进行图像投影, 将图像坐标与经纬度进行对应。然后使用Refined Lee滤波处理来降低由于雷达波相干作用产生的斑点噪声。最后还用一个基于入射角和后向散射系数之间关系的经验线性方程对入射角进行校正, 降低由于入射角的变化引起的图像在距离向灰度的变化[14,24]。

2.2 GLCM参数选取

纹理是所有物体表面先天具有的特性, 反应了物体表面结构组织排列的方式以及它们与周围环境的关系。如图4所示, 灰度共生矩阵是用图像灰度级之间联合条件概率密度函数(,,,)表示纹理, 该函数表示在给定空间距离和方向时, 灰度为起点, 出现灰度级为的概率[25]。所有的()用矩阵的形式表示, 称为灰度共生矩阵, 用来量化图像灰度的空间关系。

表2是Haralick定义的部分纹理特征值, 根据矩阵便可以计算出多个纹理特征的统计值[19]。

由于C波段同极化对地物表面粗糙程度的变化比交叉极化更加敏感[26]。本文使用HH极化方式的图像进行纹理提取。在利用灰度共生矩阵提取纹理分析时, 首先要确定计算灰度共生矩阵时所需的4个参数:窗口大小､方向､位移量和灰度量化级别[25]｡在实际计算GLCM时, 改变以上4个参数的任何一个均会影响GLCM中的元素数值, 这就会直接影响纹理特征值的大小, 使得提取的纹理特征不能真实地反映物体的纹理特征, 进而影响与海冰厚度的相关性[27]｡但这4个参数往往难以确定, 目前没有足够的理论支持参数应该如何选取, 多是凭借个人经验。以往的文献表明方向和灰度量化级别对海冰的纹理特征影响很小[28], 方向一般选择为0°、45°、90°、135°四个离散方向, 研究中多采用其中某一方向或者四个方向的纹理平均值[29]。本文选取灰度量化级别为64, 四个方向取纹理平均值｡还需要确定窗口大小和位移, 这里用5×5、7×7、9×9、11×11的窗口分别以1、3、5、7的位移进行纹理计算, 以确定最优组合。

图4 GLCM的计算示意及矩阵形式

Fig.4. Diagram of GLCM and form of matrix

表2 GLCM纹理特征值表

根据厚度数据中的经纬度提取纹理特征, 以厚度数据的经纬度值为中心点, 利用3×3的窗口从SAR图像中获取该点的平均纹理特征值, 并与对应的冰厚进行匹配｡由于同一海冰厚度对应的每个纹理值会在一定范围内出现波动, 用直接得到的纹理值和海冰厚度进行分析, 很难得到两者之间的相关性。以往的研究中都采用了不同方式的分段处理, 文献[[14]]中采用K-means聚类方法对图像进行分割, 将关于像素的冰厚数据转换成关于图像片段的冰厚数据, 文献[[15]]中采用边缘分割的方法对图像进行预处理, 然后对分割片段内数据进行平均来分析相关性。以上方法比较适合大量的连续数据, 可以保证在每一片段内的数据都可以达到足够的数量, 以便于提取有意义的统计数据。由于本文使用的厚度数据分布很不规律, 直接采用分段平均处理的方法[16], 将零散的厚度数据分组, 转换成厚度段数据, 在定性分组的基础上进行定量的分析, 即在实测数据范围内, 每0.1 m的范围分为一段, 并计算每一厚度段内的平均纹理值作为这一厚度段的纹理值, 便于之后的冰厚与纹理特性之间的相关分析｡

图5是不同参数组合下纹理特征和厚度的相关系数(<0.01), 横坐标中5.1代表5×5窗口, 位移为1, 以此类推。可以看出对比度、相异性、同质性3个特征在小窗口下比较好, 在大位移下受影响很大; 二阶矩、能量、最大概率、熵对窗口大小不敏感, 大位移对其有一定影响; 均值、相关性、方差3个特征在大窗口下略好, 受位移影响很小。综上, 考虑应采用中等窗口小位移来提取纹理, 以最大程度上保证各特征值的相关性, 本文选择9×9窗口、=1来进行后续分析。

图5 不同纹理参数组合下纹理与冰厚的相关系数(绝对值)

Fig.5. Correlation coefficient (absolute value) between texture feature and ice thickness depending on the window size and the displacement

表3 相关系数汇总

2.3 相关分析

GLCM的参数确定为窗口大小9×9, 位移量为1, 灰度量化级别64, 取四个方向纹理平均值。在此基础上各个纹理特征与海冰厚度的相关系数如表3所示。同时还沿袭传统方法, 提取了图像的后向散射系数HH, 并分析其与海冰厚度的相关性, 来与纹理反演厚度进行对比分析。由于GLCM的某些纹理之间存在公式上的相似, 彼此相关, 在多纹理相关分析时应该避开, 以及考虑纹理特征的波动情况, 最终选取熵和相关性两个纹理特征用于厚度反演。纹理值和厚度之间呈现出比较明显的线性关系, 利用第一组匹配好的实测数据确定待定系数, 从而得到海冰厚度与散射特征之间的经验方程。

得到三个不同的海冰厚度的线性经验方程, 其中方程(1)是只依靠后向散射系数, 是以往研究中最常用的反演厚度方法, 方程(2)是本文使用的纹理反演方法, 方程(3)是结合后向散射系数和纹理的反演方法。

=43.10·HH+730.22 (1)

=217.58·+1307.71·r–3277.68 (2)

=33.37·HH+254.28·+46.89·–2060.68 (3)

2.4 反演结果

将相关分析得到的三个经验方程回代到第二组匹配好的厚度分层数据中来验证规律的可靠性, 并对三个反演方程进行对比, 如图6—图8所示, 其中横坐标为每一分层内的所有测量点的平均厚度, 纵坐标为相应的平均反演厚度。从图6中可以看出, 在使用传统方法HH反演厚度上出现了比较明显的厚冰厚度被低估的现象, 以往的研究中提过此问题[15], 但并未做出解释, 原因可能是海冰厚度达到一定程度后其表面的粗糙程度变化不够明显。图7、图8显示利用纹理信息对冰厚进行反演在已有数据范围内获得了很好的表现。从表4中可以看出有纹理参与的反演方式获得的反演结果误差明显优于传统方法, 且纹理和HH共同反演的结果和纹理单独反演相差很小, 可能是由于纹理中已经包含了HH中的大部分信息。

后向散射系数HH反演结果如图6所示:

图6 σHH反演海冰厚度结果

Fig.6. Comparison between measured and retrieval ice thickness depending onHHmodel

纹理(熵和相关性)反演结果如图7所示:

Fig.7. Comparison between measured and retrieval ice thickness depending on texture feature model

三特征要素(HH、熵和相关性)反演结果如图8所示:

图8 σHH和纹理反演海冰厚度结果

Fig.8. Comparison between measured and retrieval ice thickness depending onHHand texture feature model

表4 不同反演方式及反演精度

3 结论

本文利用七景SAR数据和对应的平整海冰厚度数据, 利用灰度共生矩阵, 首先研究了纹理的四个参数值对海冰厚度和纹理特征之间的相关性的影响, 确定了最适合反演海冰厚度的纹理参数和纹理特征, 并得到了利用纹理信息计算海冰厚度的经验公式。利用经验公式对SAR图像上的海冰进行了厚度的反演, 并与传统的0厚度反演方式进行了对比。主要结论如下:

1. GLCM的参数选择上, 不同的纹理特征和冰厚的相关性对窗口大小和位移量的敏感度不同, 最终确定选用9×9的窗口, 位移为1来计算矩阵纹理反演厚度最合适, 同时也表明北极海冰的纹理尺度较大, 且比较细致, 这应该是由北极地区复杂的环境导致。

2. 在选择纹理特征值时, 根据公式及意义将纹理分为三组, 可以看到在以往的研究中应用比较多的对比度和二阶矩在厚度方面的相关性并不占优势, 而且数据的方差很大, 很不稳定, 并不适合用来反演计算厚度, 本文选择了熵和相关性来反演厚度。

3. 经过对比验证, 本文的方法相比于传统的只使用后向散射系数的反演方法精度更高, 表明仅用SAR直接呈现出的后向散射系数并不能完全体现出海冰厚度的变化, 而且也证实了纹理在海冰厚度的反演方面有着很高的潜力。

另外, 由于本文所使用的厚度数据来自船侧翻冰测量, 由于此方法的局限性, 海冰厚数据多来自平整的中一年冰, 80%以上的数据处于80—130 cm, 薄冰和厚变形冰部分样本很少, 对回归方程的函数类型及其待定系数会有一定的影响, 而且不足50 cm区间冰厚的反演情况完全空白, 本方法在薄冰和厚变形冰的应用情况还要验证。此外, 本方法对海冰密集度要求较高, 海水的存在会对纹理值产生很大的影响。以往的研究也表明, 海冰SAR图像的散射特征和海冰厚度的相关性在不同的地点、不同的海冰形态上的表现并不一致[30], 目前还没有发现普遍适用的规律, 本文的方法在其他地区的合理性还有待验证。而且, 如果不对纹理信息进行分段统计处理, 直接提取图上纹理信息进行厚度反演的结果并不理想。其中很重要的一个原因是SAR图像受固有的相干作用影响[31], 信噪比不高, 单一像素所体现的信息只有部分反映了地物的真实信息, 其他部分来自无关的随机噪声[32], 使得SAR图像一直以来广泛用于定性分析, 直接定点定量地获取地物信息目前还是困难重重。因此, 如果要实现SAR图像定点定量地获取信息, 可以考虑从两个方面入手: 一个是对滤波降噪算法的改进[33], 提高图像的信噪比, 另一个是结合多种数据对地物信息进行限制以获取更加准确的信息。想要在分段反演的基础上进一步提高精确度, 逐步达到定点反演厚度还需要更深入的研究。接下来的计划中, 考虑使用如卫星测高等更加连续密集的实测厚度数据, 应用图像聚类分析对海冰进行分类, 由于纹理信息更加侧重考虑像素与其周围像素之间的关系, 分类片段和纹理共同作用下的结果值得期待。

1 STROEVE J C, SERREZE M C, HOLLAND M M, et al. The Arctic’s rapidly shrinking sea ice cover: a research synthesis[J]. Climatic Change, 2012, 110(3): 1005—1027.

2 胡鞍钢, 张新, 张巍. 开发“一带一路一道(北极航道)”建设的战略内涵与构想[J].清华大学学报(哲学社会科学版), 2017,(3): 15—22.

3 季青, 庞小平, 许苏清,等. 极地海冰厚度探测方法及其应用研究综述[J]. 极地研究, 2016, 28(4): 431—441.

4 曹梅盛. 冰冻圈遥感[M]. 北京: 科学出版社, 2006.

5 国巧真, 陈云浩, 李京, 等. 遥感技术在我国海冰研究方面的进展[J]. 海洋预报, 2006, 23(4): 95—103.

6 张晰, 张杰, 孟俊敏, 等. 基于极化散射特征的极化合成孔径雷达海冰分类方法研究: 以渤海海冰分类为例[J]. 海洋学报, 2013, 35(5): 95—101.

7 KARVONEN J. Baltic sea ice concentration estimation based on C-band HH-polarized SAR data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(6): 1874—1884.

8 LEHTIRANTA J, SIIRIÄ S, KARVONEN J. Comparing C- and L-band SAR images for sea ice motion estimation[J]. Cryosphere, 2015, 9(1): 357—366.

9 KARVONEN J, CHENG B, VIHMA T, et al. A method for sea ice thickness and concentration analysis based on SAR data and a thermodynamic model[J]. Cryosphere Discussions, 2012, 6(3): 1507—1526.

10 ZHANG X, ZHANG J, MENG J, et al. Analysis of multi-dimensional SAR for determining the thickness of thin sea ice in the Bohai Sea[J]. Journal of Oceanology and Limnology, 2013, 31(3): 681—698.

11 KIM J W, KIM D J, HWANG B J. Characterization of Arctic sea ice thickness using high-resolution spaceborne polarimetric SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(1): 13—22.

12 郭井学. 基于电磁感应理论的极地海冰厚度探测研究[D]. 吉林大学, 2007.

13 孔毅, 陈曦, 艾未华, 等. 基于GLCM和小波特征的SAR海冰分类[J]. 解放军理工大学学报自然科学版, 2015,(1): 74—79.

14 SHI L, KARVONEN J, CHENG B, et al. Sea ice thickness retrieval from SAR imagery over Bohai sea[C]. Geoscience and Remote Sensing Symposium. IEEE, 2014: 4864—4867.

15 KARVONEN J, SIMILA M, HALLIKAINEN M, et al. Estimation of equivalent deformed ice thickness from Baltic Sea ice SAR imagery[C]. IEEE International Geoscience and Remote Sensing Symposium. IEEE, 2005: 5165—5167.

16 刘眉洁, 戴永寿, 张杰, 等. 拉布拉多海一年平整冰厚度SAR反演算法[J]. 中国石油大学学报自然科学版, 2014, 38(3): 186—192.

17 郭德军, 宋蛰存. 基于灰度共生矩阵的纹理图像分类研究[J]. 林业机械与木工设备, 2005, 33(7): 21—23.

18 邹亚荣, 林明森, 马腾波, 等. 基于GLCM的SAR溢油纹理特征分析[J]. 海洋通报, 2010, 29(4): 455—458.

19 刘俊, 毕华兴, 朱沛林, 等. 基于ALOS遥感数据纹理及纹理指数的柞树蓄积量估测[J]. 农业机械学报, 2014, 45(7): 245—254.

20 赵进平, 任敬萍. 从航空数字影像提取北极海冰形态参数的方法研究[J]. 遥感学报, 2000, 21(4): 271—278.

21 卢鹏. 基于图像分析的海冰几何参数和拖曳系数参数化研究[D]. 大连理工大学, 2007.

22 卢鹏, 李志军, 董西路, 等. 基于遥感影像的北极海冰厚度和密集度分析方法[J]. 极地研究, 2004, 16(4): 317—323.

23 KARVONEN J. A sea ice concentration estimation algorithm utilizing radiometer and SAR data[J]. Cryosphere, 2014, 8(5): 2213—2241.

24 郑杨龙, 卢鹏, 李志军, 等. 基于Radarsat-2 SAR图像分类与HY-2微波辐射计反演获取北极海冰密集度的比较研究[J]. 极地研究, 2016, 28(3): 413—423.

25 薄华, 马缚龙, 焦李成. 图像纹理的灰度共生矩阵计算问题的分析[J]. 电子学报, 2006, 34(1): 155—158.

26 孙少波, 车涛, 王树果, 等. C波段SAR山区积雪面积提取研究[J]. 遥感技术与应用, 2013, 28(3): 444—452.

27 梁小祎, 张杰, 孟俊敏. 溢油SAR图像分类中的纹理特征选择[J]. 海洋科学进展, 2007, 25(3): 346—354.

28 张晰, 张杰, 纪永刚. 基于纹理特征分析的辽东湾SAR影像海冰检测[J]. 海洋科学进展, 2008, 26(3): 386—393.

29HAVERKAMP D, SOH L K, TSATSOULIS C. A comprehensive, automated approach to determining sea ice thickness from SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1995, 33(1): 46—57.

30 KARVONEN J,CHENG B, SIMILÄ M. Ice thickness charts produced by C-band SAR imagery and HIGHTSI thermodynamic ice model[C]. Geoscience and Remote Sensing Symposium, 2008.

31皮亦鸣. 合成孔径雷达成像原理[M]. 成都: 电子科技大学出版社, 2007.

32 周晓光, 匡纲要, 万建伟. 多极化SAR图像斑点抑制综述[J]. 中国图象图形学报,2008,(03): 377—385.

33 LEE J S, POTTIER E. Polarimetric radar imaging: from basics to applications[M]. Boca Raton: Francis Grou, 2009.

ARCTIC SEA ICE THICKNESS RETRIEVAL BASED ON SAR IMAGE TEXTURE FEATURE

Yu Miao1, Lu Peng1, Li Zhijun1, Shi Lijian2,3

(1State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China;2National Satellite Ocean Application Service, Beijing 100081, China;3Key Laboratory of Space Ocean Remote Sensing and Application, Beijing 100081, China)

A novel method to estimate ice thickness based on seven Arctic Radarsat-2 SAR images and data for level ice thickness collected from the 6th Chinese National Arctic Research Expedition (2014) is presented. Suitable values of texture parameters are confirmed using a gray level co-occurrence matrix (GLCM). Relationships between sea-ice thickness and texture features are analyzed, and the possibility of determining sea-ice thickness based on ice texture features is discussed. Curve-fitting equations can be determined by optimum texture features. Estimated thickness agrees well with in-situ measurements, with an average relative error of 13.7%. Compared with another commonly used method that depends on a backscattering coefficient, this new method has less error, indicating texture features can be reliably used to determine ice thickness.

Arctic, sea-ice thickness, gray level co-occurrence matrix(GLCM), texture feature

2017年8月收到来稿, 2018年1月收到修改稿

国家重点研发计划专项(2016YFC1402702)、国家自然科学基金面上项目(41676187, 41376186)、国家国际科技合作专项(2011DFA22260)资助

于淼, 男, 1993年生。硕士研究生, 主要从事海冰遥感研究。E-mail: yumiao2411@mail.dlut.edu.cn

卢鹏, E-mail:lupeng@dlut.edu.cn

10. 13679/j. jdyj. 20170035