专利名称：Gis支持下的卫星遥感数字图像的地形变换的制作方法
技术领域：
本发明涉及一种卫星遥感数字图像处理方法，特别是GIS支持下的卫星遥感数字图像的地形变换方法。
背景技术：
地形对卫星遥感数字图像的影响和山体阴影、云层阴影的消除，一直都是卫星遥感数字图像处理与应用中的难题。现有卫星遥感数字图像处理与应用技术与正射影像处理技术，均没有从根本上解决这个重要的问题。经过查新搜索国内外均没有相关文献和专利报道，本项发明在GIS支持下的卫星遥感数字图像的地形变换从理论和实践上较好地解决了这一问题。
发明内容
本发明的目的在于针对现有技术存在的缺陷，提供GIS支持下的卫星遥感数字图像的地形变换方法。
为达到上述发明目的，本发明的构思是本项发明是在GIS支持下利用自然中太阳直射光、天空散射光照度、地形与卫星遥感数字图像相互之间的定量关系，解析原始卫星遥感数字图像为太阳直射光和天空散射光遥感数字图像。在此基础上利用数字地形图计算像元太阳直射光变换(改正)系数，天空散射光变换(改正)系数和像元地面辐射参数(散射/直射比 散射/辐射比)进行太阳直射光和天空散射光遥感数字图像的水平地形变换(如将起伏地形上的卫星遥感数字图像转换为水平地面上相应卫星遥感数字图像)。
同时进行地形阴影和云层阴影的消除，完成GIS支持下的卫星遥感数字图像的地形变换和原始卫星遥感数字图像的任意地形变换，使每个像元遥感值均严格定量满足数字图像的点成像原理。
根据上述的发明构思，本发明采用下述技术方案GIS支持下的卫星遥感数字图像的地形变换方法，其特征在于首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像生成太阳直射光遥感图像和天空散射光遥感图像、分别进行地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像的地形变换。本发明的理论、方法较传统遥感数字处理先进、科学、合理，可广泛应用于遥感技术的定量研究和应用；其具体操作步骤如下(1)、卫星遥感数字图像与数字地形图配准；(2)、卫星遥感数字图像的大气修正；(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算；(4)、数字地形图上计算卫星遥感数字图像上各点地面(坡面)太阳直射光地形改正系数和地面(坡面)天空散射光地形改正系数(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比(6)、解析卫星遥感图像；生成太阳直射光遥感图像和天空散射光遥感图像；(7)、卫星遥感数字图像上山体与云层阴影的判断和卫星遥感数字图像上山体与云层阴影的消除；(8)、太阳直射光遥感图像和天空散射光遥感图像的水平地形变换；(9)、无阴影遥感正射数字图像的计算机生成；(10)、遥感图像的任意地形变换。
上述的骤(1)中的卫星遥感数字图像与数字地形图配准是GIS常规操作。
上述的步骤(2)中的卫星遥感数字图像的大气修正的步步骤如下(1)利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，按下式计算可以得到离散分布的像元大气程辐射遥感值。
DN(a)＝γ1*(DN1-DN2)/(γ1-γ2) (1)或者对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析也可以得到离散分布的像元大气程辐射遥感值；(2)对上述结果，再采用插值算法得到每个像元大气程辐射遥感值；对卫星遥感数字图像上各像元点遥感值DNij进行减去大气程辐射遥感值DNAij的运算，即DNij-DNAij，便完成了卫星遥感数字图像的大气修正。
其中，DNij代表卫星遥感数字图像上各像元点遥感值，DNAij代表大气程辐射遥感值，DNij-DNAij代表卫星遥感数字图像的大气修正值。
上述的步骤(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)-DNa)/(DN(x)-DN(y))(6)上述的步骤(6)中的计算卫星遥感数字图像上各点地面直射光地形影响系数按下式计算Fij＝1-tgαij·ctgθij·cosωij，ωij＝ALij-Aij(7)其中，αij代表像元地表坡角，ALij代表地表坡向，θij代表卫片注记中星下点太阳高度角，Aij代表太阳方位角，ωij代表太阳方位角与地面坡向之夹角。
上述的步骤(7)中的计算卫星遥感数字图像上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比按下式计算Gij=1-2nπΣk=1nβk,]]>式中n＝2π/Δt (8)
其中，βk为第K个方位上的最大平均遮蔽角即立体角。Δt为方位角步长，为360°/n。n为方位划分数，即n为按计算精度要求将360°划分的范围个数。
上述的步骤(8)中的解析卫星遥感图像；生成太阳直射光遥感图像和天空散射光遥感图像，方法如下太阳直射光遥感图像表达式DNSij=(DNij-DNAij)·Fij(Fij+Gij·Lij)---(9)]]>天空散射光遥感图像表达式DNDij=(DNij-DNAij)·Gij·LijFij+Gij·Lij---(10)]]>i＝1，2…m；j＝1，2…n其中，DNSij代表像元遥感数值中直射光分量，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Fij代表数字地形图上各点地面直射光地形影响系数，Gij代表数字地形图上各点地面天空散射光地形影响系数，DNDij代表像元遥感数值中散射光分量，LSij代表地面光谱反射亮度，LDij代表地面辐射照度。
上述的步骤(9)中的太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算DN′Sij=DNSij/Fij=(DNij-DNAij)/(Fij+Gij·Lij)---(11)]]>DN′Dij=DNDij/Gij=(DNij-DNAij)·Lij/(Fij+Gij·Lij)---(12)]]>i＝1，2…m；j＝1，2…n其中，DNij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，DNDij′代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
上述的步骤(10)中的卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数字运算
DN′ij=(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)---(13)]]>其中，DNij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNAij代表大气程辐射遥感值。
上述的步骤(11)中的卫星遥感数字图像无阴影正射影像的计算机生成，使用以下数字运算完成DNij′=DNSij′+DNDij′---(14)]]>DN′ij=(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)---(15)]]>其中，DNij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，DNDij′代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
上述的像元大气程辐射辐射DNAij，通过数字图像中相关信息中提取或地面观测数据计算得出；上述的像元地面太阳高度θij、方位角ALij由卫片注记给出的相关信息计算得出；上述的像元地面坡度αij、坡向Aij以及像元经、纬度由卫片对应DTM给出。上述的步骤(12)中的GIS支持下的卫星遥感数字图像的任意地形变换，使用以下数字运算完成XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DAij)×F`ij/(Fij+LijGij) (16)其中虚拟的坡面像元太阳直射光水平变换系数Fij′＝(1-tgαij′·ctgθij′·cosωij′)利用公式XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DAij)×F`ij/(Fij+LijGij)i＝1，2…m；j＝1，2…n (17)逐点计算便可完成太阳直射光遥感数字图像地形和辐射的虚拟与仿真。
其中，XDNSij代表像元遥感数值中直射光分量，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，Fij代表数字地形图上各点地面直射光地形影响系数，Fij′代表经大气修正后的数字地形图上各点地面直射光地形影响系数，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Gij代表数字地形图上各点地面天空散射光地形影响系数，αij′代表修正后的像元地表坡角，θij′代表修正后的卫片注记中星下点太阳高度角，ωij′代表修正后的太阳方位角与地面坡向之夹角。
本发明与现有技术相比，具有如下显而易见的突出实质性特点和显著优点由于卫星遥感成像是在自然光条件之下，自然光成份包括太阳直射光(方向性)分量和天空散射光(基本上是各向同性)，其相对强度以及在起伏地表上的再分配方式，对卫星遥感数据的影响有很大差异。将大气修正后卫星遥感图像解析，分离成太阳直射光遥感图像和天空散射光遥感图像，分别进行太阳直射光和天空散射光的变换，形成对应的太阳直射光和天空散射光的地形变换遥感数字图像，最后将二者合成自然光条件下的地形变换卫星遥感数字图像。其理论、方法较传统遥感数字图像处理(自然光为整体)先进、科学、合理。同时，所生成的太阳直射光和天空散射光遥感数字图像以及它们的地形变换数字图像都是自然界无法直接获取的，在理论上和实践上均有重要价值的遥感信息。
由于遥感数字图像的地形变换消除了地形起伏带来的地表光辐射差异以及云层与山体阴影的影响，在遥感信息中突出了地物波谱特征，使其更具可比性。从根本上改善了遥感数字图像的质量具有重要的科研和应用前景。
本项发明可广泛应用于遥感技术的变量研究与各种应用，如土地利用研究与调查、军事中的伪装与反伪装、阴影地物的识别、森林调查研究、卫片的模式识别以及遥感数字图像的虚拟与仿真等方面。
图1为本发明GIS支持下的卫星遥感数字图像的地形变换计算机生成流程图。
图2为福建龙岩地区1998.12.8经过大气修正后的TM遥感数字图像(1024×1024)；图3为与该卫片对应数字地形图(立体)；图4，5分别为卫星遥感数字图像上各点的太阳直射光和天空散射光地形改正系数及其可视化图像；图6，7分别为该卫片解析后生成的太阳直射光和天空散射光遥感数字图像；图8，9分别为图5，6的水平地形变换图像；
图10为图2的无阴影遥感正射数字图像；图11为图2任意地形变换后的遥感数字图像。
具体实施例方式
本发明的一个优选实施例结合附图详细说明如下本例为福建龙岩地区1998.12.8的卫片以GIS支持下的卫星遥感数字图像的地形变换。参见图1，本GIS支持下的卫星遥感数字图像的地形变换方法是首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像生成太阳直射光遥感图像和天空散射光遥感图像、分别进行地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像的地形变换。本发明的理论、方法较传统遥感数字处理先进、科学、合理，可广泛应用于遥感技术的定量研究和应用；其具体操作步骤如下(1)、卫星遥感数字图像与数字地形图配准；(2)、卫星遥感数字图像的大气修正；(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算；(4)、数字地形图上计算卫星遥感数字图像上各点地面(坡面)太阳直射光地形改正系数和地面(坡面)天空散射光地形改正系数(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比(6)、解析卫星遥感图像；生成太阳直射光遥感图像和天空散射光遥感图像；(7)、卫星遥感数字图像上山体与云层阴影的判断和卫星遥感数字图像上山体与云层阴影的消除；(8)、太阳直射光遥感图像和天空散射光遥感图像的水平地形变换；(9)、无阴影遥感正射数字图像的计算机生成；(10)、遥感图像的任意地形变换。
卫星遥感数字图像的地形变换在有多点实时辐射观测资料(散射/辐射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。上述方案中各步骤说明如下1、卫星遥感数字图像与数字地形图配准
卫星遥感数字图像与数字地形图(TEM)配准系GIS常规操作。
2、星遥感数字图像的大气修正获得像元大气程辐射遥感值的方法有两种利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，进以下计算DN(a)＝γ1*(DN1-DN2)/(γ1-γ2) (1)可以得到离散分布的像元大气程辐射遥感值。
对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析也可以得到离散分布的像元大气程辐射遥感值。
对上述结果，再采用插值算法可得到每个像元大气程辐射遥感值。对卫星遥感数字图像上各像元点遥感值进行减去大气程辐射遥感值的运算，即DNij-DNAij，便完成了卫星遥感数字图像的大气修正。
其中，DNij代表卫星遥感数字图像上各像元点遥感值，DNAij代表大气程辐射遥感值，DNij-DNAij代表卫星遥感数字图像的大气修正值。
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)-DNa)/(DN(x)-DN(y)) (6)6、计算数字地形图上各点地面直射光地形影响系数Fij＝1-tgαij·ctgθij·cosωij，ωij＝ALij-Aij(7)其中，αij代表像元地表坡角，ALij代表地表坡向，θij代表卫片注记中星下点太阳高度角，Aij代表太阳方位角，ωij代表太阳方位角与地面坡向之夹角。
7、计算数字地形图上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比Gij=1-2nπΣk=1nβk,]]>式中n＝2π/Δt (8)其中，βk为第K个方位上的最大平均遮蔽角即立体角。Δt为方位角步长，为360°/n。n为方位划分数，即n为按计算精度要求将360°划分的范围个数。
8、解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像，方法如下太阳直射光遥感图像表达式DNSij=(DNij-DNAij)·Fij(Fij+Gij·Lij)---(9)]]>天空散射光遥感图像表达式DNDij=(DNij-DNAij)·Gij·LijFij+Gij·Lij---(10)]]>i＝1，2…m；j＝1，2…n其中，DNSij代表像元遥感数值中直射光分量，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Fij代表数字地形图上各点地面直射光地形影响系数，Gij代表数字地形图上各点地面天空散射光地形影响系数，DNDij代表像元遥感数值中散射光分量，LSij代表地面光谱反射亮度，LDij代表地面辐射照度。
9、太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算DN′Sij=DNSij/Fij=(DNij-DNAij)/(Fij+Gij·Lij)---(11)]]>
DN′Dij=DNDij/Gij=(DNij-DNAij)·Lij/(Fij+Gij·Lij)---(12)]]>i＝1，2…m；j＝1，2…n其中，DDij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，DNDij′代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
10、GIS支持下的卫星遥感数字图像的地形变换，使用以下数字运算完成DNij′=DNSij′+DNDij′---(13)]]>DN′ij=(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)---(14)]]>其中，DNij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，DNDij′代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
11、卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数字运算DN′ij=(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)---(15)]]>其中，DNij代表像元遥感数值，DNij′代表经大气修正后的卫星像元遥感数据，DNAij代表大气程辐射遥感值。
上述的像元大气程辐射辐射DNAij，通过数字图像中相关信息中提取或地面观测数据计算得出；上述的像元地面太阳高度θij、方位角ALij由卫片注记给出的相关信息计算得出；上述的像元地面坡度αij、坡向Aij以及像元经、纬度由卫片对应DTM给出。
12、卫星遥感数字图像的任意地形变换，使用以下数字运算完成XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DAij)×F`ij/(Fij+LijGij) (16)其中虚拟的坡面像元太阳直射光水平变换系数Fij′＝(1-tgαij′·ctgθij′·cosωij′)利用公式XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DAij)×F`ij/(Fij+LijGij)i＝1，2…m；j＝1，2，…，n (17)逐点计算便可完成太阳直射光遥感数字图像地形和辐射的虚拟与仿真。
其中，XDNSij代表像元遥感数值中直射光分量，DNSij′代表经大气修正后的卫星像元遥感数据中直射光分量，Fij代表数字地形图上各点地面直射光地形影响系数，Fij′代表经大气修正后的数字地形图上各点地面直射光地形影响系数，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Gij代表数字地形图上各点地面天空散射光地形影响系数，αij′代表修正后的像元地表坡角，θij′代表修正后的卫片注记中星下点太阳高度角，ωij′代表修正后的太阳方位角与地面坡向之夹角。
卫星遥感数字图像的地形变换在有多点实时辐射观测资料(散射/辐射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。
由于卫星遥感成像是在自然光条件之下，自然光成份包括太阳直射光(方向性)分量和天空散射光(基本上是各向同性)，其相对强度以及在起伏地表上的再分配方式，对卫星遥感数据的影响有很大差异。将大气修正后卫星遥感图像解析，分离成太阳直射光遥感图像和天空散射光遥感图像，分别进行太阳直射光和天空散射光的变换，形成对应的太阳直射光和天空散射光的地形变换遥感数字图像，最后将二者合成自然光条件下的地形变换卫星遥感数字图像。其理论、方法较传统遥感数字图像处理(自然光为整体)先进、科学、合理。同时，所生成的太阳直射光和天空散射光遥感数字图像以及它们的地形变换数字图像都是自然界无法直接获取的，在理论上和实践上均有重要价值的遥感信息。
遥感数字图像的地形变换消除了地形起伏带来的地表光辐射差异以及云层与山体阴影的影响，在遥感信息中突出了地物波谱特征，使其更具可比性。从根本上改善了遥感数字图像的质量具有重要的科研和应用前景。
图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的合成)。
图11为图2任意地形变换后的遥感数字图像。
卫星遥感数字图像的地形变换在有多点同步实时辐射观测资料(散射/直射比)、卫星遥感数字图像高精度大气修正、典型地物反射率和相应高精度的数字地形图等条件下有最佳的地形变换效果。如果上述条件不能满足，仅仅利用一般精度的数字地形图进行卫星遥感数字图像地形变换也可以取得令人满意的较好效果。
权利要求
1.GIS支持下的卫星遥感数字图像的地形变换，其特征在于首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像生成太阳直射光遥感图像和天空散射光遥感图像、分别进行地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像的地形变换；其具体操作步骤如下(1)、卫星遥感数字图像与数字地形图配准；(2)、卫星遥感数字图像的大气修正；(3)、卫星遥感数字图像上各点太阳高度角、方位角的计算；(4)、数字地形图上计算卫星遥感数字图像上各点地面或坡面太阳直射光地形改正系数和地面或坡面天空散射光地形改正系数；(5)、观察测量或者计算卫星遥感数字图像上各点水平地面上的散射/辐射比；(6)、解析卫星遥感图像；生成太阳直射光遥感图像和天空散射光遥感图像；(7)、卫星遥感数字图像上山体与云层阴影的判断和卫星遥感数字图像上山体与云层阴影的消除；(8)、太阳直射光遥感图像和天空散射光遥感图像的水平地形变换；(9)、无阴影遥感正射数字图像的计算机生成；(10)、遥感图像的任意地形变换。
2.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(2)中的卫星遥感数字图像的大气修正的步骤如下(1)利用卫片相同地面上相邻两非同类像元DN1和DN2的地面反射率γ1和γ2，按下式计算可以得到离散分布的像元大气程辐射遥感值DN(a)＝γ1*(DN1-DN2)/(γ1-γ2)，或者对卫星遥感数字图像上山体和云层阴影中的水体像元遥感值进行分析得到离散分布的像元大气程辐射遥感值；(2)对上述结果，再采用插值算法得到每个像元大气程辐射遥感值；(3)对卫星遥感数字图像上各像元点遥感值DNij进行减去大气程辐射遥感值DNAij的运算，即DNij-DNAij，便完成了卫星遥感数字图像的大气修正，其中，DNij代表卫星遥感数字图像上各像元点遥感值，DNAij代表大气程辐射遥感值，DNij-DNAij代表卫星遥感数字图像的大气修正值。
3.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(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所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(4)中的卫星遥感数字图像上山体与云层阴影的判断方法如下利用与卫星遥感数字图像配准的数字地形图和卫片像元的太阳的位置参数，即高度角、方位角，进行山体与云层阴影的判断，其判断的准则是在太阳照射方向该像元点的最大地形、云层、高度角等于或者大于该像元点的太阳高度角即DH(i，j)≥θij，则该像元点为阴影；反之则不是阴影。
5.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(5)中的观测或计算卫星遥感数字图像上各点水平地面上的直、散射辐射比方法如下若卫片水平地面上有相邻2个同类像元X和Y，像元Y位于阴影中，DN(y)和DN(x)、DNa分别是其遥感值和大气程辐射遥感值，则该处水平地面上的散、直射辐射比L＝(DN(y)-DNa)/(DN(x)-DN(y))。
6.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(6)中的计算数字地形图上各点地面直射光地形影响系数按下式计算Fij＝1-tgαij·ctgθij·cosωij，ωij＝ALij-Aij，其中，αij代表像元地表坡角，ALij代表地表坡向，θij代表卫片注记中星下点太阳高度角，Aij代表太阳方位角，ωij代表太阳方位角与地面坡向之夹角。
7.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(7)中的计算数字地形图上各点地面天空散射光地形影响系数，即卫片像元地面天空散射光照射立体角与水平地面上的天空散射光照射立体角2π之比按下式计算Gij=1-2nπΣk=1nβk,]]>式中n＝2π/Δt，其中，βk为第K个方位上的最大平均遮蔽角即立体角，Δt为方位角步长，为360°/n，n为方位划分数，即n为按计算精度要求将360°划分的范围个数。
8.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(8)中的解析卫星遥感图像，生成太阳直射光遥感图像和天空散射光遥感图像，方法如下太阳直射光遥感图像表达式DNSij=(DNij-DNAij)·Fij(Fij+Gij·Lij),]]>天空散射光遥感图像表达式DNDij=(DNij-DNAij)·Gij·LijFij+Gij·Lij,]]>i＝1，2…m，j＝1，2…n，其中，DNSij代表像元遥感数值中直射光分量，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Fij代表数字地形图上各点地面直射光地形影响系数，Gij代表数字地形图上各点地面天空散射光地形影响系数，DNDij代表像元遥感数值中散射光分量，LSij代表地面光谱反射亮度，LDij代表地面辐射照度。
9.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(9)中的太阳直射光遥感数字图像和天空散射光遥感数字图像的地形变换，采用以下公式计算DN′Sij＝DNSij/Fij＝(DNij-DNAij)/(Fij+Gij·Lij)，DN′Dij＝DNDij/Gij＝(DNij-DNAij)·Lij/(Fij+Gij·Lij)，i＝1，2…m，j＝1，2…n，其中，DNij代表像元遥感数值，DN′ij代表经大气修正后的卫星像元遥感数据，DN′Sij代表经大气修正后的卫星像元遥感数据中直射光分量，DN′Dij代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
10.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(10)中的GIS支持下的卫星遥感数字图像的地形变换，使用以下数字运算完成DN′ij＝DN′Sij+DN′Dij，DN′ij＝(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)，其中，DNij代表像元遥感数值，DN′ij代表经大气修正后的卫星像元遥感数据，DN′Sij代表经大气修正后的卫星像元遥感数据中直射光分量，DN′Dij代表经大气修正后的卫星像元遥感数据中散射光分量，DNDij代表像元遥感数值中散射光分量，DNSij代表像元遥感数值中直射光分量，DNAij代表大气程辐射遥感值。
11.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(11)中的卫星遥感数字图像上山体与云层阴影的消除，采用卫星遥感数字图像上的阴影像元进行以下数字运算DN′ij＝(DNij-DNAij)·(1+Lij)/(Fij+Gij·Lij)，其中，DNij代表像元遥感数值，DN′ij代表经大气修正后的卫星像元遥感数据，DNAij代表大气程辐射遥感值。
12.根据权利要求2或3或6所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的像元大气程辐射辐射DNAij，通过数字图像中相关信息中提取或地面观测数据计算得出；所述的像元地面太阳高度θij、方位角ALij由卫片注记给出的相关信息计算得出；所述的像元地面坡度αij、坡向Aij以及像元经、纬度由卫片对应DTM给出。
13.根据权利要求1所述的GIS支持下的卫星遥感数字图像的地形变换，其特征在于所述的步骤(12)中的卫星遥感数字图像的任意地形变换，使用以下数字运算完成XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DNAij)×F`ij/(Fij+LijGij)，其中虚拟的坡面像元太阳直射光水平变换系数F′ij＝(1-tgα′ij·ctgθ′ij·cosω′ij)利用公式XDNsij＝DN`sij×F`ij＝DNsij×F`ij/Fij＝(DNij-DNAij)×F`ij/(Fij+LijGij)i＝1，2，…，m；j＝1，2，…，n逐点计算便可完成太阳直射光遥感数字图像地形和辐射的虚拟与仿真，其中，XDNSij代表像元遥感数值中直射光分量，DN′Sij代表经大气修正后的卫星像元遥感数据中直射光分量，Fij代表数字地形图上各点地面直射光地形影响系数，F′ij代表经大气修正后的数字地形图上各点地面直射光地形影响系数，DNij代表像元遥感数值，DNAij代表大气程辐射遥感值，Lij代表地面反射亮度，Gij代表数字地形图上各点地面天空散射光地形影响系数，α′ij代表修正后的像元地表坡角，θ′ij代表修正后的卫片注记中星下点太阳高度角，ω′ij代表修正后的太阳方位角与地面坡向之夹角。
全文摘要
本发明涉及GIS支持下的卫星遥感数字图像的地形变换。本方法是首先，建立卫星遥感信息与地形、太阳直射光和天空散射光之间定量关系的数理模型；继之解析卫星遥感数字图像生成太阳直射光遥感图像和天空散射光遥感图像、分别进行地形变换，消除起伏地形上的辐射差异以及这种差异对卫星遥感数字图像的影响；然后判断和消除遥感数字图像的山体与云层阴影；在此基础上最后完成卫星遥感数字图像的地形变换。本发明的理论、方法较传统遥感数字处理先进、科学、合理，可广泛应用于遥感技术的定量研究和应用。
文档编号G06T5/00GK101034472SQ20071003872
公开日2007年9月12日 申请日期2007年3月29日 优先权日2007年3月29日
发明者李先华, 师彪 申请人:上海大学