CN101876700B - 一种基于辐射度的复杂地形区域辐射传输模拟方法 - Google Patents

Abstract

本发明涉及一种基于辐射度的复杂地形区域辐射传输模拟方法，步骤如下：(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备；(2)输入太阳入射位置和能量分布状况数据，计算每个地表网格面元接收到的初始太阳辐射通量密度；(3)使用计算机图形学的方法计算地表网格面元之间的可视因子；(4)输入地表网格面元的波谱特性数据，建立该复杂地形区域的辐射度方程；(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的辐射通量密度；(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和地表面元辐亮度，实现复杂地形区域辐射传输过程的精确模拟和计算。

Description

一种基于辐射度的复杂地形区域辐射传输模拟方法

技术领域

本发明涉及一种基于辐射度的复杂地形区域辐射传输模拟方法，适合于对山区复杂地形区域条件下的遥感成像过程精确模拟，特别是定量分析地形对遥感影像的影响和影像地形校正方面具有重要价值，属于遥感成像模拟技术领域。

背景技术

在地形起伏情况下，地物将接收到周围物体反射的太阳直射光和大气散射光，导致地物实际接收辐射能量的增加。在雪地、植被近红外波段以及陡峭坡面等情况下，周围环境对目标物体反射的影响无法忽略。地物在接收周围物体反射能量的同时，也将部分能量重新反射回周围物体，循环往复，形成复杂的地表多次散射过程，导致地物光谱信息的变化，影响了山区遥感图像的信息提取精度，成为了限制山区遥感图像应用水平的主要障碍。地表多次散射是地表太阳能计算以及遥感仿真中的一个难点问题。

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 effects in 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)中，地表网格面元(以下简称面元)接收到的初始太阳辐射通量密度(F_s)由如下公式计算：

F_s(i)＝[r_s/cos(θ_s)]|n_i·s_d|a(i，d)

其中，r_s为太阳入射辐射通量密度，θ_s为太阳天顶角，n_i为面元i的法线单位向量，s_d为太阳方向单位向量，a(i，d)为该面元可以直接被太阳照射到的面积比例，F_s为该面元接收到的初始辐射通量密度。

其中，步骤(3)中利用计算机图形学中的半立方体算法计算面元之间的可视因子，以计算场景中任意两个面元A_i和A_j(i≠j时)之间的可视因子为例，具体计算过程如下：

a.在面元A_i上微面元dA_i中心处沿其法向量方向放置一个虚拟的半立方体，该半立方体的五个表面被剖分成均匀正方形网格(通常剖分为100×100)，每一网格均对应dA_i朝向半球空间的一个微立体角，从而形成一个半空间立方体角查找表。预先计算，并存贮半立方体中心微面元dA_i对各表面网格的微可视因子，设q为半立方体表面上一个网格，微面元dA_i对q的微可视因子为F_q；

b.以dA_i为中心将面元A_j投影到半立方体的五个表面上，计算投影区域包含的半立方体表面网格，对表面网格对应的微可视因子进行叠加，计算出dA_i对A_j的可视因子，设Q为投影区域包含的半立方体表面网格集合，则dA_i对A_j的可视因子可以表示为：

$F_{d{E}_i,E_j}=\underset{q\∈Q}{\mathrm{\Σ}}{F}_q$

c.上述步骤(a和b)计算中没有考虑面元之间的可视问题，对有遮挡的情形，若半立方体表面网格被两个或者两个以上面元的投影区域所覆盖，则通过比较这些面元离半立方体中心微面元的距离来决定在该像素处可见的面元，距离近的面元可见，距离远的面元不可见；

d.以a-c方法计算面元A_i上所有微面元dA_i到A_j的可视因子，则

$F_{\mathrm{ij}}=F_{A_i-A_j}=\frac{1}{A_i}{\∫}_{A_i}{\∫}_{A_j}\frac{\cos {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_{j}d{A}_{i}d{A}_j}{\mathrm{\π}{r}^2}$

在步骤(3)中，假设地表为朗伯反射，则面元反射率(ρ_i)乘以面元接收到的初始辐射通量密度(F_s)，可得该面元初始散射的辐射通量密度(E_i)：

E_i＝[r_s/cos(θ_s)]|n_i·s_d|a(i，d)·ρ_i

于是，可构造辐射度方程如下：

$B_i=E_i+{\mathrm{\ρ}}_i\underset{j=1,j\≠i}{\overset{N}{\mathrm{\Σ}}}{B}_j{F}_{\mathrm{ij}}$ (i＝1，...，N).

其中B_i，B_j分别为面元间达到辐射能量平衡后的面元i和j的辐射通量密度。以矩阵形式表示为：

B＝E+CB

其中，B和E是1×N维向量，C是N×N维矩阵：

$B=\left(\begin{array}{c}{B}_1\\ B_2\\ .\\ .\\ .\\ B_N\end{array}\right)$ $E=\left(\begin{array}{c}{E}_1\\ E_2\\ .\\ .\\ .\\ E_N\end{array}\right)$ $C=\left(\begin{array}{cccc}0& {\mathrm{\ρ}}_1{F}_{12}& ...& {\mathrm{\ρ}}_1{F}_{1N}\\ .& & & .\\ .& & & .\\ .& & & .\\ {\mathrm{\ρ}}_N{F}_{N1}& {\mathrm{\ρ}}_N{F}_{N2}& ...& 0\end{array}\right).$

步骤(5)中求解辐射度方程采用Gauss-Seidel迭代法，其基本过程为：首先给定一个初值B⁰，带入辐射度方程右边计算，得到一个新的值B¹，然后再带入辐射度方程右边计算，又可以得到一个新值B²，如此反复迭代，直到|Bⁿ⁺¹-Bⁿ|小于一个允许的误差范围为止。计算时，以初始散射的辐射通量密度(E_i)作为初值B⁰，误差范围一般选择|Bⁿ⁺¹-Bⁿ|≤10^-6即可满足精度要求。迭代收敛后，可准确计算出三维场景中每个面元的辐射通量密度。

步骤(6)中指定传感器观测方向，即可按下式计算该观测方向上的反射率因子(Directional Reflectance Factor，DRF)：

$\mathrm{DRF}\left(v\right)=\frac{\underset{i}{\mathrm{\Σ}}{B}_i/\mathrm{\π}|n_i\·s_v|a(i,v)\mathrm{area}\left(i\right)}{\underset{i}{\mathrm{\Σ}}1/\mathrm{\π}|n_i\·s_v|a(i,v)\mathrm{area}\left(i\right)}$

v为观测方向(以方向矢量s_v表示)，B_i为面元i(以方向矢量n_i表示，其面积为area(i))的辐射通量密度，a(i，v)为该面元对观测方向上可见的面积比例。

本发明是一种基于辐射度的复杂地形区域辐射传输模拟方法，与现有技术相比的优点在于：

(1)本发明基于严格的辐射传输过程，采用计算机图形学中经典的辐射度方法来对崎岖山区进行遥感反射率的精确模拟，纠正了现有技术中采用的各种近似和简化计算带来的误差。

(2)本发明对地形没有特殊的假定和要求，原则上任意复杂的地形区域，只要提供该区域的数字高程模型数据和地表反射率数据，即可精确计算传感器测量得到的反射率，大大扩展了山区辐射建模的应用范围，具有广泛的应用前景。

附图说明

图1为本发明一种基于辐射度的复杂地形区域辐射传输模拟方法的流程图。

具体实施方式

以西藏驱龙地区遥感场景模拟过程为例，如图1所示，本发明的具体实施方法如下：

(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备。

使用西藏驱龙地区地面分辨率为30m的数字高程模型和地表反射率数据。模拟计算时间设置为2002年4月3日上午10:6:56，区域中心地理坐标为29.98°N，36.37°E。根据数字高程模型数据提供的像元坐标信息，计算地表网格面元相邻两个边构成的向量的法向量，确定该法向量的天顶角和方位角，即可求得该面元的朝向。

(2)输入太阳入射位置和能量分布状况等数据，计算每个地表网格面元接收到的初始辐射通量密度。

假设太阳入射能量由直射光构成，天顶角(θ_s)为63°，方位角

为260°，其直射辐射通量密度r_s为1瓦/平方米。每个地表网格面元接收到的初始辐射通量密度计算按照：

F_s(i)＝[r_s/cos(θ_s)]|n_i·s_d|a(i，d)＝1/cos(63*pi/180)|n_i·s_d|a(i，d)

其中，n_i为面元i的法线单位向量，s_d为太阳方向单位向量

a(i，d)为该面元可以直接被太阳照射到的面积比例，此值通过将所有面元向太阳方向投影，然后进行深度排序计算得到。

(3)输入地表网格面元的反射率数据，使用计算机图形学的方法计算地表网格面元之间的可视因子。以计算面元A_i和A_j(i≠j时)之间的可视因子为例，具体计算过程如下：

a.在面元A_i上微面元dA_i中心处沿其法向量方向放置一个虚拟的半立方体，该半立方体的五个表面被剖分成均匀正方形网格(通常剖分为100×100)，每一网格均对应dA_i朝向半球空间的一个微立体角，从而形成一个半空间立方体角查找表。预先计算，并存贮半立方体中心微面元dA_i对各表面网格的微可视因子，设q为半立方体表面上一个网格，微面元dA_i对q的微可视因子为F_q；

b.以dA_i为中心将面元A_j投影到半立方体的五个表面上，计算投影区域包含的半立方体表面网格，对表面网格对应的微可视因子进行叠加，计算出dA_i对A_j的可视因子，设Q为投影区域包含的半立方体表面网格集合，则dA_i对A_j的可视因子可以表示为：

$F_{d{E}_i,E_j}=\underset{q\∈Q}{\mathrm{\Σ}}{F}_q$

c.上述步骤(a和b)计算中没有考虑面元之间的可视问题，对有遮挡的情形，若半立方体表面网格被两个或者两个以上面元的投影区域所覆盖，则通过比较这些面元离半立方体中心微面元的距离来决定在该像素处可见的面元，距离近的面元可见，距离远的面元不可见；

d.以a-c计算面元A_i上所有微面元dA_i到A_j的可视因子，则

$F_{\mathrm{ij}}=F_{A_i-A_j}=\frac{1}{A_i}{\∫}_{A_i}{\∫}_{A_j}\frac{\cos {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_{j}d{A}_{i}d{A}_j}{\mathrm{\π}{r}^2}$

(4)由面元初始散射的辐射通量密度和可视因子，构建该复杂地形区域的辐射度方程。假设地表面元i反射率为ρ_i，可得初始散射的辐射通量密度(E_i)：

E(i)＝F_s(i)·ρ_i＝[r_s/cos(θ_s)]|n_i·s_d|a(i，d)·ρ_i

于是构建辐射度方程如下：

$B_i=E_i+{\mathrm{\ρ}}_i\underset{j=1,j\≠i}{\overset{N}{\mathrm{\Σ}}}{B}_j{F}_{\mathrm{ij}}$

其中B_i，B_j分别为面元间达到辐射能量平衡后的面元i和j的辐射通量密度，是待求解量。以矩阵形式表示为：

B＝E+CB

其中：

$B=\left(\begin{array}{c}{B}_1\\ B_2\\ .\\ .\\ .\\ B_N\end{array}\right)$ $E=\left(\begin{array}{c}{E}_1\\ E_2\\ .\\ .\\ .\\ E_N\end{array}\right)$ $C=\left(\begin{array}{cccc}0& {\mathrm{\ρ}}_1{F}_{12}& ...& {\mathrm{\ρ}}_1{F}_{1N}\\ .& & & .\\ .& & & .\\ .& & & .\\ {\mathrm{\ρ}}_N{F}_{N1}& {\mathrm{\ρ}}_N{F}_{N2}& ...& 0\end{array}\right).$

(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的辐射通量密度。

采用Gauss-Seidel迭代法求解辐射度方程，以初始散射的辐射通量密度(E_i)作为初值B⁰，带入辐射度方程右边计算，得到一个新的值B¹，然后再带入辐射度方程右边计算，又可以得到一个新值B²，如此反复迭代，直到|Bⁿ⁺¹-Bⁿ|≤10^-6为止。迭代收敛后，可准确计算出三维场景中每个面元的辐射通量密度B_i。

(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和辐亮度，实现复杂地形区遥感过程的精确模拟和计算。

假设传感器垂直观测，观测天顶角(θ_v)为0°，观测方位角(

)为0°，即可按下式计算该观测方向上的反射率因子(DRF)：

$\mathrm{DRF}=\frac{\underset{i}{\mathrm{\Σ}}{B}_i/\mathrm{\π}|n_i\·s_v|a(i,v)\mathrm{area}\left(i\right)}{\underset{i}{\mathrm{\Σ}}1/\mathrm{\π}|n_i\·s_v|a(i,v)\mathrm{area}\left(i\right)}$

s_v为观测方向矢量

B_i为面元i(以方向矢量n_i表示，其面积为area(i))的辐射通量密度，a(i，v)为该面元对观测方向上可见的面积比例(同样通过将面元向观测方向投影，然后进行深度排序计算得到)。

Claims (1)

1.一种基于辐射度的复杂地形区域辐射传输模拟方法，其特征在于包含以下步骤： 

(1)输入研究区的数字高程模型，计算地表网格面元的朝向，进行数据准备； 

(2)输入太阳入射位置和能量分布状况数据，计算每个地表网格面元接收到的初始太阳辐射通量密度； 

(3)使用计算机图形学的方法计算地表网格面元之间的可视因子； 

(4)输入地表网格面元的波谱特性数据，建立该复杂地形区域的辐射度方程； 

(5)采用数值迭代的方法求解辐射度方程，得到辐射平衡后每个地表网格面元的辐射通量密度； 

(6)根据传感器采样位置计算该复杂地形区域的方向反射率因子和地表面元辐亮度，实现复杂地形区域辐射传输过程的精确模拟和计算； 

其中，所述步骤(2)中，地表网格面元，以下简称面元，接收到的初始辐射通量密度F_s由太阳入射角度和能量和面元i朝向决定，具体计算按照： 

F_s(i)=[r_s/cos(θ_s)]|n_i·s_d|a(i,d) 

其中，r_s为太阳入射辐射通量密度，θ_s为太阳天顶角，n_i为面元i的法线单位向量，s_d为太阳方向单位向量，a(i，d)为该面元直接被太阳照射到的面积比例； 

其中，所述步骤(3)中可视因子定义为每两个面元之间的可视面积比例，用以描述面元间能量交换的比例，面元i与面元j之间的可视因子F_ij表示为： 

其中，θ_i和θ_j分别为面元i与面元j各自法线与两面元连线的夹角，r为两面元之间的距离，A_i为面元i的面积，A_j为面元j的面积； 

其中，所述步骤(4)中建立辐射度方程，即面元间达到辐射平衡后的能量方程： 

B_i，B_j分别为辐射平衡后面元i和j的辐射通量密度，E_i为面元i的初始辐射出射度，ρ_i为面元反射率；E_i由面元接收到的初始太阳辐射通量密度F_s和面元反射率ρ_i决定： 

E(i)=F_s(i)·ρ_i=[r_s/cos(θ_s)]|n_i·s_d|a(i,d)·ρ_i

对场景中所有面元，共N个，建立辐射度方程，构成线性方程组，以矩阵形式表示为： 

B=E+CB 

其中，B和E是N×1维向量，C是N×N维矩阵： 

其中，所述步骤(5)中采用Gauss-Seidel迭代法求解辐射度方程，得到场景中每个面元的辐射通量密度； 

其中，所述步骤(6)中，计算该复杂地形区域的方向反射率因子是根据指定的传感器观测方向，将所有该观测方向上能够观测到的面元的辐射通量累加，得到该方向上场景总的辐射通量，该值与相同条件下理想漫反射面的辐射通量的比值，即为该场景的方向反射率因子。 

Images