CN104809704B - 基于地形改正的卫星影像虚假地形感知现象改正方法 - Google Patents

Abstract

本发明公开了一种基于地形改正的卫星影像虚假地形感知现象改正方法，从改变光照方向的角度出发，引入地形改正中的C校正模型，通过建立原始影像与虚假地形感知现象改正后影像经过地形改正后的辐射关系，获得虚假地形感知现象改正前后影像辐射值的对应关系，进而推导出改正后影像灰度值。本发明方法能够在改正虚假地形感知现象的同时，较好地保持原始影像的空间细节及光谱特性，可显著提高影像解译准确度。

Description

基于地形改正的卫星影像虚假地形感知现象改正方法

技术领域

本发明属于测绘产品之一的正射影像产品自动化生产技术领域，尤其涉及一种基于地形改正的卫星影像虚假地形感知现象改正方法。

背景技术

近年来，随着数字地球的发展，众多国际国内公司、企业及政府部门相继推出了各具特色的地理信息服务平台。此类平台中使用的底图数据大部分来源于可见光遥感影像，如LANDSAT系列、SPOT系列、IKONOS、Quick Bird等。可见光遥感卫星一般使用太阳同步回归轨道，成像时间一般为过境地方时的上午，此时太阳位于东南方向，高大山脊等地物产生的阴影会出现在它的西北方向^[1-4]。由于地图通常采用上北、下南、左西、右东的方式表达地理方位，在可见光遥感影像地图上，阴影就出现在了山脊的左上方。人类日常生活中，太阳光线永远位于物体的上方，并在物体的下方产生阴影，导致人类在大脑中形成阴影永远位于物体下方的定势思维^[5,6]。因此在对遥感影像地图进行判读时，人们认为阴影来源于地图上阴影上方的物体，并由此判断阴影上方区域为高程突出区域，实际情况却是高程突出区域位于阴影的下方，导致人们把山谷误认为山脊，反之亦然，这就是可见光遥感影像的虚假地形感知现象(False Topographic Perception Phenomenon,FTPP)。虚假地形感知现象会造成人对影像中真实地形的错误理解，若不进行恰当处理，会导致公众对影像产品的真实性和可靠性产生质疑，甚至导致在重大决策中出现错误。因此，对虚假地形感知现象进行改正是非常有必要的。

形成虚假地形感知现象的原因大致可归为观察者视角、地表特征、光照方向3类。据此，可将现有改正方法分为以下几类：1)基于观察者视角的改正方法。以Saraf^[7]将影像旋转180°的方法为代表，通过变换观察者视角来改正虚假地形感知现象。2)基于地表特征的改正方法。通过将原始地形反转或添加正向地形信息来获得正确的立体视觉效果，如灰度反转法^[7]、亮度分量反转法^[1]、阴影模型替代法^[2]、小波变换融合法^[8]、阴影模型叠加法^[9]等。3)基于光照方向的改正方法。章皖秋等^[10]引入地形改正中的cosine模型，将光照由东南方改为竖直方向，再由竖直方向改为西北方向，达到改正反立体的目的。4)其它改正方法。以人工处理法为代表，主要依靠图像处理软件对影像进行处理，如Patterson^[11]使用的Photoshop中仿制图章工具等。

文中涉及的参考文献如下：

[1]Saraf,A.K.,Sinha,S.T.,Ghosh,P.,&Choudhury,S.(2007).A new techniqueto remove false topographic perception phenomenon and its impacts in imageinterpretation.International Journal of Remote Sensing,28(5),811-821.

[2]Saraf,A.K.,Ghosh,P.,Sarma,B.,&Choudhury,S.(2005).Cover:Developmentof a new image correction technique to remove false topographic perceptionphenomena.International journal of remote sensing,26(8),1523-1529.

[3]Patterson,T.(2002).Getting real:Reflecting on the new look ofNational Park Service maps.Cartographic perspectives,(43),43-56.

[4]Rudnicki,W.(2000).The new approach to the relief shading appliedin the satellite image maps.Proceedings of High Mountain Cartography,18,105-106.

[5]Ramachandran,V.S.(1988).Perception of shape from shading.Nature.

[6]Murray,J.(1994).Some perspectives on visual depth perception.ACMSIGGRAPH Computer Graphics,28(2),155-157.

[7]Saraf,A.K.,Das,J.D.,Agarwal,B.,&Sundaram,R.M.(1996).Falsetopography perception phenomena and its correction.International Journal ofRemote Sensing,17(18),3725-3733.

[8]Wu,B.,Li,H.,Gao,Y.(2013).Investigation and remediation of falsetopographic perception phenomena observed on Chang’E-1lunar imagery.Planetaryand Space Science,75,158-166.

[9]Gil M L,Arza M,Ortiz J,et al.DEM shading method for the correctionof pseudoscopic effect on multi-platform satellite imagery[J].GIScience&Remote Sensing,2014,51(6):630-643.

[10]章皖秋,袁华,徐天蜀.(2010).遥感影像正立体化研究.北京林业大学学报,(3),41-46.

[11]Patterson,T.,Kelso,N.V.(2004).Hal Shelton revisited:Designing andproducing natural-color maps with satellite land cover data.CartographicPerspectives,(47),28-55.

发明内容

针对现有的FTPP改正技术中存在的影像信息损失问题，本发明提供了一种基于地形改正的卫星影像虚假地形感知现象改正方法，该方法在获得正确的立体视觉效果的同时，还能最大程度保持原始影像的空间信息与光谱特性。

为解决上述技术问题，本发明采用如下的技术方案：

一、基于地形改正的卫星影像虚假地形感知现象改正方法，包括步骤：

步骤1，基于原始影像对应的太阳方位角及与数字高程模型，获得原始影像中各像素在东南光照方向与西北光照方向的太阳入射角；

步骤2，基于地形改正中的C模型，根据原始影像各像素的起伏地面辐射值与东南光照方向的太阳入射角，获得原始影像的C校正参数；

步骤3，基于地形改正原理获得虚假地形感知现象改正后影像(以下简称“改正后影像”)的C校正参数，本步骤进一步包括：

3.1基于原始影像和改正后影像的地形改正结果相同的原则，依据C模型获得原始影像和改正后影像的起伏地面辐射值关系；

3.2基于原始影像和改正后影像的起伏地面辐射值的关系，以原始影像和改正后影像所有像素的起伏地面辐射值均值相等为约束，利用原始影像像素在东南光照方向与西北光照方向的太阳入射角、以及原始影像C校正参数，获得改正后影像的C校正参数；

步骤4，基于原始影像和改正后影像的地形改正结果相同的原则，利用改正后影像的C校正参数获得改正后影像像素的灰度值。

步骤2进一步包括子步骤：

2.1依据C模型，根据原始影像各像素的起伏地面辐射值L_T与东南方向光照的太阳入射角i进行直线拟合，获得线性回归方程L_T＝mcos(i)+b的截距b及斜率m；

2.2根据获得原始影像的C校正参数c。

步骤3中所述的改正后影像的C校正参数c₂为：

其中，n为影像像素总数，k表示像素编号，k＝0,1,...,n-1；和分别为原始影像像素k在东南方向光照和西北方向光照的太阳入射角；为原始影像像素k的起伏地面辐射值；sz为太阳天顶角；c₁为原始影像的C校正参数。

步骤4进一步包括：

4.1基于原始影像和改正后影像的地形改正结果相同的原则获得改正后影像像素k的起伏地面辐射值其中，为原始影像像素k的起伏地面辐射值，sz为太阳天顶角，c₁、c₂分别为原始影像和改正后影像的C校正参数，和分别为原始影像像素k在东南方向光照和西北方向光照的太阳入射角；

4.2根据改正后影像像素起伏地面辐射值获得像素灰度值。

二、基于地形改正的卫星影像虚假地形感知现象改正系统，包括：

第一模块，用来基于原始影像对应的太阳方位角及与数字高程模型，获得原始影像中各像素在东南光照方向与西北光照方向的太阳入射角；

第二模块，用来基于地形改正中的C模型，根据原始影像各像素的起伏地面辐射值与东南光照方向的太阳入射角，获得原始影像的C校正参数；

第三模块，用来基于地形改正原理获得改正后影像的C校正参数，本模块进一步包括：

第一子模块，用来基于原始影像和改正后影像的地形改正结果相同的原则，依据C模型获得原始影像和改正后影像的起伏地面辐射值关系；

第二子模块，用来基于原始影像和改正后影像的起伏地面辐射值的关系，以原始影像和改正后影像所有像素的起伏地面辐射值均值相等为约束，利用原始影像像素在东南光照方向与西北光照方向的太阳入射角、以及原始影像C校正参数，获得改正后影像的C校正参数；

第四模块，用来基于原始影像和改正后影像的地形改正结果相同的原则，利用改正后影像的C校正参数获得改正后影像像素的灰度值。

光学影像中，由于地形起伏的影响，地表向阳面受到较多光照，背阳面受到光照不足，从而导致影像分类过程中出现错误，即向阳面和背阳面的同一种地物被分为两类。为改正这种现象，人们提出了地形改正的概念。地形改正，也叫做地形归一化或地形光照改正，是指通过建立一系列模型来改正由地形起伏引起的地表光照不均的现象。地形改正通过将光照方向改为垂直入射来消除影像上的地形效应，使原本不平整的地表看起来像一个平面。

本发明基于虚假地形感知现象成因，以光照方向作为出发点，提出了一种基于地形改正的卫星影像虚假地形感知现象自动改正方法，该方法以地形改正为中介，由于对同一地表在光照强度和太阳高度不变的情况下进行地形改正的结果相同，得以建立原始影像与虚假地形感知现象改正后影像在地形改正后的辐射关系，进而推导出虚假地形感知现象改正前后影像辐射值的对应关系。

与现有技术相比，本发明具有如下特点：

本发明从改变光照角度出发，引入地形改正模型，将光照方向由东南向变换为西北向；并建立原始影像与虚假地形感知现象改正后影像的定量关系，进而获得虚假地形感知现象改正后影像的灰度值。本发明采用C模型作为地形改正模型，可有效避免由亮度较低处的地形过度改正而引起的虚假地形现象改正错误。

采用本发明方法改正后的影像，能获得正确的立体视觉效果；另外，本发明能在改正虚假地形感知现象的同时，较好地保持原始影像的空间细节及光谱特性，从而显著提高影像解译准确度。

附图说明

图1为本发明方法的具体流程图；

图2为太阳高度角与太阳方位角示意图；

图3为本发明原理示意图，图(a)～(c)所示地表相同但光照方向不同，其中，图(a)为东南光照方向，图(b)为西北光照方向，图(c)为竖直光照方向；

图4为地形效果对比图，图(a)为原始多光谱影像，图(b)为相应DEM在西北向光照下的阴影模型；

图5为不同方法的改正结果，其中图(a)为阴影模型替代法改正结果，图(b)为亮度分量反转法改正结果，图(c)为阴影模型叠加法改正结果，图(d)为本发明方法改正结果；

图6为不同改正方法结果的细节对比图。

具体实施方式

本发明从改变光照方向角度出发，以地形改正为中介推导出原始影像以及FTPP改正后影像在地形改正后的关系，进而获得原始影像与FTPP改正后影像辐射值的对应关系，获得FTPP改正后影像灰度值，从而实现虚假地形感知现象的自动改正。

为了更好理解本发明的技术方案和技术效果，下面将结合附图对本发明做进一步说明。见图1，本发明方法具体步骤如下：

步骤1，利用原始影像对应的太阳方位角及与原始影像配准的数字高程模型(Digital Elevation Model,DEM)，获取原始影像中各像素在东南方向光照与西北方向光照的太阳入射角。

原始影像成像时的太阳方位角可由影像元数据获得。若影像元数据中不包含太阳方位角，也可由影像成像的时间及经纬度，通过天文学公式计算得到，见公式(1)～(3)。如图2所示，图中h_S为太阳高度角，A_S为太阳方位角，太阳高度角和太阳方位角的计算公式如下：

sin h_S＝sinφsinδ+cosφcosδcosΩ (1)

式(1)～(2)中，φ为当地纬度；δ为成像当日太阳赤纬角；Ω为太阳时角。以上所有数据均可通过成像时刻及星下点经纬度数据换算得到。

太阳入射角i的计算方法如下：

式(3)中，S为坡度；sz为太阳天顶角，与太阳高度角h_S互为余角；为坡向角。

原始影像对应的DEM的格网间距需要与原始影像分辨率一致，通常通过基础测绘获得。将经过外业控制点精确定向后的卫星影像利用DEM进行正射纠正处理，得到正射影像以及与之配套的DEM。

步骤2，基于地形改正中的C模型，将原始影像各像素的辐射值与东南方向光照的太阳入射角进行最小二乘拟合，求解原始影像的C校正参数。

地形改正中，本发明采用一种经典且普适的非朗伯体反射模型—C模型进行试验。许多实验证明，C模型能够在可见光范围内比较完整地保持原始影像的光谱特性，并且可以提高影像分类精度。C模型由Teillet提出，其基本原理如下：

式(4)中，L_H表示水平面辐射值；L_T表示起伏地面辐射值；c即C校正参数；sz为太阳天顶角；i为太阳入射角。

c的计算方法如下：

b和m分别是线性回归方程(6)的截距和斜率：

L_T＝m cos(i)+b (6)

本步骤具体包括：

2.1依据C模型，将原始影像各像素的起伏地面辐射值与东南方向光照下太阳入射角i的余弦进行直线拟合，求得线性回归方程(6)的截距b及斜率m，起伏地面辐射值可由像素灰度值计算获得。具体实施时，拟合方法可采用最小二乘拟合法、或整体最小二乘法中的SVD法、迭代法等。

2.2根据公式(5)，利用线性回归方程(6)的截距b及斜率m，获得原始影像的C校正参数。

步骤3，利用步骤1获得的东南方向光照与西北方向光照的太阳入射角和步骤2获得的C校正参数，获得虚假地形感知现象改正后影像的C校正参数。

形成FTPP的因素主要有光照、地形、观察者视角三类。在不改变观察者视角的前提下，可以从光照方向出发探求虚假地形感知现象的改正方法。通过将光照方向由东南向转变为西北向，使原始影像的光照-阴影关系发生相应变化，从而达到改正虚假地形感知现象的目的。由于直接转变光照方向比较困难，此处引入地形改正方法，以其作为中介来推导FTPP改正前后影像辐射值的关系。如图3所示，对于完全相同的地表，图3(a)是对东南方向光照，即原始影像光照方向的模拟；图3(b)是对西北方向光照，即FTPP改正后影像光照的模拟，其中，东南方向光照与西北方向光照对应的太阳高度角相同，西北方向光照的太阳方位角由东南方向光照的太阳方位角加180°获得；图3(c)是对竖直方向光照，即地形改正后影像光照的模拟。

利用C模型，分别对东南光照方向及西北方向光照下的影像进行地形改正，可以得到：

式(7)～(8)中，L_H1和L_H2分别为东南方向光照和西北方向光照下水平面辐射值；L_T1和L_T2分别为东南方向光照和西北方向光照下起伏地面辐射值；i₁和i₂分别为东南方向光照和西北方向光照下太阳入射角；c₁和c₂分别为东南方向光照和西北方向光照下的C校正参数，c₁即步骤2获得的C校正参数；sz为太阳天顶角。

由于地形改正的作用是将倾斜光照改为竖直光照，因此不难推断，对于同样的地表，在光照强度和太阳高度不变的情况下，将东南方向光照与西北方向光照下的影像均改为竖直方向光照后的结果相同，即：

L_H1＝L_H2 (9)

由式(9)即可获得FTPP改正前后影像起伏地面辐射值的关系：

由式(10)可知，对于一幅像素总数为n的影像，可以列出以下n个方程：

其中，与分别代表第0、1、...n-1个像素在东南方向光照和西北方向光照下起伏地面辐射值；与分别代表第0、1、...n-1个像素在东南方向光照和西北方向光照下的太阳入射角。

由式(11)可知，方程组中含有n个方程、(n+1)个未知数方程组的解不唯一。因此，需要添加一定约束条件将方程组的解唯一化。出于保持FTPP改正前后影像整体色调一致的原则，本发明中将约束条件设定为FTPP改正前后影像灰度均值相等。由于FTPP改正前后影像像素总数相等，因此约束条件等同于：

将方程组中单独移到等号一侧，再将所有方程相加，整理可得：

由式(12)～(13)可得：

对式(14)进行整理，即可求得改正后影像的C校正参数c₂的值：

其中n为影像像素总数，k表示像素编号，k＝0,1,...,n-1；sz为太阳天顶角，与太阳高度角互为余角；为原始影像像素k在东南向光照的太阳入射角；为原始影像像素k在西北向光照的太阳入射角；c₁为原始影像的C校正参数；为原始影像像素k的起伏地面辐射值。

将步骤1和步骤2获得的已知量代入上式，即可求得改正后影像的C校正参数c₂。

步骤4，利用步骤1～3获得参数，计算改正后影像的灰度值，即完成影像的虚假地形感知现象改正。

由改正后影像的C校正参数c₂，依据C模型，即可获得FTPP改正后影像各像素起伏地面辐射值：

将步骤1～3获得的数据，即c₁,c₂,sz,和代入式(16)，即完成FTPP改正。

根据像素起伏地面辐射值即可获得像素灰度值。

图4为地形效果对比图，图(a)为原始多光谱影像，图(b)为相应DEM在西北向光照下的阴影模型。从图中可以明显看出，多光谱影像中存在地形反转现象，原始影像中的山脊在阴影模型中变成了山谷，反之亦然。

图5为不同方法的改正结果，其中图(a)为阴影模型替代法(简称为“方法1”)改正结果，图(b)为亮度分量反转法(简称为“方法2”)改正结果，图(c)为阴影模型叠加法(简称为“方法3”)改正结果，图(d)为本发明方法(简称为“方法4”)改正结果。

从图5中不难看出，几种方法对于FTPP均有不同程度的改正。为了便于观察及说明，选取3处细节进行对比，见图6，01～03分别表示原始影像中细节1、细节2、细节3，11～13分别表示方法1改正后的细节1、细节2、细节3，21～23分别表示方法2改正后的细节1、细节2、细节3，31～33分别表示方法3改正后的细节1、细节2、细节3，41～43分别表示方法4改正后的细节1、细节2、细节3。可以看出，仅就反立体的改正效果而言，几种方法中，本发明方法、晕渲图叠加法及阴影模型替代法的效果较好，从改正后影像中能够明显感受出与原始影像相反的立体，且改正后影像立体感较强。其次是亮度分量反转法，在改正后影像中能够感受到原始光照-阴影位置关系的变化，进而感受到立体的反转，但改正后影像在立体感上略逊于其他方法。

就细节信息的保持而言，本发明方法与亮度分量反转法效果最优，能够在改正反立体的同时最大程度地保留原始影像中的细节信息。例如，图6中细节1中线状的道路、位于细节2上部的田地，这些原始影像中的细节信息在改正后影像中仍然能够完好地体现出来。其次为晕渲图叠加法与阴影模型替代法，原始影像中的部分细节在改正后变得模糊不清，有些细节甚至难以辨认，如细节1中的道路等。

在改正反立体的同时保持原始影像的光谱特性并不容易。由图5可知，利用不同方法改正后的影像与原始影像相比均存在一定程度的色差。从整体上来看，本发明方法能够使改正前后影像的颜色改变最少。与原始影像相比，改正后影像在色调上虽然有所差异，但仍然能保持整体上的一致性。其次为晕渲图叠加法，能够保持原始影像的部分光谱信息，但改正后影像的对比度和饱和度有较大程度的降低。阴影模型替代法与亮度分量反转法在光谱特性的保持上不是很理想，两种方法在整体上出现了较明显的偏色情况，原始影像的光谱信息发生了严重扭曲。

综上，本发明方法能够在有效改正反立体的同时，最大程度地保持原始影像的空间细节及光谱特性。

Claims (5)

1.基于地形改正的卫星影像虚假地形感知现象改正方法，其特征在于，包括步骤：

步骤1，基于原始影像对应的太阳方位角及与数字高程模型，获得原始影像中各像素在东南光照方向与西北光照方向的太阳入射角；

步骤2，基于地形改正中的C模型，根据原始影像各像素的起伏地面辐射值与东南光照方向的太阳入射角，获得原始影像的C校正参数；

步骤3，基于地形改正原理获得改正后影像的C校正参数，本步骤进一步包括：

3.1基于原始影像和改正后影像的东南方向光照及西北方向光照下地形改正结果相同的原则，依据C模型获得原始影像和改正后影像的起伏地面辐射值关系；

3.2基于原始影像和改正后影像的起伏地面辐射值的关系，以原始影像和改正后影像所有像素的起伏地面辐射值均值相等为约束，利用原始影像像素在东南光照方向与西北光照方向的太阳入射角、以及原始影像C校正参数，获得改正后影像的C校正参数；

步骤4，基于原始影像和改正后影像的地形改正结果相同的原则，利用改正后影像的C校正参数获得改正后影像像素的灰度值。

2.如权利要求1所述的基于地形改正的卫星影像虚假地形感知现象改正方法，其特征在于：

步骤2进一步包括子步骤：

2.1依据C模型，根据原始影像各像素的起伏地面辐射值L_T与东南方向光照的太阳入射角i进行直线拟合，获得线性回归方程L_T＝m cos(i)+b的截距b及斜率m；

2.2根据获得原始影像的C校正参数c。

3.如权利要求1所述的基于地形改正的卫星影像虚假地形感知现象改正方法，其特征在于：

步骤3中所述的改正后影像的C校正参数c₂为：

$c_2=\frac{{\Σ}_{k=0}^{n-1}cos\left(i_{2_k}\right)L_{T{1}_k}\frac{cos\left(sz\right)+c_1}{cos\left(i_{1_k}\right)+c_1}-cos\left(sz\right){\Σ}_{k=0}^{n-1}{L}_{T{1}_k}}{{\Σ}_{k=0}^{n-1}{L}_{T{1}_k}-{\Σ}_{k=0}^{n-1}{L}_{T{1}_k}\frac{cos\left(sz\right)+c_1}{cos\left(i_{1_k}\right)+c_1}}$

其中，n为影像像素总数，k表示像素编号，k＝0,1,...,n-1；和分别为原始影像像素k在东南方向光照和西北方向光照的太阳入射角；为原始影像像素k的起伏地面辐射值；sz为太阳天顶角；c₁为原始影像的C校正参数。

4.如权利要求1所述的基于地形改正的卫星影像虚假地形感知现象改正方法，其特征在于：

步骤4进一步包括：

4.1基于原始影像和改正后影像的地形改正结果相同的原则获得改正后影像像素k的起伏地面辐射值其中，为原始影像像素k的起伏地面辐射值，sz为太阳天顶角，c₁、c₂分别为原始影像和改正后影像的C校正参数，和分别为原始影像像素k在东南方向光照和西北方向光照的太阳入射角；

4.2根据改正后影像像素起伏地面辐射值获得像素灰度值。

5.基于地形改正的卫星影像虚假地形感知现象改正系统，其特征在于，包括：

第一模块，用来基于原始影像对应的太阳方位角及与数字高程模型，获得原始影像中各像素在东南光照方向与西北光照方向的太阳入射角；

第二模块，用来基于地形改正中的C模型，根据原始影像各像素的起伏地面辐射值与东南光照方向的太阳入射角，获得原始影像的C校正参数；

第三模块，用来基于地形改正原理获得改正后影像的C校正参数，本模块进一步包括：

第一子模块，用来基于原始影像和改正后影像的东南方向光照及西北方向光照下地形改正结果相同的原则，依据C模型获得原始影像和改正后影像的起伏地面辐射值关系；

第二子模块，用来基于原始影像和改正后影像的起伏地面辐射值的关系，以原始影像和改正后影像所有像素的起伏地面辐射值均值相等为约束，利用原始影像像素在东南光照方向与西北光照方向的太阳入射角、以及原始影像C校正参数，获得改正后影像的C校正参数；

第四模块，用来基于原始影像和改正后影像的地形改正结果相同的原则，利用改正后影像的C校正参数获得改正后影像像素的灰度值。