留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于正则化的地表反射特性参数反演方法

崔生成 朱文越 杨世植 李学彬

崔生成, 朱文越, 杨世植, 李学彬. 基于正则化的地表反射特性参数反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
引用本文: 崔生成, 朱文越, 杨世植, 李学彬. 基于正则化的地表反射特性参数反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
CUI Shengcheng, ZHU Wenyue, YANG Shizhi, LI Xuebin. Regularization-Based Retrieval Method for Surface Reflective Property Parameters[J]. Geomatics and Information Science of Wuhan University, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
Citation: CUI Shengcheng, ZHU Wenyue, YANG Shizhi, LI Xuebin. Regularization-Based Retrieval Method for Surface Reflective Property Parameters[J]. Geomatics and Information Science of Wuhan University, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244

基于正则化的地表反射特性参数反演方法

doi: 10.13203/j.whugis20160244
基金项目: 

国家自然科学基金 41305019

安徽省自然科学基金 1308085QD70

详细信息
    作者简介:

    崔生成, 博士, 副研究员, 主要从事地表反射特性定量遥感理论与反演方法研究。csc@aiofm.ac.cn

    通讯作者: 李学彬, 博士, 副研究员。xbli@aiofm.ac.cn
  • 中图分类号: P237;P407.9

Regularization-Based Retrieval Method for Surface Reflective Property Parameters

Funds: 

The National Natural Science Foundation of China 41305019

the Anhui Provincial Natural Science Foundation 1308085QD70

More Information
    Author Bio:

    CUI Shengcheng, PhD, associate professor, specializes in the theories and methods of surface reflective property retrieval and applications. E-mail: csc@aiofm.ac.cn

    Corresponding author: LI Xuebin, PhD, associate professor. E-mail: xbli@aiofm.ac.cn
  • 摘要: 区域和全球尺度下研究气候变化,需要准确描述大气边界层处下垫面地表的光反射特性。为此,提出基于正则化的约束反演方法,其核心在于正则化参数的选取。通过L曲线拐点位置的确定,求出最合适的正则化约束参数,以稳定陆表反射特性参数的定量反演。京津唐地区数值反演试验与熵减分析结果表明,可见光红通道和近红外通道的信息指数分别为11.682 2和10.072 6;基于L曲线的地物反射特性正则化反演前后,两个通道的平均信息指数增量分别为0.440 0和0.354 6,最大增量为2.467 2和2.290 5。L曲线正则化方法的优势在于其不依赖于地表参数的先验知识,可有效解决卫星观测信息不足时地表参数的反演问题。
  • 图  1  地表反射特性正则化反演方法流程

    Figure  1.  Flowchart for Inverse Method of Regularized Surface Reflection Properties

    图  2  研究区域的地类分布

    Figure  2.  Land Cover Distribution over Study Area

    图  3  正则化参数的L曲线求解

    Figure  3.  Regularized Parameter Obtained Using L-Curve

    图  4  近红外波段正则化参数直方图分布

    Figure  4.  Histogram of Regularized Parameters in Near-Infrared Band

    图  5  近红外通道数据累积期间地表黑空反照率和白空反照率反演结果

    Figure  5.  Retrieved Near-Infrared Black-and White-Sky Albedos During Data Accumulation Period

    图  6  白空反照率直方图统计结果

    Figure  6.  Histogram of Near-Infrared White-Sky Albedo

    图  7  BRDF模型反演前后信息熵减情况

    Figure  7.  Information Entropy Reduction Before and After BRDF Model Inversion

    表  1  反演区域地物覆盖情况

    Table  1.   Land Cover of Inversion Area

    编号 MODIS分类名称 像元数 占总像元数的百分比/%
    0 水体 63 365 17.601 4
    1 常绿针叶林 1 256 0.348 9
    2 常绿阔叶林 0 0
    3 落叶针叶林 221 0.061 4
    4 落叶阔叶林 18 0.005 0
    5 混交林 2 267 0.629 7
    6 郁闭灌丛 2 564 0.712 2
    7 开放灌丛 697 0.193 6
    8 多树草原 1 532 0.425 6
    9 稀树草原 115 0.031 9
    10 草地 6 751 1.875 3
    11 永久湿地 1 020 0.283 3
    12 耕地 251 873 69.964 7
    13 城市与建筑用地 19 749 5.485 8
    14 耕地与自然植被混交地 2 550 0.708 3
    15 雪与冰 7 0.001 9
    16 裸地与稀疏植被 6 015 1.670 8
    下载: 导出CSV
  • [1] 肖登攀, 陶福禄, Moiwo J P.全球变化下地表反照率研究进展[J].地球科学进展, 2011, 26(11):1217-1224 http://www.cnki.com.cn/Article/CJFDTOTAL-DXJZ201111015.htm

    Xiao Dengpan, Tao Fulu, Moiwo J P. Research Progress on Surface Albedo Under Global Change[J]. Advances in Earth Science, 2011, 26(11):1217-1224 http://www.cnki.com.cn/Article/CJFDTOTAL-DXJZ201111015.htm
    [2] 刘荣高, 刘洋, 刘纪远. MODIS科学数据处理进展[J].自然科学进展, 2009, 19(2):141-147 doi:  10.3724/SP.J.1047.2016.00564

    Liu Ronggao, Liu Yang, Liu Jiyuan. Advances in MODIS Scientific Data Processing[J]. Progress in Natural Science, 2009, 19(2):141-147 doi:  10.3724/SP.J.1047.2016.00564
    [3] Li X, Gao F, Wang J, et al. A Priori Knowledge Accumulation and Its Application Linear BRDF Model Inversion[J]. J Geophys Res, 2001, 106(D11):11925-11936 doi:  10.1029/2000JD900639
    [4] 王彦飞, Ma Shiqian, 杨华, 等.带先验约束的地表参数提取的有效反演方法[J].中国科学D辑:地球科学, 2009, 39(3):360-369 http://www.oalib.com/paper/4153214

    Wang Yanfei, Ma Shiqian, Yang Hua, et al. On the Effective Inversion by Imposing a Priori Information for Retrieval of Land Surface Parameters[J]. Science in China:Series D Earth Sciences, 2009, 39(3):360-369 http://www.oalib.com/paper/4153214
    [5] 崔生成, 杨世植, 乔延利, 等.基于偏差原则和正则化方法反演晴空地表BRDF和反照率[J].武汉大学学报·信息科学版, 2011, 36(4):498-503 http://ch.whu.edu.cn/CN/Y2011/V36/I4/498

    Cui Shengcheng, Yang Shizhi, Qiao Yanli, et al. Retrieval of Clear-Sky Land Surface BRDF and Albedo Based on Discrepancy Principle and Regularization Method[J]. Geomatics and Information Science of Wuhan University, 2011, 36(4):498-503 http://ch.whu.edu.cn/CN/Y2011/V36/I4/498
    [6] Qin J, Yan G, Liu S, et al. Application of Ensemble Kalman Filter to Geophysical Parameters Retrieval in Remote Sensing:A Case Study of Kernel-Driven BRDF Model Inversion[J]. Science in China:Series D:Earth Sciences, 2006, 49(6):632-640 doi:  10.1007/s11430-006-0632-x
    [7] 陈爱军, 卞林根, 刘玉洁, 等. MODIS反演多云地区地表反照率的一种新方法及其反演试验[J].大气科学, 2007, 31(5):855-862 http://www.cqvip.com/qk/91836X/200705/25381846.html

    Chen Aijun, Bian Lingen, Liu Yujie, et al. A New Way to Retrieve Surface Albedo over Cloudy Areas with MODIS Data and the Retrieval Experiment[J]. Chinese Journal of Atmospheric Sciences, 2007, 31(5):855-862 http://www.cqvip.com/qk/91836X/200705/25381846.html
    [8] 王开存, 王建凯, 王普才, 等.用MODIS反演北京城市地区地表反照率精度以及算法改进[J].大气科学, 2008, 32(1):67-74 http://www.oalib.com/paper/4367960

    Wang Kaicun, Wang Jiankai, Wang Pucai, et al. Inversion of Accuracy of MODIS Albedo over Beijing Urban Area and Its Algorithm Improvement[J]. Chinese Journal of Atmospheric Sciences, 2008, 32(1):67-74 http://www.oalib.com/paper/4367960
    [9] 崔生成, 杨世植, 乔延利, 等.正则化滤波在反演地表光谱反照率中的应用[J].武汉大学学报·信息科学版, 2012, 37(6):653-657 http://ch.whu.edu.cn/CN/Y2012/V37/I6/653

    Cui Shengcheng, Yang Shizhi, Qiao Yanli, et al. Application of Regularized Filtering Technique to Retrieval of Land Surface Spectral Albedo[J]. Geomatics and Information Science of Wuhan University, 2012, 37(6):653-657 http://ch.whu.edu.cn/CN/Y2012/V37/I6/653
    [10] Pokrovsky I, Pokrovsky O, Roujean J L. Land Surface Albedo Retrieval via Kernel-Based BRDF Mo-deling. Part Ⅰ:Statistical Inversion Method and Mo-del Comparison[J]. Remote Sens Environ, 2002, 84:100-119
    [11] Luo Y, Trishchenko A P, Latifovic R, et al. Surface Bidirectional Reflectance and Albedo Properties Derived Using a Land Cover-Based Approach with Moderate Resolution Imaging Spectroradiometer Observations[J]. J Geophys Res, 2005, 110:D01106, doi:10.1029l2004JD004741
    [12] Quaife T, Lewis P. Temporal Constraints on Linear BRF Model Parameters[J]. IEEE Tran Geosci Remote Sens, 2010, 48(5):2445-2450 doi:  10.1109/TGRS.2009.2038901
    [13] Roujean J L, Leroy M, Deschanps P Y. A Bidirectional Reflectance Model of the Earth's Surface for the Correction of Remote Sensing Data[J]. J Geophys Res, 1992, 97:20455-20468 doi:  10.1029/92JD01411
    [14] Tikhonov A. On the Solution of Incorrectly Stated Problems and a Method of Regularization[J]. Dokl Akad Nauk SSSR, 1963, 151:501-504
    [15] Cui S, Yang S, Qiao Y, et al. Adaptive Regularized Filtering for BRDF Model Inversion and Land Surface Albedo Retrieval Based on Spectrum Cutoff Technique[J]. Optik, 2012, 123:250-256 doi:  10.1016/j.ijleo.2011.04.008
    [16] Hansen P, O'Leary D. The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[J]. SIAM J Sci Comput, 1993, 14(6):1487-1503 doi:  10.1137/0914086
    [17] Cui S, Wang Z, Yang S. Parameterization of Land Surface Albedo[J]. Chin Opt Lett, 2014, 12(11):1-5
    [18] Wang Y, Li X, Nashed Z, et al. Regularized Kernel-Based BRDF Model Inversion Method for Ill-Posed Land Surface Parameter Retrieval[J]. Remote Sens Environ, 2007, 111:36-50 doi:  10.1016/j.rse.2007.03.007
  • [1] 徐新强, 赵俊.  航空重力向下延拓的多参数正则化方法 . 武汉大学学报 ● 信息科学版, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
    [2] 李建, 周屈, 陈晓玲, 田礼乔, 李亭亭.  近岸/内陆典型水环境要素定量遥感空间尺度问题研究 . 武汉大学学报 ● 信息科学版, 2018, 43(6): 937-942. doi: 10.13203/j.whugis20160174
    [3] 杨梦诗, 廖明生, 秦晓琼, 史绪国.  C和L波段SAR数据在填海新区的应用及特性分析 . 武汉大学学报 ● 信息科学版, 2017, 42(9): 1300-1305. doi: 10.13203/j.whugis20150356
    [4] 曾小牛, 刘代志, 李夕海, 苏娟, 陈鼎新, 齐玮.  一种改进的病态问题奇异值修正法 . 武汉大学学报 ● 信息科学版, 2015, 40(10): 1349-1353. doi: 10.13203/j.whugis20130709
    [5] 崔生成, 杨世植, 赵强, 王震.  正则化滤波在反演地表光谱反照率中的应用 . 武汉大学学报 ● 信息科学版, 2012, 37(6): 653-657.
    [6] 崔生成, 杨世植, 乔延利, 赵强.  基于偏差原则和正则化方法反演晴空地表BRDF和反照率 . 武汉大学学报 ● 信息科学版, 2011, 36(4): 498-503.
    [7] 时向勇, 李先华, 郑成建.  基于椭圆曲线密码体制的遥感图像加密算法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1309-1313.
    [8] 王泽民, 邱蕾, 孙伟, 王海军.  GPS现代化后L_2载波的定位精度研究 . 武汉大学学报 ● 信息科学版, 2008, 33(8): 779-782.
    [9] 赵红蕊, 唐中实, 李小文.  非线性不适定遥感反演中正则化参数的定量确定 . 武汉大学学报 ● 信息科学版, 2008, 33(6): 577-580.
    [10] 闫利, 岳昔娟, 崔晨风.  一种定量确定遥感融合图像空间分辨率的方法 . 武汉大学学报 ● 信息科学版, 2007, 32(8): 667-670.
    [11] 赵红蕊, 唐中实, 李小文.  线性正则化遥感反演中正则化参数的确定方法 . 武汉大学学报 ● 信息科学版, 2007, 32(6): 531-535.
    [12] 刘耀林, 罗志军.  基于GIS的小流域水土流失遥感定量监测研究 . 武汉大学学报 ● 信息科学版, 2006, 31(1): 35-38.
    [13] 吕恒, 李新国, 曹凯.  基于BP神经网络模型的太湖悬浮物浓度遥感定量提取研究 . 武汉大学学报 ● 信息科学版, 2006, 31(8): 683-686.
    [14] 陈雪, 马建文.  遥感影像边缘增强的L_2向量空间范数算法 . 武汉大学学报 ● 信息科学版, 2006, 31(6): 489-491.
    [15] 孙春生, 吴军.  基于“刀刃”曲线特性的航空影像纹理视觉改善研究 . 武汉大学学报 ● 信息科学版, 2005, 30(3): 255-259.
    [16] 王振杰, 欧吉坤, 曲国庆, 韩保民.  L-曲线法确定半参数模型中的平滑因子 . 武汉大学学报 ● 信息科学版, 2004, 29(7): 651-653.
    [17] 王振杰, 欧吉坤.  L-曲线法确定岭估计中的岭参数 . 武汉大学学报 ● 信息科学版, 2004, 29(3): 235-238.
    [18] 徐天河, 杨元喜.  抗差Tikhonov正则化方法及其应用 . 武汉大学学报 ● 信息科学版, 2003, 28(6): 719-722.
    [19] 童小华, 刘大杰.  道路曲线数字化数据的联合平差模型 . 武汉大学学报 ● 信息科学版, 2001, 26(1): 64-69.
    [20] 刘文宝, 黄幼才.  地图曲线数字化误差估计 . 武汉大学学报 ● 信息科学版, 1994, 19(4): 352-358.
  • 加载中
图(7) / 表(1)
计量
  • 文章访问数:  999
  • HTML全文浏览量:  89
  • PDF下载量:  287
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-04-10
  • 刊出日期:  2018-08-05

基于正则化的地表反射特性参数反演方法

doi: 10.13203/j.whugis20160244
    基金项目:

    国家自然科学基金 41305019

    安徽省自然科学基金 1308085QD70

    作者简介:

    崔生成, 博士, 副研究员, 主要从事地表反射特性定量遥感理论与反演方法研究。csc@aiofm.ac.cn

    通讯作者: 李学彬, 博士, 副研究员。xbli@aiofm.ac.cn
  • 中图分类号: P237;P407.9

摘要: 区域和全球尺度下研究气候变化,需要准确描述大气边界层处下垫面地表的光反射特性。为此,提出基于正则化的约束反演方法,其核心在于正则化参数的选取。通过L曲线拐点位置的确定,求出最合适的正则化约束参数,以稳定陆表反射特性参数的定量反演。京津唐地区数值反演试验与熵减分析结果表明,可见光红通道和近红外通道的信息指数分别为11.682 2和10.072 6;基于L曲线的地物反射特性正则化反演前后,两个通道的平均信息指数增量分别为0.440 0和0.354 6,最大增量为2.467 2和2.290 5。L曲线正则化方法的优势在于其不依赖于地表参数的先验知识,可有效解决卫星观测信息不足时地表参数的反演问题。

English Abstract

崔生成, 朱文越, 杨世植, 李学彬. 基于正则化的地表反射特性参数反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
引用本文: 崔生成, 朱文越, 杨世植, 李学彬. 基于正则化的地表反射特性参数反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
CUI Shengcheng, ZHU Wenyue, YANG Shizhi, LI Xuebin. Regularization-Based Retrieval Method for Surface Reflective Property Parameters[J]. Geomatics and Information Science of Wuhan University, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
Citation: CUI Shengcheng, ZHU Wenyue, YANG Shizhi, LI Xuebin. Regularization-Based Retrieval Method for Surface Reflective Property Parameters[J]. Geomatics and Information Science of Wuhan University, 2018, 43(8): 1264-1270. doi: 10.13203/j.whugis20160244
  • 气候模型需要在不同尺度下准确描述大气边界层处下垫面地表对太阳入射辐射的反射作用[1]。地表反射特性参数的定量遥感关键在于地物二向性反射分布函数(bidirectional reflectance distribution function, BRDF)模型的反演。在太阳反射波段,云严重限制了地物方向性反射信息的累积,再加上遥感信号通常存在不同程度的相关性[2-3],使得地物晴空有效观测的数目不足, 角度分布较弱,线性核驱动BRDF模型会表现出离散的不适定性[4-5]

    针对BRDF模型的不适定反演问题,Li等[3]基于贝叶斯理论,使用先验地物波谱构建“潜数据”,首次实现了单次观测反演地物BRDF和反照率。Qin等[6]利用地物波谱数据产生的先验知识,成功使用扩展卡曼滤波算法降低了中分辨率成像光谱仪(moderate-resolution imaging spectro radiometer,MODIS)观测信息不足时BRDF模型反演的不确定性。但基于先验知识的反演方法对波谱数据样本的地类代表性有一定的要求,其反演结果的精度依赖于对波谱知识的更新。陈爱军等[7]利用地表分类信息和同步的归一化植被指数(normolized difference vegetation index, NDVI)数据,通过同类像元的方向性观测构建缺测像元的多角度观测,有效实现了MODIS晴空观测不足时青藏高原地区地表反照率的反演。该方法考虑了地类空间异质性的特点,在短时相内样本充足时其精度有可靠保障,但逐像元寻找同类地物的晴空有效观测,则存在效率方面的问题。王开存等[8]利用归一化的BRDF模型,实现了观测量不少于3个时城市地表反照率的稳定反演,但不适用于少于3个观测量的欠定情况。王彦飞等[4]提出了使用正则化手段来克服由于MODIS观测不足引起的BRDF模型的离散不适定性,对陆表参数的稳定反演起到了促进作用。针对正则化的反演算法,双参数模型的采用[5]以及依据传感器信噪比方法的提出[9],有效解决了正则化参数的初始取值问题。

    国外研究中,Pokrovsky等[10]使用欧洲第二代地球静止轨道气象卫星MSG(meteosat second generation)数据,统计了Roujean BRDF模型在全球不同纬度带反演的病态指数,建立了基于优化观测数据的统计反演方法。Luo等[11]研究了不同植被的BRDF和NDVI之间的关系,建立了基于NDVI信息辅助下的短时相MODIS BRDF反演方法。Quaife和Lewis[12]使用拉格朗日乘子实现了任意时间窗口内地物BRDF的反演,并实现了初步区域应用。

    无论从观测数据本身建立的优化策略或同化多种数据构造多角度观测,还是从算法本身建立的贝叶斯方法、卡尔曼滤波方法以及正则化方法,目的都在于增加数据反演所需的有效信息量。正则化方法是解决卫星数据信息量不足时地表反射特性遥感反演的一种有效途径,其算法的核心在于正则化参数的确立。目前地表参数反演过程中正则化参数的选取主要基于经验选取[12]或需要事先知道数据的噪声水平[9]。前者过于粗糙, 后者则通常难以知晓。因此,数据噪声水平未知情况下的正则化算法更具有普适性。结合国内外BRDF的研究现状和发展趋势,本文提出基于L曲线法反演地表反射特性参数的思路,建立卫星观测不足时地表BRDF模型正则化反演方法,并通过区域反演实例分析反演前后信息熵减变化情况,以验证构建方法的有效性。

    • 地表对入射光场的反射特性可用BRDF模型描述[13]

      $$ \rho \left( {{\theta _i},{\theta _o},\varphi ;\lambda } \right) = \sum\limits_k {{f_k}\left( \lambda \right){K_k}\left( {{\theta _i},{\theta _o},\varphi } \right)} $$ (1)

      式中, λ为中心波长; 相对方位角φ=|φo -φi |; (θi, φi)、(θo, φo)分别表示入射和出射方向; K表示散射核函数; 下标k=0, 1, 2分别对应各向同性散射核、几何光学散射核以及体散射核; f为模型参数。

      针对观测样本数量不足、角度分布较弱的情况,Phillips-Tikhonov(PT)正则化方法通过引入最小化模型参数f的范数‖f‖作为边约束条件,实现模型参数的稳定求解[14]

      $$ {\mathit{\boldsymbol{f}}_\gamma } = \min \left( {{{\left\| {\mathit{\boldsymbol{Kf}} - \mathit{\boldsymbol{\rho }}} \right\|}^2} + \gamma {{\left\| \mathit{\boldsymbol{f}} \right\|}^2}} \right) $$ (2)

      式中,K为BRDF模型核函数矩阵算子;ρ为观测向量;γ为正则化参数,调节遥感模型反演过程中最小二乘和‖ f ‖范数最小化两者之间信息流的平衡问题。当卫星观测数据的有效信息量不足以支撑模型参数的成功反演时,正则化方法通过约束机制加入额外的信息,促进信息流向模型成功反演的方向,保证在解空间中找到最合适的模型参数。

      基于对核矩阵K的奇异值分解(singular value decomposition,SVD), PT方法的正则化解为[15]

      $$ {\mathit{\boldsymbol{f}}_\gamma } = \sum {\frac{{\sigma _i^2}}{{\sigma _i^2 + {\gamma ^2}}}\frac{{\left( {\mathit{\boldsymbol{u}}_i^{\rm{T}}\mathit{\boldsymbol{\rho }}} \right)}}{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i}} $$ (3)

      式中,σi为SVD分解得到的奇异值; uivi分别是构成测量空间和模型参数空间的基。当卫星观测数据的角度分布弱, 观测样本不足时,最小二乘解的代数谱会出现不连续的现象,表现为集合{σi}中存在趋于0的微小元素簇,导致最小二乘解对数据信息量不足表现出来的噪声非常敏感,而PT方法则通过滤波算子$ {\psi _i}\left( \gamma \right){\rm{ }} = \frac{{\sigma _i^2}}{{\sigma _i^2 + {\gamma ^2}}}$来抑制数据噪声向模型参数的传播。因此正则化参数γ的作用非常关键。

    • L曲线法实质上是寻找残差范数与解范数的“拐点”,解范数和残差范数分别为:

      $$ \xi \left( \gamma \right) = {\left\| {{\mathit{\boldsymbol{f}}_\gamma }} \right\|^2} $$ (4)
      $$ \xi \left( \gamma \right) = {\left\| {\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{f}}_\gamma } - \mathit{\boldsymbol{\rho }}} \right\|^2} $$ (5)

      L曲线拐点的确定需要计算曲线各点的曲率,曲率最大的点即为所求的“拐点”[16]

      $$ l\left( \gamma \right) = \frac{{\hat \xi '\hat \delta '' - \hat \xi ''\hat \delta '}}{{{{\left( {{{\hat \xi '}^2} + {{\hat \delta '}^2}} \right)}^{3/2}}}} $$ (6)

      式中,$\hat \xi = {\rm{ln}}\xi \left( \gamma \right){\rm{ }}, \hat \delta = {\rm{ln}}\delta \left( \gamma \right) $。为了易于数值实现,此处给出L曲线曲率的计算式:

      $$ \begin{array}{l} l\left( \gamma \right) = \frac{{\xi \left( \gamma \right)\delta \left( \gamma \right)}}{{\xi '\left( \gamma \right)}} \times \\ \frac{{{\gamma ^2}\xi '\left( \gamma \right)\delta \left( \gamma \right) + 2\gamma \xi \left( \gamma \right)\delta \left( \gamma \right) + {\gamma ^4}\xi \left( \gamma \right)\xi '\left( \gamma \right)}}{{{{\left( {{\gamma ^4}{\xi ^2}\left( \gamma \right) + {\delta ^2}\left( \gamma \right)} \right)}^{3/2}}}} \end{array} $$ (7)

      从作用机制上来说,L曲线将作用空间划分为两个区域,正则化的解只能位于L曲线上或曲线上方。当正则化约束很小时,数据噪声是整体误差的主要来源,此即欠约束现象,对应于L曲线的垂直部分;而当引入大量的正则化约束时,光滑约束误差是总误差的主要来源,称为过约束,对应于L曲线的水平部分。因此,最优化的正则化参数必然在L曲线拐点附近。

    • 将BRDF在所有出射方向上积分,可得到黑空反照率αBS;进一步将黑空反照率在所有入射方向上积分,便可得到白空反照率αWS

      $$ \begin{array}{l} {\alpha _{{\text{BS}}}}\left( {\lambda ,{\theta _s}} \right) = \sum\limits_k {{f_k}\left( \lambda \right){h_k}\left( {{\theta _s}} \right)} = \hfill \\ \sum\limits_k {{f_k}\left( \lambda \right)\iint\limits_{2{\text{\pi }}} {{K_k}\left( {{\theta _s},{\theta _o},\varphi } \right)\cos {\theta _o}\sin {\theta _o}{\text{d}}{\theta _o}{\text{d}}\varphi }} \hfill \\ \end{array} $$ (8)
      $$ \begin{array}{l} {\alpha _{{\rm{WS}}}}\left( {\lambda ,{\theta _s}} \right) = \sum\limits_k {{f_k}\left( \lambda \right){H_k}} = \\ \sum\limits_k {{f_k}\left( \lambda \right)\int_0^{\pi /2} {{h_k}\left( {{\theta _i}} \right)\cos {\theta _i}\sin {\theta _i}{\rm{d}}{\theta _i}} } \end{array} $$ (9)

      对于线性核驱动BRDF模型[3, 13],可按式(8)~(9)预先计算核函数的角度积分,按固定角度间隔步长生成查找表[17],从而避免逐像元四重积分计算,大大提高了运算效率。

    • 基于L曲线的正则化反演方法,其整体技术流程如图 1所示。

      图  1  地表反射特性正则化反演方法流程

      Figure 1.  Flowchart for Inverse Method of Regularized Surface Reflection Properties

      算法实现步骤如下。

      1) 使用MODIS MOD09A1数据,获取样本累积窗区内的有效观测;

      2) 根据太阳和卫星观测天顶角与方位角信息,分别计算不同观测时间的BRDF模型核函数,生成核矩阵;

      3) 从MOD09A1的状态标记(state flag,SF)数据获取质量评估(quality assessment,QA)信息,逐像元提取每次观测的质量标识,计算单个像元的综合质量指示因子,判识有效观测像元;

      4) 运用L曲线确定正则化参数,同时确定正则化光滑矩阵,从而实现利用有效观测稳定反演BRDF模型,获取BRDF模型参数;

      5) 对BRDF核函数进行角度积分运算,生成BRDF积分核查找表。利用插值方法,获得观测时刻的地表反照率。

    • 利用2008年11月16日和24日(分别对应年积日(day of year, DOY)321和329)2 d的MOD09A系列数据,对京津唐部分地区进行了陆表反射特性参数的反演试验。MOD09A1是经过去云处理与大气校正后的地表反射率产品。研究区域地表覆盖类型分布情况如图 2所示,基于MODIS MCD12Q1地类数据产品,相应的分类统计信息见表 1

      图  2  研究区域的地类分布

      Figure 2.  Land Cover Distribution over Study Area

      表 1  反演区域地物覆盖情况

      Table 1.  Land Cover of Inversion Area

      编号 MODIS分类名称 像元数 占总像元数的百分比/%
      0 水体 63 365 17.601 4
      1 常绿针叶林 1 256 0.348 9
      2 常绿阔叶林 0 0
      3 落叶针叶林 221 0.061 4
      4 落叶阔叶林 18 0.005 0
      5 混交林 2 267 0.629 7
      6 郁闭灌丛 2 564 0.712 2
      7 开放灌丛 697 0.193 6
      8 多树草原 1 532 0.425 6
      9 稀树草原 115 0.031 9
      10 草地 6 751 1.875 3
      11 永久湿地 1 020 0.283 3
      12 耕地 251 873 69.964 7
      13 城市与建筑用地 19 749 5.485 8
      14 耕地与自然植被混交地 2 550 0.708 3
      15 雪与冰 7 0.001 9
      16 裸地与稀疏植被 6 015 1.670 8

      图 3给出了Ross Thick-Roujean BRDF模型[13]参数的L曲线求解过程,包括水平和垂直两部分。水平部分对应正则化参数γ较大时的情况。γ减小时,残差相应变小。当L曲线达到垂直部分后,继续减小γ会对解的范数产生较大影响。当算法沿垂直方向和水平方向寻优到交汇点时的正则化参数是最优的,即所有观测数据中的有效信息均充分利用且噪声影响最小。利用此时的目标函数得到的解即为最佳的模型参数解。

      图  3  正则化参数的L曲线求解

      Figure 3.  Regularized Parameter Obtained Using L-Curve

      图 4以直方图形式统计了研究区域正则化参数的分布情况。可以看出,BRDF模型反演过程中,边约束权重多分布在0.2附近;边约束权重分布在1.6附近对应的卫星观测数据, 因卫星观测信息不足,需要更多的信息进行光滑约束。不同地物材质、结构和形状等不同,其BRDF形状不同,利用卫星数据反演BRDF模型所需的观测信息具有空间异质性的特点,导致边约束引入的附加信息也不是空间均一的。因此,对整景卫星图像如只使用相同的正则化参数,实质上忽视了地物空间异质性的特征,不可避免地会给反演带来不同程度的误差。

      图  4  近红外波段正则化参数直方图分布

      Figure 4.  Histogram of Regularized Parameters in Near-Infrared Band

      图 5给出了MODIS近红外通道地表反射特性参数反演结果,由于黑空反照率存在强烈的太阳天顶角依赖特性,图 5(a)5(b)分别对应DOY 321、DOY 329的MODIS/Terra过境时刻(地方时10:30)的黑空地表反照率,图 5(c)为数据累积期间地表的白空反照率。表 1中,研究区域涵盖了15种MODIS国际地圈-生物圈计划(international geosphere-biosphere program, IGBP)地物,其中以耕地为主。结合图 6的白空反照率直方图分布情况不难发现,耕地在可见-近红外谱段的白空反照率分别分布在0.03和0.05附近。由于试验区域位于我国北方,冬季太阳高度角小,加上卫星观测数不足(≤3),无论是基于最小二乘的MODIS全反演法,还是基于全球先验BRDF模型参数储备库的量反演法,均无法给出反演结果,而基于L曲线的正则化反演方法有效实现了模型参数的稳定反演。

      图  5  近红外通道数据累积期间地表黑空反照率和白空反照率反演结果

      Figure 5.  Retrieved Near-Infrared Black-and White-Sky Albedos During Data Accumulation Period

      图  6  白空反照率直方图统计结果

      Figure 6.  Histogram of Near-Infrared White-Sky Albedo

    • 基于信息熵的概念,此处引入信息指数(information index, Ⅱ),用于衡量反演过程中的有效信息量。

      $$ {\rm{II}} = \ln \left( {\prod\limits_i {{\sigma _i}} } \right) - \ln \delta = \sum\limits_i {\ln {\sigma _i} - \ln \delta } $$ (10)

      信息熵减则可用于分析模型反演前后或不同方法之间有效信息量的变化情况:

      $$ \Delta {\rm{II}} = {\rm{II}} - {\rm{I}}{{\rm{I}}_{{\rm{REF}}}} $$ (11)

      鉴于谱截断奇异值分解(spectral truncated single value decomposition, STSVD)方法反演陆表反射特性参数的有效性[9, 18],以STSVD方法的信息指数作为参考ⅡREF,分析本文提出的L曲线PT正则化算法(L-curve Phillips-Tikhonov regularization, LPTR)的信息熵变化情况。

      图 7给出了Roujean模型反演前后的信息指数。作为一个算例,对同一像元分别使用STSVD和LPTR方法。结果表明,STSVD算法在红色可见光和近红外波段的信息指数分别为9.450 1和7.837 4,信息量增益为2.232 1;而LPTR算法则分别为11.682 2和10.072 6,相应的信息量增益为2.235 2。

      图  7  BRDF模型反演前后信息熵减情况

      Figure 7.  Information Entropy Reduction Before and After BRDF Model Inversion

      从模型反演前后信息熵变化的角度出发,LPTR方法对欠定模型的反演具有优势,主要原因在于:①相对于经验取值法,LPTR正则化参数的L曲线优化选取充分考虑了像元地类分布的空间非均一性(图 4);②即使是同一类地物像元,其观测角度的几何分布不同以及通过质量控制标识可能会过滤掉部分观测样本,导致正则化参数的优化程度不同。针对存在两次观测的像元,利用L曲线和广义交叉验证(generalized cross validation, GCV)方法得到的正则化参数分别为0.085 7和0.045 8。这意味着单纯从正则化参数无法整体考量算法的作用。不同的正则化参数选择方法的直接效果是体现在地表反射特性参数反演前后信息量的增减情况上。因此,需要从分析BRDF模型反演前后信息熵的变化来评价整个算法。对整个研究区域而言,红色可见光波段与近红外波段的平均信息指数增量分别为0.440 0和0.354 6,最大增量可达2.467 2和2.290 5。

      图 7可以发现,研究区域内有两个像元在反演前后信息熵减较大,分别为-50.76和-58.17。经分析发现,这两个像元的地类分别是耕地与城市建筑用地,但根据质量标识发现这两个像元在卫星过境时刻均为短暂的水(ephemeral water)像元,且气溶胶含量分别为低和平均程度。造成这种现象的一种可能原因是像元累积期间降雨致地表覆盖表层短期积水,为了达到与MODIS算法结果比较的一致性,本文算法处理的样本中包含了临时水像元观测。

      综上可见,LPTR算法通过边约束的加入,有效平衡了边约束过程中信息的增加和额外噪声引入的信息损失,能有效解决卫星观测信息不足带来的不稳定反演问题。

    • 为了解决地表反射特性参数反演时卫星观测信息量不足的共性难题,本文针对数据噪声水平事先未知的实际情况,提出了基于L曲线的地表反射特性参数正则化反演方法。通过对BRDF模型残差范数与解范数“拐点”的寻优来确立最佳的正则化参数。同时充分顾及了地表分布的空间异质性特征,根据下垫面地表的类型、反射的方向性和光谱信息,自适应调节BRDF模型解约束的权重,使边约束信息和实际的观测信息向有利于模型稳定求解的方向流动。

      该方法通过对模型核矩阵算子的代数谱结构进行优化改善,有效增加了地表反射特性参数成功反演所需的信息量,在处理病态反演问题时,具有独特的算法优势,从而为缺乏信息支持下地表BRDF模型参数的稳定求解问题提供了新的解决途径。后续的工作重点将围绕利用单次观测数据反演地表反射特性参数展开,以满足单视场高空间分辨率卫星数据处理的业务需求。

参考文献 (18)

目录

    /

    返回文章
    返回