发明名称 一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法
摘要 本发明涉及一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法。通过应用BRDF原型对多角度观测数据进行循环迭代拟合,求出最小均方根误差(RMSE),提出基于BRDF原型反演地表反照率的新方法,以黑河地区机载红外广角双模式成像仪(WiDAS)实验数据为例,应用BRDF原型作为二向性反射先验知识,约束大观测天顶角下二向性反射形状的剧烈变化,克服极端反射率变化所引起的反照率的不确定性。该方法不同于目前业务化运行的基于核驱动模型反演地表反照率的反演算法,为空间采样能力不足的传感器提供了一种反演地表反照率的解决方案。本发明在空间信息技术领域,尤其在定量遥感方面具有重要的应用价值。
申请公布号 CN103413014A 申请公布日期 2013.11.27
申请号 CN201310075882.3 申请日期 2013.03.11
申请人 北京师范大学 发明人 焦子锑;张虎
分类号 G06F19/00(2011.01)I 主分类号 G06F19/00(2011.01)I
代理机构 代理人
主权项 一种基于二向性反射分布函数(BRDF)原型反演地表反照率的方法,步骤包括:步骤一:针对各原型求模型和观测的最小拟合误差对于稀疏多角度采样数据,BRDF空间采样不足会导致多角度观测信息量不足,因此,用全模型反演地表反照率,由于在大角度观测方向缺少约束,会高估或低估反照率,BRDF原型提供了一种有效约束,通过对稀疏采样数据用BRDF原型进行约束拟合,可综合考虑传感器实际观测和地表各向异性反射变化,拟合过程可表示如下:假定一组有n个观测的多角度观测数据B:B=ρ0(Ωi,0;Ωr,0),ρ1(Ωi,1;Ωr,1)...ρn‑1(Ωi,n‑1;Ωr,n‑1),若与其对应的BRDF原型为BRDF′,对应相同观测位置的原型反射率数据表示为B′:B′=ρ′0(Ωi,0;Ωr,0),ρ′1(Ωi,;Ωr,1)...ρ′n‑1(Ωi,n‑1;Ωr,n‑1),取调整系数a,通过平移BRDF原型,使得B和aB′的差异最小,采用最小二乘拟合法最小化拟合误差e2 <mrow> <msup> <mi>e</mi> <mn>2</mn> </msup> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>&rho;</mi> <mi>i</mi> </msub> <mo>-</mo> <mi>a</mi> <msub> <msup> <mi>&rho;</mi> <mo>&prime;</mo> </msup> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>求调整系数a <mrow> <mi>a</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>&rho;</mi> <mi>j</mi> </msub> <mo>&times;</mo> <msub> <msup> <mi>&rho;</mi> <mo>&prime;</mo> </msup> <mi>j</mi> </msub> <mo>/</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <msub> <msup> <mi>&rho;</mi> <mo>&prime;</mo> </msup> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>这样,数据B的BRDF可表示为BRDF=a×BRDF′       (3)拟合误差RMSE可以表示为 <mrow> <mi>RMSE</mi> <mo>=</mo> <msqrt> <mfrac> <mn>1</mn> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>&rho;</mi> <mi>k</mi> </msub> <mo>-</mo> <mi>a</mi> <msub> <msup> <mi>&rho;</mi> <mo>&prime;</mo> </msup> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>步骤二:选取最优原型求其光谱固有反照率用6个原型分别拟合观测数据,计算各原型的拟合误差RMSE,选择拟合观测数据RMSE最小的原型作为数据的先验BRDF,求该原型的方向‑半球和半球‑半球积分如下: <mrow> <msub> <mi>h</mi> <mrow> <mi>BRDF</mi> <mo>&prime;</mo> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mi>&pi;</mi> </mfrac> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </msubsup> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mrow> <mi>&pi;</mi> <mo>/</mo> <mn>2</mn> </mrow> </msubsup> <mi>BRDF</mi> <mo>&prime;</mo> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>v</mi> </msub> <mo>,</mo> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>&Delta;&phi;</mi> <mo>)</mo> </mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>v</mi> </msub> <mo>)</mo> </mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>v</mi> </msub> <mo>)</mo> </mrow> <mi>d</mi> <msub> <mi>&theta;</mi> <mi>v</mi> </msub> <mi>d&phi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>对hBRDF′(θs)在太阳入照半球进一步积分,得半球‑半球积分 <mrow> <msub> <mi>H</mi> <mrow> <mi>BRDF</mi> <mo>&prime;</mo> </mrow> </msub> <mrow> <mo>(</mo> <mi>&theta;</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>2</mn> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mrow> <mi>&pi;</mi> <mo>/</mo> <mn>2</mn> </mrow> </msubsup> <msub> <mi>h</mi> <mrow> <mi>BRDF</mi> <mo>&prime;</mo> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mi>sin</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mi>cos</mi> <mrow> <mo>(</mo> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mi>d</mi> <msub> <mi>&theta;</mi> <mi>s</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>通过式(3)确定调整系数a,则黑天空反照率(Black Sky Albedo,BSA)和白天空反照率(White SkyAlbedo,WSA)与选定原型的上述积分值之间的关系可以表示为:BSA(θ,λ)=a×hBRDF′(θ,λ)      (7)WSA(λ)=a×HBRDF′(λ)      (8)步骤三:产生地表宽波段真实反照率根据已获得的光谱固有反照率,即白天空反照率(WSA)和黑天空反照率(BSA),考虑到大气状态,通过测量天空散射光和直接光的比值,获得光谱真实反照率αblue(θi,λ)=SKYL(θi,λ)×WSA(λ)+(1‑SKYL(θi,λ))×BSA(θi,λ)    (9)通过建立WiDAS传感器波段与一个标准波段(如MODIS前7个波段)之间的关系,将WiDAS传感器的光谱反照率转换成标准波段的光谱返归率,通过窄波段到宽波段的转换公式,获得地表宽波段的真实反照率αshort=0.160α1+0.291α2+0.243α3+0.116α4(10)+0.112α5+0.081α7‑0.0015最后,通过地面验证和全模型无约束反演进行比较,验证算法精度。
地址 100875 北京市海淀区新街口外大街19号科技处