CN101034475A - 无阴影卫星遥感正射数字图像的计算机生成方法 - Google Patents

Abstract

本发明涉及一种无阴影卫星遥感正射数字图像的计算机生成方法。本方法是：首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像：生成太阳直射光遥感图像和天空散射光遥感图像、计算像元太阳直射光变换系数，天空散射光变换系数和像元地面辐射参数，即散射/直射比，散射/辐射比；建立太阳直射光和天空散射光遥感数字图像的水平地形变换模型，分别进行太阳直射光和天空散射光遥感数字图像的水平地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像无阴影正射影像的计算机生成。本发明的理论、方法较传统遥感数字处理先进、科学、合理，可广泛应用于遥感技术的定量研究和应用。

Description

无阴影卫星遥感正射数字图像的计算机生成方法

技术领域：

本发明涉及一种卫星遥感数字图像处理方法，特别是一种遥感数字图像上阴影消除和阴影中像元遥感值恢复的原理方法。

背景技术：

地形对遥感信息的影响与改正(包括山体阴影、云层阴影的消除)一直都是卫星遥感数字图像处理与应用中的难题。现有卫星遥感数字图像处理与应用技术与正射影像处理技术，均没有从根本上解决这些重要的问题。经过查新搜索国内外均没有相关文献和专利报道，本项发明从理论和实践上较好地解决了这一问题。

发明内容：

本发明的目的在于针对现有技术存在的缺陷，提供一种无阴影卫星遥感正射数字图像的计算机生成方法。

为达到上述发明目的，本发明的构思是：

本项发明是在GIS支持下建立和利用自然中太阳直射光、天空散射光照度、地形与卫星遥感数字图像相互之间的定量关系，分解原始卫星遥感数字图像为太阳直射光和天空散射光遥感数字图像。在此基础上利用数字地形图计算像元太阳直射光变换系数，天空散射光变换系数和像元地面辐射参数(散射/直射比，散射/辐射比)进行太阳直射光和天空散射光遥感数字图像的水平地形变换(如将起伏地形上的卫星遥感数字图像转换为水平地面上相应卫星遥感数字图像)。

同时进行地形阴影和云层阴影的判断与消除，最后完成原始卫星遥感数字图像的无阴影正射影像的计算机生成。

卫星遥感数字图像无阴影正射影像的计算机生成必须使生成的该图像的每个像元遥感值均严格定量满足卫星遥感图像的点成像原理。

根据上述的发明构思，本发明采用下述技术方案：

一种无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于：首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像：生成太阳直射光遥感图像和天空散射光遥感图像、计算像元太阳直射光变换系数，天空散射光变换系数和像元地面辐射参数，即散射/直射比，散射/辐射比；建立太阳直射光和天空散射光遥感数字图像的水平地形变换模型，分别进行太阳直射光和天空散射光遥感数字图像的水平地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像无阴影正射影像的计算机生成。

其具体操作步骤如下：

(1)、卫星遥感数字图像与数字地形图配准。

(2)、卫星遥感数字图像的大气修正。

(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算。

(4)、卫星遥感数字图像上山体与云层阴影的判断。

(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比。

(6)、计算卫星遥感数字图像上各点地面太阳直射光地形影响系数。

(7)、计算卫星遥感数字图像上各点地面天空散射光地形影响系数。

(8)、解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像。

(9)、太阳直射光遥感数字图像和天空散射光遥感数字图像的水平地形变换。

(10)、卫星遥感数字图像上山体与云层阴影的消除。

(11)、卫星遥感数字图像无阴影正射影像的计算机生成。

上述的骤(1)中的卫星遥感数字图像与数字地形图配准是GIS常规操作。

上述的步骤(2)中的卫星遥感数字图像的大气修正的步步骤如下：

利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，按下式计算可以得到离散分布的像元大气程辐射遥感值。

                 DN(a)＝γ1*(DN1-DN2)/(γ1-γ2)                  (1)

或者对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析也可以得到离散分布的像元大气程辐射遥感值；

对上述结果，再采用插值算法得到每个像元大气程辐射遥感值；

对卫星遥感数字图像上各像元点遥感值DN_ij进行减去大气程辐射遥感值DN_Aij的运算，即DN_ij-DN_Aij，便完成了卫星遥感数字图像的大气修正。

其中，DN_ij代表卫星遥感数字图像上各像元点遥感值，DN_Aij代表大气程辐射遥感值，DN_ij-DN_Aij代表卫星遥感数字图像的大气修正值。

上述的步骤(3)中的卫星遥感数字图像上各点太阳高度角、方位角的计算按下式进行计算：

θ_ij＝arcsin(sinφ*sinδ+cosφ*cosδ*cost(i，j))                   (2)

A(i，j)＝arcsin(sinθ_ij*sinφ-sinδ)/cosθ_ij*cosφ)                (3)

δ＝arcsin(sinθ*sinφ-cosθ*cosφ*cosA)                           (4)

t＝arcsin(cosθ_ij*sinA/cosδ)+Δλ                                 (5)

θ，A、δ：卫片注记中星下点太阳高度角、方位角和太阳赤角；λ、φ：分别为星下点地理经纬度。Δλ为像元点对星下点的经度增量。

上述的步骤(4)中的卫星遥感数字图像上山体与与云层阴影的判断方法如下：

利用与卫星遥感数字图像配准的数字地形图和卫片像元的太阳的位置参数，即高度角、方位角，进行山体与云层阴影的判断，其判断的准则是：在太阳照射方向该像元点的最大地形、云层、高度角等于或者大于该像元点的太阳高度角即：DH(i，j)≥θ_ij，则该像元点为阴影；反之则不是阴影。

上述的步骤(5)中的观察测量或者计算卫星遥感数字图像上各点水平地面上的直、散射辐射比方法如下：

若卫片水平地面上有相邻2个同类像元X和Y，像元Y位于阴影中，DN(y)和DN(x)、DN_a分别是其遥感值和大气程辐射辐射遥感值，则该处水平地面上的散、直射辐射比：

      L＝((DN(y)-DN_a)/(DN(x)-DN(y))                                (6)

上述的步骤(6)中的计算卫星遥感数字图像上各点地面直射光地形影响系数按下式计算：

      F_ij＝1-tgα_ij·ctgθ_ij·cosω_ij，ω_ij＝AL_ij-A_ij              (7)

其中，α_ij代表像元地表坡角，AL_ij代表地表坡向，θ_ij代表卫片注记中星下点太阳高度角，A_ij代表太阳方位角，ω_ij代表太阳方位角与地面坡向之夹角。

上述的步骤(7)中的计算卫星遥感数字图像上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比按下式计算：

$G_{\mathrm{ij}}=1-\frac{2}{\mathrm{n\π}}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\β}}_k,$ 式中n＝2π/Δt                              (8)

其中，β_k为第K个方位上的最大平均遮蔽角即立体角。Δt为方位角步长，为360°/n。n为方位划分数，即n为按计算精度要求将360°划分的范围个数。

上述的步骤(8)中的解析卫星遥感图像；生成太阳直射光遥感图像和天空散射光遥感图像，方法如下：

太阳直射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Sij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·F_{\mathrm{ij}}}{(F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}})}---\left(9\right)$

天空散射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Dij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·G_{\mathrm{ij}}\·L_{\mathrm{ij}}}{F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}}}---\left(10\right)$

i＝1，2…m；j＝1，2…n

其中，DN_Sij代表像元遥感数值中直射光分量，DN_ij代表像元遥感数值，DN_Aij代表大气程辐射遥感值，L_ij代表地面反射亮度，F_ij代表数字地形图上各点地面直射光地形影响系数，G_ij代表数字地形图上各点地面天空散射光地形影响系数，DN_Dij代表像元遥感数值中散射光分量，L_Sij代表地面光谱反射亮度，L_Dij代表地面辐射照度。

上述的步骤(9)中的太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算：

        DN_Sij′＝DN_Sij/F_ij＝(DN_ij-DN_Aij)/(F_ij+G_ij·L_ij)                    (11)

        DN_Dij′＝DN_Dij/G_ij＝(DN_ij-DN_Aij)·L_ij/(F_ij+G_ij·L_ij)               (12)

i＝1，2…m；j＝1，2…n

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Sij′代表经大气修正后的卫星像元遥感数据中直射光分量，DN_Dij′代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

上述的步骤(10)中的卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数字运算：

     DN_ij′＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)               (13)

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Aij代表大气程辐射遥感值。

上述的步骤(11)中的卫星遥感数字图像无阴影正射影像的计算机生成，使用以下数字运算完成：

      DN_ij′＝DN_Sij′+DN_Dij′                                  (14)

      DN_ij′＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)             (15)

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Sij′代表经大气修正后的卫星像元遥感数据中直射光分量，DN_Dij′代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

上述的像元大气程辐射辐射D_Aij，通过数字图像中相关信息中提取或地面观测数据计算得出；上述的像元地面太阳高度θ_ij、方位角AL_ij由卫片注记给出的相关信息计算得出；上述的像元地面坡度α_ij、坡向A_ij以及像元经、纬度由卫片对应DTM给出。

本发明与现有技术相比，具有如下显而易见的突出实质性特点和显著优点：

由于卫星遥感成像是在自然光条件之下，自然光成份包括太阳直射光(方向性)分量和天空散射光(基本上是各向同性)，其相对强度以及在起伏地表上的再分配方式，对卫星遥感数据的影响有很大差异。将大气修正后卫星遥感图像解析，分离成太阳直射光遥感图像和天空散射光遥感图像，分别进行太阳直射光和天空散射光的变换，形成对应的太阳直射光和天空散射光的地形变换遥感数字图像，最后将二者合成自然光条件下的地形变换卫星遥感数字图像。其理论、方法较传统遥感数字图像处理(自然光为整体)先进、科学、合理。同时，所生成的：太阳直射光和天空散射光遥感数字图像以及它们的地形变换数字图像都是自然界无法直接获取的，在理论上和实践上均有重要价值的遥感信息。

由于遥感数字图像的地形变换消除了地形起伏带来的地表光辐射差异以及云层与山体阴影的影响，在遥感信息中突出了地物波谱特征，使其更具可比性。从根本上改善了遥感数字图像的质量具有重要的科研和应用前景。

本发明可广泛应用于遥感技术的定量研究与应用，如资源环境的遥感调查研究与成像、军事中的伪装与反伪装、阴影地物的识别与遥感数字图像的模式识别、地物波谱特性的遥感反演、遥感数字图像的虚拟与仿真等方面。

附图说明：

图1为本发明一个实施例的无阴影卫星遥感正射数字图像计算机生成流程图。

图2为福建龙岩地区1998.12.8经过大气修正后的TM遥感数字图像(1024×1024)；

图3为与图2卫片匹配的对应数字地形图(立体)；

图4为图2卫片解析后生成的太阳直射光遥感数字图像。

图5为图2卫星解析后生成的天空散射光遥感数字图像。

图6为太阳直射光地形变换系数的线性拉伸图。

图7为天空散射光地形变换系数的线性拉伸图。

图8为太阳直射光地形修正系数所成的图像。

图9为天空散射光地形修正系数所成的图像。

图10为计算机生成的无阴影卫星遥感数字正射影像。

具体实施方式

本发明的一个优选实施例结合附图详细说明如下：本例为福建龙岩地区1998.12.8的卫片以计算机生成无阴影卫星遥感正射数字图像。

参见图1，本无阴影卫星正射数字图像的计算机生成方法是：首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像：生成太阳直射光遥感图像和天空散射光遥感图像、计算像元太阳直射光变换系数，天空散射光变换系数和像元地面辐射参数，即散射/直射比，散射/辐射比；建立太阳直射光和天空散射光遥感数字图像的水平地形变换模型，分别进行太阳直射光和天空散射光遥感数字图像的水平地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像无阴影正射影像的计算机生成。

其具体操作步骤如下：

(1)、卫星遥感数字图像与数字地形图配准。

(2)、卫星遥感数字图像的大气修正。

(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算。

(4)、卫星遥感数字图像上山体与云层阴影的判断。

(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比。

(6)、计算卫星遥感数字图像上各点地面太阳直射光地形影响系数。

(7)、计算卫星遥感数字图像上各点地面天空散射光地形影响系数。

(8)、解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像。

(9)、太阳直射光遥感数字图像和天空散射光遥感数字图像的水平地形变换。

(10)、卫星遥感数字图像上山体与云层阴影的消除。

(11)、卫星遥感数字图像无阴影正射影像的计算机生成。

卫星遥感数字图像的地形变换在有多点同步实时辐射观测资料(散射/直射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。

上述方案中各步骤说明如下：

1、卫星遥感数字图像与数字地形图配准：卫星遥感数字图像与数字地形图(TEM)配准系GIS常规操作。

2、星遥感数字图像的大气修正：获得像元大气程辐射遥感值的方法有两种：

利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，进行以下计算：DN(a)＝γ1*(DN1-DN2)/(γ1-γ2)                            (1)

可以得到离散分布的像元大气程辐射遥感值。

对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析也可以得到离散分布的像元大气程辐射遥感值。

对上述结果，再采用插值算法可得到每个像元大气程辐射遥感值。对卫星遥感数字图像上各像元点遥感值进行减去大气程辐射遥感值的运算，即DN_ij-DN_Aij，便完成了卫星遥感数字图像的大气修正。

其中，DN_ij代表卫星遥感数字图像上各像元点遥感值，DN_Aij代表大气程辐射遥感值，DN_ij-DN_Aij代表卫星遥感数字图像的大气修正值。

3、卫星遥感数字图像上各点太阳高度角、方位角的计算按下式进行计算：

θ_ij＝arcsin(sinφ*sinδ+cosφ*cosδ*cost(i，j))                      (2)

A(i，j)＝arcsin(sinθ_ij*sinφ-sinδ)/cosθ_ij*cosφ)                   (3)

δ＝arcsin(sinθ*sinφ-cosθ*cosφ*cosA)                              (4)

t＝arcsin(cosθ_ij*sinA/cosδ)+Δλ                                    (5)

θ，A、δ：卫片注记中星下点太阳高度角、方位角和太阳赤角；λ、φ：分别为星下点地理经纬度。Δλ为像元点对星下点的经度增量。

4、卫星遥感数字图像上山体与云层阴影的判断

利用与卫星遥感数字图像配准的数字地形图和卫片像元的太阳的位置参数，即高度角、方位角，进行山体与云层阴影的判断，其判断的准则是：在太阳照射方向该像元点的最大地形、云层、高度角等于或者大于该像元点的太阳高度角即：DH(i，j)≥θ_ij，则该像元点为阴影；反之则不是阴影。

5、观测或计算卫星遥感数字图像上各点水平地面上的直、散射辐射比：若卫片水平地面上有相邻2个同类像元X和Y，像元Y位于阴影中，DN(y)和DN(x)、DNa分别是其遥感值和大气程辐射辐射遥感值，则该处水平地面上的散、直射辐射比：

     L＝((DN(y)-DN_a)/(DN(x)-DN(y))                                    (6)

6、计算数字地形图上各点地面直射光地形影响系数：

           F_ij＝1-tgα_ij·ctgθ_ij·cosω_ij，ω_ij＝AL_ij-A_ij            (7)

其中，α_ij代表像元地表坡角，AL_ij代表地表坡向，θ_ij代表卫片注记中星下点太阳高度角，A_ij代表太阳方位角，ω_ij代表太阳方位角与地面坡向之夹角。

7、计算数字地形图上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比：

$G_{\mathrm{ij}}=1-\frac{2}{\mathrm{n\π}}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\β}}_k,$ 式中n＝2π/Δt                     (8)

其中，β_k为第K个方位上的最大平均遮蔽角即立体角。Δt为方位角步长，为360°/n。n为方位划分数，即n为按计算精度要求将360°划分的范围个数。

8、解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像，方法如下：

太阳直射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Sij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·F_{\mathrm{ij}}}{(F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}})}---\left(9\right)$

天空散射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Dij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·G_{\mathrm{ij}}\·L_{\mathrm{ij}}}{F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}}}---\left(10\right)$ i＝1，2…m；j＝1，2…n

其中，DN_Sij代表像元遥感数值中直射光分量，DN_ij代表像元遥感数值，DN_Aij代表大气程辐射遥感值，L_ij代表地面反射亮度，F_ij代表数字地形图上各点地面直射光地形影响系数，G_ij代表数字地形图上各点地面天空散射光地形影响系数，DN_Dij代表像元遥感数值中散射光分量，L_Sij代表地面光谱反射亮度，L_Dij代表地面辐射照度。

9、太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算：

       DN_Sij′＝DN_Sij/F_ij＝(DN_ij-DN_Aij)/(F_ij+G_ij·L_ij)                (11)

       DN_Dij′＝DN_Dij/G_ij＝(DN_ij-DN_Aij)·L_ij/(F_ij+G_ij·L_ij)           (12)

i＝1，2…m；j＝1，2…n

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Sij′代表经大气修正后的卫星像元遥感数据中直射光分量，DN_Dij′代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

10、卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数学运算：

            DN_ij′＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)              (13)

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Aij代表大气程辐射遥感值。

11、卫星遥感数字图像无阴影正射影像的计算机生成，使用以下数字运算完成：

     DN_ij′＝DN_Sij′+DN_Dij′                                        (14)

     DN_ij′＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)                   (15)

其中，DN_ij代表像元遥感数值，DN_ij′代表经大气修正后的卫星像元遥感数据，DN_Sij′代表经大气修正后的卫星像元遥感数据中直射光分量，DN_Dij′代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

上述的像元大气程辐射辐射DN_Aij，通过数字图像中相关信息中提取或地面观测数据计算得出；上述的像元地面太阳高度θ_ij、方位角AL_ij由卫片注记给出的相关信息计算得出；上述的像元地面坡度α_ij、坡向A_ij以及像元经、纬度由卫片对应DTM给出。卫星遥感数字图像的地形变换在有多点实时辐射观测资料(散射/辐射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。

由于卫星遥感成像是在自然光条件之下，自然光成份包括太阳直射光(方向性)和天空散光(基本上是各向同性)分量，其相对强度以及在起伏地表上的再分配方式以及对卫星遥感数据的影响都有很大差异。将大气修正后卫星遥感图像解析，分离成太阳直射光遥感图像和天空散射光遥感图像，分别进行太阳直射光和天空散射光的变换以形成对应的太阳直射光和天空散射光的地形变换遥感数字图像，最后将二者合成自然光条件下的地形变换卫星遥感数字图像。其理论、方法较传统遥感数字图像处理(自然光为整体)更为先进、科学、合理。同时，所生成的：太阳直射光和天空散射光遥感数字图像以及它们的地形变换数字图像都是自然界无法直接获取的，在理论上和实践上均有重要价值的遥感信息。

遥感数字图像的地形变换消除了地形起伏带来的地表光辐射差异以及云层与山体阴影的影响，在遥感信息中突出了地物波谱特征，使其更具可比性。从根本上改善了遥感数字图像的质量具有重要的科研和应用前景，与现有技术相比具有的显著的优点与积极的效果。

图2为福建龙岩地区1998.12.8经过大气修正后的TM遥感数字图像(1024×1024)；

图3为与该卫片匹配的对应数字地形图(立体)；

图4、图5分别为该卫片解析后生成的太阳直射光和天空散射光遥感数字图像：太阳直射光遥感数字图像(图4)为地形影响的太阳直射光下的遥感数字图像上有明显的阴影分布、山脊两边直射光辐射有明显差异；

天空散射光遥感数字图像(图5)为地形影响的低照度天空散射光下的遥感数字图像。

图6、图7分别地形直射光地形变换系数和天空散射光地形变换系数的可视化，它们有着十分清晰的物理和地学意义：

图6：为地形直射光地形变换系数的线性拉伸。表达成像瞬间太阳直射光辐射在起伏地面上的真实归一化(水平地面上太阳直射光辐射照度为1)分布。

图7：为天空散射光地形变换系数的线性拉伸。表达成像瞬间天空散射光辐射在起伏地面上的真实归一化(水平地面上天空散射光辐射照度为1)分布。直射光地形变换系数和天空散射光地形变换系数及其可视化在数字沙盘、地图晕染和数字可视化等方面有广泛应用。

图8、图9分别为太阳直射光和天空散射光遥感数字图像所对应的水平地形变换的图像：

太阳直射光(非偏振光)遥感数字图像的地形变换图像(图8)上阴影已经消除，原来隐藏在阴影中的水体清晰可见，山脊两边直射光辐射差异已不明显，起伏山地已变成平原。与原来的遥感数字图像相比太阳直射光遥感数字图像的地形变换图像显然有利于计算机的自动识别与分类，其图像质量有重大的改善和提高。

天空散射光遥感数字图像的地形变换图像(图9)为消除了地形影响在各向均匀的天空散射光(偏振光)条件下的成像的遥感数字图像。

图10为计算机生成的无阴影卫星遥感数字正射影像(图8和图9的合成)。

卫星遥感数字图像的地形变换在有多点同步实时辐射观测资料(散射/直射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。

Claims (12)

1.一种无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于：

首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像：生成太阳直射光遥感图像和天空散射光遥感图像、计算像元太阳直射光变换系数，天空散射光变换系数和像元地面辐射参数，即散射/直射比，散射/辐射比；建立太阳直射光和天空散射光遥感数字图像的水平地形变换模型，分别进行太阳直射光和天空散射光遥感数字图像的水平地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像无阴影正射影像的计算机生成；

其具体操作步骤如下：

(1)、卫星遥感数字图像与数字地形图配准，

(2)、卫星遥感数字图像的大气修正，

(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算，

(4)、卫星遥感数字图像上山体与云层阴影的判断，

(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比，

(6)、计算卫星遥感数字图像上各点地面太阳直射光地形影响系数，

(7)、计算卫星遥感数字图像上各点地面天空散射光地形影响系数，

(8)、解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像，

(9)、太阳直射光遥感数字图像和天空散射光遥感数字图像的水平地形变换，

(10)、卫星遥感数字图像上山体与云层阴影的消除，

(11)、卫星遥感数字图像无阴影正射影像的计算机生成。

2.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(2)中的卫星遥感数字图像的大气修正的步骤如下：

(1)利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，

按下式计算可以得到离散分布的像元大气程辐射遥感值：

               DN(a)＝γ1*(DN1-DN2)/(γ1-γ2)；

或者对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析得到离散分布的像元大气程辐射遥感值；

(2)对上述结果，再采用插值算法得到每个像元大气程辐射遥感值；

(3)对卫星遥感数字图像上各像元点遥感值DN_ij进行减去大气程辐射遥感值DN_Aij的运算，即DN_ij-DN_Aij，便完成了卫星遥感数字图像的大气修正，

其中，DN_ij代表卫星遥感数字图像上各像元点遥感值，DN_Aij代表大气程辐射遥感值，DN_ij-DN_Aij代表卫星遥感数字图像的大气修正值。

3.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(3)中的卫星遥感数字图像上各点太阳高度角、方位角的计算按下式进行计算：

      θ_ij＝arcsin(sinφ*sinδ+cosφ*cosδ*cost(i，j))，

      A(i，j)＝arcsin(sinθ_ij*sinφ-sinδ)/cosθ_ij*cosφ)，

           δ＝arcsin(sinθ*sinφ-cosθ*cosφ*cosA)，

            t＝arcsin(cosθ_ij*sinA/cosδ)+Δλ，

θ，A、δ：卫片注记中星下点太阳高度角、方位角和太阳赤角；λ、φ：分别为星下点地理经纬度，Δλ为像元点对星下点的经度增量。

4.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(4)中的卫星遥感数字图像上山体与云层阴影的判断方法如下：

利用与卫星遥感数字图像配准的数字地形图和卫片像元的太阳的位置参数，即高度角、方位角，进行山体与云层阴影的判断，其判断的准则是：在太阳照射方向该像元点的最大地形、云层、高度角等于或者大于该像元点的太阳高度角即：DH(i，j)≥θ_ij，则该像元点为阴影；反之则不是阴影。

5.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(5)中的观测或计算卫星遥感数字图像上各点水平地面上的直、散射辐射比方法如下：

若卫片水平地面上有相邻2个同类像元X和Y，像元Y位于阴影中，DN(y)和DN(x)、DN_a分别是其遥感值和大气程辐射遥感值，则该处水平地面上的散、直射辐射比：

              L＝(DN(y)-DN_a)/(DN(x)-DN(y))。

6.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(6)中的计算数字地形图上各点地面直射光地形影响系数按下式计算：

        F_ij＝1-tgα_ij·ctgθ_ij·cosω_ij，ω_ij＝AL_ij-A_ij，

其中，α_ij代表像元地表坡角，AL_ij代表地表坡向，θ_ij代表卫片注记中星下点太阳高度角，A_ij代表太阳方位角，ω_ij代表太阳方位角与地面坡向之夹角。

7.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(7)中的计算数字地形图上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比按下式计算：

$G_{\mathrm{ij}}=1-\frac{2}{\mathrm{n\π}}\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\β}}_k,$ 式中n＝2π/Δt，

其中，β_k为第K个方位上的最大平均遮蔽角即立体角，Δt为方位角步长，为360°/n，n为方位划分数，即n为按计算精度要求将360°划分的范围个数。

8.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(8)中的解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像，方法如下：

          太阳直射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Sij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·F_{\mathrm{ij}}}{(F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}})},$

          天空散射光遥感图像表达式： ${\mathrm{DN}}_{\mathrm{Dij}}=\frac{({\mathrm{DN}}_{\mathrm{ij}}-{\mathrm{DN}}_{\mathrm{Aij}})\·G_{\mathrm{ij}}\·L_{\mathrm{ij}}}{F_{\mathrm{ij}}+G_{\mathrm{ij}}\·L_{\mathrm{ij}}},$

          i＝1，2…m，j＝1，2…n，

其中，DN_Sij代表像元遥感数值中直射光分量，DN_ij代表像元遥感数值，DN_Aij代表大气程辐射遥感值，L_ij代表地面反射亮度，F_ij代表数字地形图上各点地面直射光地形影响系数，G_ij代表数字地形图上各点地面天空散射光地形影响系数，DN_Dij代表像元遥感数值中散射光分量，L_Sij代表地面光谱反射亮度，L_Dij代表地面辐射照度。

9.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(9)中的太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算：

           DN′_Sij＝DN_Sij/F_ij＝(DN_ij-DN_Aij)/(F_ij+G_ij·L_ij)，

           DN′_Dij＝DN_Dij/G_ij＝(DN_ij-DN_Aij)·L_ij/(F_ij+G_ij·L_ij)，

           i＝1，2…m，j＝1，2…n，

其中，DN_ij代表像元遥感数值，DN′_ij代表经大气修正后的卫星像元遥感数据，DN′_Sij代表经大气修正后的卫星像元遥感数据中直射光分量，DN′_Dij代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

10.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(10)中的卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数字运算：

            DN′_ij＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)，

其中，DN_ij代表像元遥感数值，DN′_ij代表经大气修正后的卫星像元遥感数据，DN_Aij代表大气程辐射遥感值。

11.根据权利要求1所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的步骤(11)中的卫星遥感数字图像无阴影正射影像的计算机生成，使用以下数字运算完成：

            DN′_ij＝DN′_Sij+DN′_Dij，

            DN′_ij＝(DN_ij-DN_Aij)·(1+L_ij)/(F_ij+G_ij·L_ij)，

其中，DN_ij代表像元遥感数值，DN′_ij代表经大气修正后的卫星像元遥感数据，DN′_Sij代表经大气修正后的卫星像元遥感数据中直射光分量，DN′_Dij代表经大气修正后的卫星像元遥感数据中散射光分量，DN_Dij代表像元遥感数值中散射光分量，DN_Sij代表像元遥感数值中直射光分量，DN_Aij代表大气程辐射遥感值。

12.根据权利要求2或6所述的无阴影卫星遥感正射数字图像的计算机生成方法，其特征在于所述的像元大气程辐射辐射DN_Aij，通过数字图像中相关信息中提取或地面观测数据计算得出；所述的像元地面太阳高度θ_ij、方位角AL_ij由卫片注记给出的相关信息计算得出；所述的像元地面坡度α_ij、坡向A_ij以及像元经、纬度由卫片对应DTM给出。

Images