专利名称：一种基于辐射度的复杂地形区域辐射传输模拟方法
技术领域：
本发明涉及一种基于辐射度的复杂地形区域辐射传输模拟方法，适合于对山区复 杂地形区域条件下的遥感成像过程精确模拟，特别是定量分析地形对遥感影像的影响和影 像地形校正方面具有重要价值，属于遥感成像模拟技术领域。
背景技术：
在地形起伏情况下，地物将接收到周围物体反射的太阳直射光和大气散射光，导 致地物实际接收辐射能量的增加。在雪地、植被近红外波段以及陡峭坡面等情况下，周围 环境对目标物体反射的影响无法忽略。地物在接收周围物体反射能量的同时，也将部分能 量重新反射回周围物体，循环往复，形成复杂的地表多次散射过程，导致地物光谱信息的变 化，影响了山区遥感图像的信息提取精度，成为了限制山区遥感图像应用水平的主要障碍。 地表多次散射是地表太阳能计算以及遥感仿真中的一个难点问题。Proy等人(见引用文献1)提出了一个精确计算周围地物反射的方法，该算法仅考 虑单次散射影响，认为地物接收到的周围物体反射由周围所有可视物体的反射辐射叠加而 成。每个物体的反射辐射与其自身辐亮度、到目标地物的距离以及物体法线与两者连线的 夹角等因素有关。Proy算法不仅计算量大，而且周围物体辐亮度本身也是未知数，无法得到 广泛应用。Dozier和Frew(见引用文献2)利用周围物体反射的平均辐射亮度和地形可 视因子近似计算周围地物反射，该方法面临着与计算量大，周围物体辐亮度未知的问题。 Sandimerier等人(见引用文献3)提出了一个近似计算周围地物反射的公式，在该公式中， 周围地物反射仅与周围物体平均反射率、地形可视因子和水平地面上接收的总辐照度(包 括太阳直射和大气散射)有关。在假设坡面为无限长斜坡的时候，地形可视因子无需经过 繁杂的计算，该方法精度不高，但在实际中得到广泛的应用。计算复杂地形区域传感器接收到的辐射能量涉及众多因素，现存技术可归纳为两 大类一是缺乏严格物理意义的经验、半经验方法，特点是参数少、简便高效，缺点在于精度 低，适用性差；二是基于辐射传输原理的物理解析模型，其理论基础完善，模型参数具有明 确的物理意义，模拟精度较高，但模型往往对地形做了大量的简化近似，适用的范围有限。 复杂地形区域太阳光和地表相关作用是个复杂的非线性过程，现存技术忽略部分影响因素 或简化地形结构，造成模拟结果精度不高，致使模拟结果与遥感图像反射率之间出现较大 出入，严重影响了山区遥感图像的信息提取。引用文献1，Dozier, J.，James Frew. Rapid calculation of terrain parameters for radiationmodeling from digital elevation data[J] (由数字高禾呈数据快速计算用于 辐射建模的地形参数)IEEE Transactions on Geoscience and Remote Sensing，1990， 28(5) :963-969.2，Proy, C.，Tanre D.，Deschampls P. Y. Evaluation of topographic effectsin remotelysensed data[J].(对遥感数据中的地形效应进行评价分析)Remote Sense ofEnvironment,1989, 30 21.3, Sandmeier, S. , and Itten, K.I. A physically-based model to correct atmospheric andillumination effects in optical satellite data of rugged terrain [J],(一个用于对崎岖地形区光学遥感数据中的大气和太阳辐照效应进行校正的 物理模型)IEEETransactions on Geoscience and Remote Sensing, 1997, 35 :708_717.
发明内容
本发明的目的在于提供一种实现复杂地形状况下遥感辐射传输过程精确计算的 方法，克服现有技术的不足，基于复杂地形状况下辐射传输的物理过程，采用计算机模拟技 术，实现山区复杂地形条件下各部分辐射贡献的精确计算。本发明的技术解决方案是利用数字高程模型对地表进行建模，计算地表网格面 元的朝向，并计算地表网格面元之间的可视因子，根据太阳入射位置和能量分布状况计算 地表面元接收到的辐射通量密度，输入地表网格面元的波谱特性，建立研究区域的辐射度 方程，采用数值迭代的方法求解辐射度方程，得到每个地表网格面元的辐射通量密度，最后 根据传感器采样位置计算该复杂地形区域的方向反射率因子和辐亮度。本发明一种基于辐射度的复杂地形区域辐射传输模拟方法，具体步骤如下(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备；(2)输入太阳入射位置和能量分布状况等数据，计算每个地表网格面元接收到的 初始辐射通量密度；(3)使用计算机图形学的方法计算地表网格面元之间的可视因子；(4)输入地表网格面元的波谱特性数据，建立该复杂地形区域的辐射度方程；(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的 辐射通量密度；(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子，实现复杂地形 区遥感过程的精确模拟和计算。其中步骤(2)中，地表网格面元(以下简称面元)接收到的初始太阳辐射通量密 度(Fs)由如下公式计算Fs⑴=[rs/cos ( 0 s) ] sd| a(i, d)其中，rs为太阳入射辐射通量密度，0 s为太阳天顶角，叫为面元i的法线单位向 量，^为太阳方向单位向量，a(i, d)为该面元可以直接被太阳照射到的面积比例，Fs为该 面元接收到的初始辐射通量密度。其中，步骤(3)中利用计算机图形学中的半立方体算法计算面元之间的可视因 子，以计算场景中任意两个面元Ai和兴j时)之间的可视因子为例，具体计算过程如 下a.在面元&上微面元dAi中心处沿其法向量方向放置一个虚拟的半立方体，该半 立方体的五个表面被剖分成均勻正方形网格(通常剖分为100X100)，每一网格均对应dAi 朝向半球空间的一个微立体角，从而形成一个半空间立方体角查找表。预先计算，并存贮半 立方体中心微面元dAi对各表面网格的微可视因子，设q为半立方体表面上一个网格，微面元dAi对q的微可视因子为Fq ； b.以dAi为中心将面元A」投影到半立方体的五个表面上，计算投影区域包含的半 立方体表面网格，对表面网格对应的微可视因子进行叠加，计算出dAi对A」的可视因子，设 Q为投影区域包含的半立方体表面网格集合，则dAi对A」的可视因子可以表示为
F
dEi,Ej
q^Qc.上述步骤(a和b)计算中没有考虑面元之间的可视问题，对有遮挡的情形，若半 立方体表面网格被两个或者两个以上面元的投影区域所覆盖，则通过比较这些面元离半立 方体中心微面元的距离来决定在该像素处可见的面元，距离近的面元可见，距离远的面元 不可见；d.以a-c方法计算面元&上所有微面元dAi到A」的可视因子，则
cos <9,. cos 0JdAidAj
■ I i
A, Fv = FAr
aj
=-{ f :，——m
A JAi JAj 在步骤(3)中，假设地表为朗伯反射，则面元反射率(P》乘以面元接收到的初始 辐射通量密度(Fs)，可得该面元初始散射的辐射通量密度(E》
Ei = [rs/cos ( 0 s)] Irii sd|a(i, d) p 于是，可构造辐射度方程如下
N
jOi其中B” Bj分别为面元间达到辐射能量平衡后的面元i和j的辐射通量密度。以 矩阵形式表示为B = E+CB其中，B和E是1XN维向量，C是NXN维矩阵
B,)E,、B2E2f 0PXFU PAN、
B =E =c = 、BN〉、￡N)PN^N2 ...0 )步骤(5)中求解辐射度方程采用Gauss-Seidel迭代法，其基本过程为首先给定 一个初值B°，带入辐射度方程右边计算，得到一个新的值B1，然后再带入辐射度方程右边计 算，又可以得到一个新值B2，如此反复迭代，直到|Bn+1-Bn|小于一个允许的误差范围为止。 计算时，以初始散射的辐射通量密度(ED作为初值B°，误差范围一般选择|Bn+1-Bn|彡10_6 即可满足精度要求。迭代收敛后，可准确计算出三维场景中每个面元的辐射通量密度。步骤(6)中指定传感器观测方向，即可按下式计算该观测方向上的反射率因子 (Directional Reflectance Factor，DRF)
ln\ni sv\a{i,v)area{i)DRF(v) = 一w 丨 丨 、ty
2^,1/nyii. sv v)area\i)v为观测方向(以方向矢量sv表示)，Bi为面元i (以方向矢量&表示，其面积为
6area(i))的辐射通量密度，a(i，v)为该面元对观测方向上可见的面积比例。本发明是一种基于辐射度的复杂地形区域辐射传输模拟方法，与现有技术相比的 优点在于(1)本发明基于严格的辐射传输过程，采用计算机图形学中经典的辐射度方法来 对崎岖山区进行遥感反射率的精确模拟，纠正了现有技术中采用的各种近似和简化计算带 来的误差。(2)本发明对地形没有特殊的假定和要求，原则上任意复杂的地形区域，只要提供 该区域的数字高程模型数据和地表反射率数据，即可精确计算传感器测量得到的反射率， 大大扩展了山区辐射建模的应用范围，具有广泛的应用前景。
图1为本发明一种基于辐射度的复杂地形区域辐射传输模拟方法的流程图。
具体实施例方式以西藏驱龙地区遥感场景模拟过程为例，如图1所示，本发明的具体实施方法如 下(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备。使用西藏驱龙地区地面分辨率为30m的数字高程模型和地表反射率数据。模拟计 算时间设置为2002年4月3日上午10:6:56，区域中心地理坐标为29. 98° N，36.37° E。 根据数字高程模型数据提供的像元坐标信息，计算地表网格面元相邻两个边构成的向量的 法向量，确定该法向量的天顶角和方位角，即可求得该面元的朝向。(2)输入太阳入射位置和能量分布状况等数据，计算每个地表网格面元接收到的 初始辐射通量密度。假设太阳入射能量由直射光构成，天顶角(e s)为63°，方位角(<ps)为260°，其 直射辐射通量密度rs为1瓦/平方米。每个地表网格面元接收到的初始辐射通量密度计 算按照Fs(i) = [rs/cos( e s) ] sd|a(i，d) = 1/cos (63*pi/180) ^ sd|a(i, d)其中，ni为面元i的法线单位向量，sd为太阳方向单位向量(sinQs'Cos(ps, sinQs -sin(ps, cosQs ) , a(i，d)为该面元可以直接被太阳照射到的面积比例，此值通过 将所有面元向太阳方向投影，然后进行深度排序计算得到。(3)输入地表网格面元的反射率数据，使用计算机图形学的方法计算地表网格面 元之间的可视因子。以计算面元卓j时)之间的可视因子为例，具体计算过程如 下a.在面元&上微面元dAi中心处沿其法向量方向放置一个虚拟的半立方体，该半 立方体的五个表面被剖分成均勻正方形网格(通常剖分为100X100)，每一网格均对应dAi 朝向半球空间的一个微立体角，从而形成一个半空间立方体角查找表。预先计算，并存贮半 立方体中心微面元dAi对各表面网格的微可视因子，设q为半立方体表面上一个网格，微面 元dAi对q的微可视因子为Fq ；b.以dAi为中心将面元A」投影到半立方体的五个表面上，计算投影区域包含的半立方体表面网格，对表面网格对应的微可视因子进行叠加，计算出dAi对A」的可视因子，设Q为投影区域包含的半立方体表面网格集合，则dAi对A」的可视因子可以表示为 c.上述步骤(a和b)计算中没有考虑面元之间的可视问题，对有遮挡的情形，若半 立方体表面网格被两个或者两个以上面元的投影区域所覆盖，则通过比较这些面元离半立 方体中心微面元的距离来决定在该像素处可见的面元，距离近的面元可见，距离远的面元 不可见；d.以a-c计算面元&上所有微面元到A」的可视因子，则 (4)由面元初始散射的辐射通量密度和可视因子，构建该复杂地形区域的辐射度 方程。假设地表面元i反射率为P”可得初始散射的辐射通量密度(E》 于是构建辐射度方程如下其中ByB^分别为面元间达到辐射能量平衡后的面元i和j的辐射通量密度，是待 求解量。以矩阵形式表示为
其中 (5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的 辐射通量密度。采用Gauss-Seidel迭代法求解辐射度方程，以初始散射的辐射通量密度(E)作 为初值B°，带入辐射度方程右边计算，得到一个新的值B1，然后再带入辐射度方程右边计 算，又可以得到一个新值B2，如此反复迭代，直到|Bn+1-Bn|≤10_6为止。迭代收敛后，可准 确计算出三维场景中每个面元的辐射通量密度Bp(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和辐亮度，实现 复杂地形区遥感过程的精确模拟和计算。假设传感器垂直观测，观测天顶角(e v)为0°，观测方位角(<PV)为0°，即可按 下式计算该观测方向上的反射率因子(DRF) 观测方向矢量( coy<)D1/，slnOv -Sincpv, coi^J ,Bi 为面元 i (以
方向矢量&表示，其面积为areaG))的辐射通量密度，a(i, v)为该面元对观测方向上可 见的面积比例(同样通过将面元向观测方向投影，然后进行深度排序计算得到)。
权利要求
一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征在于包含以下步骤(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备；(2)输入太阳入射位置和能量分布状况数据，计算每个地表网格面元接收到的初始太阳辐射通量密度；(3)使用计算机图形学的方法计算地表网格面元之间的可视因子；(4)输入地表网格面元的波谱特性数据，建立该复杂地形区域的辐射度方程；(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的辐射通量密度；(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和地表面元辐亮度，实现复杂地形区域辐射传输过程的精确模拟和计算。
2.根据权利要求1所述的一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征 在于所述步骤(2)中，地表网格面元，以下简称面元，接收到的初始辐射通量密度(Fs)由 太阳入射角度和能量和面元朝向决定，具体计算按照(以面元i为例)Fs ⑴=[rs/cos ( Q s) ] | 叫 sd | a (i, d)其中，rs为太阳入射辐射通量密度，e s为太阳天顶角，叫为面元i的法线单位向量，sd 为太阳方向单位向量，a(i，d)为该面元可以直接被太阳照射到的面积比例。
3.根据权利要求1所述的一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征 在于所述步骤(3)中可视因子定义为每两个面元之间的可视面积比例，用以描述面元间 能量交换的比例，面元i与面元j之间的可视因子(Fu)表示为 其中，9,和^分别为面元i与面元j各自法线与两面元连线的夹角，r为两面元之 间的距离，&为面元i的面积，^为面元j的面积。
4.根据权利要求1所述的一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征 在于所述步骤(4)中建立辐射度方程，即面元间达到辐射平衡后的能量方程 B” Bj分别为辐射平衡后面元i和j的辐射通量密度，E,为面元i的初始辐射出射度， P i为面元反射率。Ei由面元接收到的初始太阳辐射通量密度Fs和面元反射率P i决定E(i) = Fs(i) P j = [rs/cos( 0 s) ] sd|a(i，d) P 4对场景中所有面元，共N个，建立辐射度方程，构成线性方程组，以矩阵形式表示为B = E+CB其中，B和E是1XN维向量，C是NXN维矩阵
5.根据权利要求1所述的一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征在于所述步骤(5)中采用Gauss-Seidel迭代法求解辐射度方程，得到场景中每个面元的辐射通量密度。
6.根据权利要求1所述的一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征 在于所述步骤(6)中，计算该复杂地形却与的方向反射率因子是根据指定的传感器观测 方向，将所有该观测方向上能够观测到的面元的辐射通量累加，得到该方向上场景总的辐 射通量，该值与相同条件下理想漫反射面的辐射通量的比值，即为该场景的方向反射率因 子。
全文摘要
本发明涉及一种基于辐射度的复杂地形区域辐射传输模拟方法，步骤如下(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备；(2)输入太阳入射位置和能量分布状况数据，计算每个地表网格面元接收到的初始太阳辐射通量密度；(3)使用计算机图形学的方法计算地表网格面元之间的可视因子；(4)输入地表网格面元的波谱特性数据，建立该复杂地形区域的辐射度方程；(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的辐射通量密度；(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和地表面元辐亮度，实现复杂地形区域辐射传输过程的精确模拟和计算。
文档编号G01S7/497GK101876700SQ200910243148
公开日2010年11月3日 申请日期2009年12月29日 优先权日2009年12月29日
发明者王亚超, 董艳芳, 贾国瑞, 赵峰, 赵慧洁 申请人:北京航空航天大学