多光谱数据“特征空间”变换法信息提取
梅安新
摘 要:多光谱和高光谱遥感数据中,包含着丰富的表示地物特性的许多物理量。本文综述了利用多光谱和高光谱遥感数据,通过特征空间变换处理,构成如亮度、绿度、湿度、温度和水体的深度等不同物理量的方法。根据不同目标特性选择不同的空间信息提取方法,可获得比较好的效果。这个方法物理概念明确,有利于与专家知识的结合,有应用潜力,也可用于高光谱数据的信息提取。
关键词:多光谱,高光谱,特征空间变换,信息提取
一、引 言
遥感的信息获取向高空间分辨率和高光谱方向发展是一个重要的趋势,光谱信息的提取方法也随着光谱数量的增加而发生变化。通常的多光谱信息提取方法中,K-L、K-T变换得到较为普遍的应用,也被认为有较好的效果。
K-L变换(Karhunen-Loveve),亦称霍特林变换,是一种正交变换,在数学上称为主成分分析。图像处理中,K-L变换的目的是把多变量信息(如多光谱信息)压缩在一个新图像上,新图像是各个波段亮度值的加权和,按方差大小排序为第一、第二、第三……主分量。在大多数情况下,前三个主分量占97%以上,以后几个分量被视为噪声,可以丢弃不计。经变换压缩得到的新图像与老图像相比起到了滤波的作用,以R、G、B赋予各不相关的三个主分量,新图像的各分量之间相关性最小,起到了一定的分类作用。
K-L变换虽可压缩部分信息,但完全不考虑各波段本身的波谱效应,单纯用数学方法来提取多波段遥感信息,在许多情况下,与实际将有很大的出入。以TM来说,有7个通道,其覆盖范围从可见光、近红外到热红外,各波段的物理含义有很大的不同,同时,不同地区的图像内容也各不相同,其第一主分量信息的内涵也随图像不同而变化。如图像中水域面积远大于陆地时,从方差分布来看,K-L变换后的第一主分量必以水域信息为主。显然,在这种情况下,K-L变换中的第一主分量就不能称为“土壤线”(土壤轴)。
K-T变换(Kauth-Thomas)原先是对植物在一个生长期内反射光谱二维空间位置变化的描述,其光谱曲线构成图形如同一个“缨帽”,故称“缨帽”(Tasseled Cap)变换,其变换矩阵与K-L变换有些相似,也是一种空间坐标旋转的线性变换。原先K-T变换广泛用于MSS,1984年Crist和Cicone针对TM数据提出K-T变换时的B值变换矩阵,至今仍被广泛应用。与K-L变换不同的是K-T变换考虑到波段效应,但用一个不变的算法来应用与地球表面不同地区千变万化的地物信息,不可能都取得好效果。此外,K-T变换的各主分量的含义仍须根据图像中包含不同内容而转移。
无论是K-L变换还是K-T变换,原先只有数学上的意义,并无与地物波谱特征相应的、明确的物理意义,其物理含义需视图像本身的信息构成而定。人们在应用中将某些局部的结果视为普遍的规律而加以推广。例如,K-L、K-T变换中的第一主分量在二维坐标上的线均被称为“亮度线”(亮度轴)或“土壤线”(土壤轴);第二主分量均被称为“绿度线”或“植被线”;K-L变换的第三主分量未予命名,而K-T变换的第三主分量被称为“湿度线”;K-L变换的第四主分量被认为是“噪声”,K-T变换的第四主分量未明确命名(也有被称为“黄度轴”或“干扰度轴”等等)。
可以认为遥感图像的信息提取,必须依图像中包含的地物目标构成的不同而采取不同的算法。也就是说,图像信息的提取要取得较好效果,方法的选择必须根据不同图像上目标的信息特征、各种地物的相互关系、地区的特点和时相特征等。
二、原理和方法
研究表明,有些多波段遥感数据,只有3~4个波段,由于波段数太少、波段之间的相关性很高,无法构建多个物理含义明确的数据子集。当波段达到一定数量时,以TM为例,有7个波段的数据,至少可以构成亮度、绿度、湿度、深度和温度等5个特征空间的数据子集,各有其明确的物理含义。这一方法可称为特征空间变换法(Feature Space Transformation)。
在通常的信息提取中,最有用的是亮度、绿度和湿度三个特征值。各特征空间的数据集构成如下:
1.亮度BRI(Brightness)
在多光谱数据中,一种最简单的方法为各波段的求和。公式(1)~(4)为不同的求和方法:
式中:BRI1——亮度指数1;
ch1、ch2、chn——构成亮度值有关波段的像元亮度值,可以是TM1—3或其他能反映地表亮度的可见光波段。
其他形式还有:
式中:BRI2、BRI3、BRI4——亮度指数构成的形式2、3、4;
x1,x2,x3,…,xn为波段函数,代表物体的可见-近红总体辐射水平。
(1)、(2)式实际上是一种最小距离求和的方法,即求ch1-chn各波段对BRI的平均贡献。(3)式是在(1)、(2)式的基础上的一种简化计算方法。(4)式是一种加权平均求和的方法。考察λ~反差情况发现,各波段的加权和不仅保留了陆地信息,同时更能进一步抑制水体的灰度,增加陆地的反差,使BRI可更好地反映陆地裸露程度。另外,专业的经验亦可在权系数中得以体现,x1,x2,…,xn可由λ~反差表达,也可将水陆的灰度值归约到同等的灰度范围为标准来选择。K-T变换中的亮度指数,可视为(4)式的一个特例。
TM图像CCT上的亮度值与传感器接收到的辐射值之间有如下关系:
式中:V——CCT像元亮度值;
Dmax——CCT上最大量化级,对于TM为255;
Rmin值——对于Landsat5的TM来说,参见表1。
表1 Landsat5各波段传感器检测最大及最小辐射值(mW/cm2/Sr)
注:TM6的R与绝对温度的关系为:R=5.129 2×102 T-10-2 T+1.602 3,0=200K,255=340K.
2.绿度GRI(Greenness),其构成形式很多,有几十种,常用的有下列几种:
式中:GRI1——绿度指数形式1,这里又称比值植被指数(TRVI);
CHNIR——近红外像元亮度值;
CHR——红光波段像元亮度值。
显然,对于健康的植被来说,在红光波段由于叶绿素的吸收有一个低反射值,在近红外波段有一个高反射峰,两者比值愈大表明绿色植被的长势愈好。(5)式对于TM来说为TM4/TM3。经过这样的处理,图像上其他的地物的GRI值很小,而健康植被的GRI值很大,使得植被很容易在背景中区分出来。
这里,(6)、(7)式就是常用的“归一化差值植被指数”(NDVI)。与(5)式不同之处在于其分母的设置,使数据的变化范围受到约定,背景噪声的干扰较小,其植被长势的差别可更好地显示。
式中:GRI3.1——“农业植被指数”(AVI)或“环境植被指数”(EVI);
x——为待定系数,一般取2.0~2.4。
与此相似的还有GRI3.2:
(8)、(9)式统称为“差值植被指数”,差值愈大表示植被的长势愈好,也有利于不同植被的区分。
或者
上述两式中:X1,X2,…,Xn——为波段函数,可“+”可“-”;
CH1,CH2,…,CHn——为波段像元亮度值;
C——为常数,在线性相关坐标图上为直线的截距。
GRI4.1也被称为“综合绿度指数”(GR)。
对于Landsat TM:
对于MSS数据Kauth Thomas变换式的GVI值为:
类似形式被称为“垂直正交植被指数”(PVI)为:
显然,(12)、(13)、(14)式都是(10)式的特例。
其他特殊应用的关于绿度GRI值的算式,不完全统计有50~60种之多,例如,在国外林业中常用:
式中,SOIL5=-0.498+0.543MSS5+0.498MSS6;
SOIL6=2.734+0.498MSS5+0.498MSS6。
还有一些其他关于绿度值的计算方法,这里不再一一列举,估计今后还会有新的计算方法出来。但是,即使形式相同算法近似,采用不同的系数,结果也各异。在众多的绿度指标中,究竟何者为好,不能一概而论,必须根据植被的生长情况和时相、背景情况而定。
3.湿度WTI(Wetness)
由于水体对红外光的强烈吸收,在TM 7个波段中,TM5和TM7的图像水陆边界十分明显,土壤和植被水分的差别也不同程度地有所显示。因此,TM5和TM7图像上的像元亮度值可单独构成绿度指标,并可综合构成湿度数据子集。湿度值以WTI表示:
式中:WTI1,WTI2——湿度指数;
X——待定系数。
对于(23)式或当K-T变换时,Crist定义为:
Thrid=0.150 9TM1+0.197TM2+0.327 9TM3+0.340 6TM4-0.711 2TM5-0.457 2TM7
对于湿度指数,目前研究不及绿度指数研究得多,但从信息提取角度来看,当地面湿度大时,(18)式偏大,WTI5算法效果较好。在地面湿度不大时,WTI1、WTI2和WTI3信息提取的效果较好。
显然,以若干波段组成的特征空间数据子集,比单波段更能反映出物体的特征,在提取地面覆盖特征时具有明确的物理含义。具体处理过程中,通过波段选择和系数的调整,可取得满意的结果。
三、实例及效果分析
1.实例1
上海市嘉定城及周围地区,1987年5月18日Landsat TM1-5、TM7数据土地信息提取,先经3次褶积处理、Laplace变换,再进行特征空间变换(FST变换),按下列方案处理:BRI4(R),GRI2.2(G),WTI1(B)获一张似真彩色图像(图3),取得如下效果:
(1)BRI愈大颜色愈红,表示地面裸露程度高。如沪嘉高速公路、建筑工地、新建水泥顶房屋、未种植庄稼的土地和其他无植被覆盖的裸露土地等。
(2)GRI愈大,绿色愈深,表示地面植被良好。如林地、果园、生长健康的农作物、绿化等。
(3)WTI愈大者,蓝色愈深,表示水体或土壤中含水量较高的地物。如河流、坑塘、鱼塘等。水稻田和较湿的土地则显示不同深度的“蓝色→青色”,从而可以区别出水稻田和种植夏熟作物(麦子、油菜)的田地。
(4)菜地,呈“绿色→黄绿色→黄色”,显示出菜地不同生长状况和地面覆盖特征。长势好、覆盖密的呈绿色,主要分布在城市周围;一个田块内部分已经收割的菜地与裸土有相同的光谱特征,呈红色,而同一田块的另一部分仍有尚未收割长势良好的蔬菜则呈绿色,两者斑块不大呈镶嵌分布组成混合像元,按色光相加原理R+G=Y,而成黄色;由于混合像元内两种不同地物的面辐射贡献和面积贡献不同,使得菜地呈现出从“绿→黄绿→黄”等不同颜色。
(5)灰/红瓦屋顶的建筑物以及其他低反射率屋顶的建筑物,它们在各个波段的反射率都很低,图像上呈“灰→灰黑色”。
图1 K-T变换图像
图2 FST变换图像
图3 K-L变换图像
从这个例子可以看出无论是从信息量或者物理意义上,FST变换的效果显然好于KL变换(图1)和K-T变换(图2)。将其中一些主要地物置于亮度(BRI)、绿度(GRI)和湿度(WTI)的三维特征空间坐标上更能显示出差别,以利于进一步分类(图4)。
图4 不同物体在特征空间上的位置
2.实例2
上海市新建城区信息提取,图像数据同实例1。
我们分析了上海市新建城区和老城区图像表明,可鉴别的主要指标是屋顶材质和灰尘覆盖对反射率的影响,新建城区主要是水泥屋顶,且时间较新,灰尘覆盖少,老城区大都是瓦屋顶灰尘覆盖多反射率低,因而采用以下方案处理:BRI4(R),GRI3.1(G),GRI2.2(B)。由于GRI3.1与GRI2.2有微小差别,CHR和CHNIR在其中的作用不一样,强调了新建城区高亮度特征(图像上呈红色),也反映出老城区(绿→黑)和主要街道(蓝白)的差别。
结果表明,除新建城区的界线非常清楚地得到显示外,老城区内的若干较大的建筑工程,如上海铁路新客站、延安东路过江隧道入口、上海商城、静安希尔顿大楼……都清楚地反映出来。
3.实例3
长江口最大泥沙浑浊带信息提取,数据同前。
曾经有人用K-L变换做过这一地区的图像处理,对长江口的流态取得较好的效果,但未能获得泥沙浓度的有关信息。我们用特征空间变换法提取长江口泥沙浓度时,选取了以下三个特征空间集:
反射因子集BRI4——R;
透射因子集TF=XCHB——G;
吸收因子AF=。
对比“K-L变换”与“特征空间变换”信息提取结果看出,后者不仅清楚地显示出长江口的流态,而且清楚地突出了水体泥沙最大浑浊带的范围。实测表明,最大浑浊带表层含沙浓度一般在0.1~0.7kg/m2,洪季有时可达10kg/m2。北支最大浑浊带的上端在青龙港附近,下端在北支口门附近,整个北支堆积旺盛,日益变狭,黄瓜沙以下一系列水下沙坝露出水面或成沙洲链;南支江流作用强,最大浑浊带在口门以外,下端已出图幅;南槽最大浑浊带也在口门以外。
4.效果分析
以上方法表明,特征变换法信息提取的效果比K-L、K-T变换好,具体说来:
(1)物理概念明确。不像K-L变换完全按照“方差”大小排序,波段作用不明显。KT变换虽然保存波段信息特征,缺点是未考虑地面实况千差万别,采用单一算法。而特征空间变换方法提取信息,首先从物理概念出发,选择特征空间数据集。每一子空间数据集都有明确的物理含义。
(2)简化运算。特征空间数据子集,有的用简单的数学运算即可完成,比K-L变换的大矩阵乘积和方差计算快得多,且无需像K-T变换那样六七个波段都参加运算。从而减少了计算量,易于实现。
(3)有利于专家知识的加入。在波段选择和特征空间向量的算法选择上,均可由专家根据具体地物特征和知识来确定,把计算机计算和专家知识结合起来,使处理效果达到最佳。
(4)这种方法尤其对波段多、信息量大的数据非常有用。今后,成像光谱技术获得大量数据,通过特征空间数据子集选择,能减少大量运算,取得好效果。
参考文献
[1]John G.Lyon,Ding Yuan,Ross S.Lunetta,and Chris D.Elvidge.A Change Detection Experiment Using Vegetation Indices.Photogrammetric Engineering &Remote Sensing,1998.
[2]Carr,J.R..A Visual Basic Program for Principal Transformation of Digital Images.Computers &Geosciences,1998,24(3).
[3]Crist,E.P.and R.J.Kauth.The Tasseled Cap De-mystified.Photogrammetric Engineering &Remote Sensing,1998,52(1).
[4]Crist,E.P.and R.C.Cicone.Application of the Tasseled Cap Concept to Simulated Thematic Mapper Data.Photogrammetric Engineering &Remote Sensing,1998,50:343-352.
[5]Clausi,D.A.and M.E.Jernigan.A Fast Method to Determine Co-occurrence Texture Features.IEEE Transactions on Geoscience and Remote Sensing,1998,36(1):298~300.
[6]陈述彭,童庆禧,郭华东.遥感信息机理研究.北京:科学出版社,1998.
[7]Didan,K..MODIS Vegetation Index Production Algorithms,MODIS Vegetation Workshop,Missoula,Montana.Terrestrial Biophysics and Remote Sensing,2002.
[8](TBRS)MODIS Team,Tucson:University of Arizona,www.ntsg.umt.edu/MODIScon/index.html.
[9]戴昌达.遥感图像应用处理与分析.北京:清华大学出版社,2004.
[10]Holden,H.and E.LeDrew.Spectral Discrimination of Healthy and Non-healthy Coral BASED ON Cluster Analysis,Principal Component Analysis,and Derivative Spectroscopy,Remote Sensing of Environment,1998,65:217-224.
[11]John R.Jensen.Introductory Digital Image Processing,A Remote Sensing Perspective,Second Edition,Prentice Hall,Upper Saddle River,New Jersey,USA,1996.
[12]姜小光,唐伶俐,王长耀.高光谱数据的光谱信息特点及面向对象的特征参数选择——以北京顺义区为例.遥感技术与应用,2002,17(2):59-65.
[13]Kauth,R.J.and G.S.Thomas.The Tasseled Cap—A Graphic Description of the Spectral-Temporal Development of Agricultural Crops as Seen by Landsat,Proceedings,Symposium on Machine Processing of Remotely Sensed Data,West Lafayette,IN:LARSA,1976,41-51.
[14]梅安新,彭望琭,秦其明,刘慧平.遥感导论.北京:高等教育出版社,2002:122-127.
[15]田庆久,闵祥军.植被指数研究进展.地球科学进展,1998,13(4):327-333.
[16]Zhao,G.and A.L Maclean.A Comparison of Canonical Discriminant Analysis and Pricipal Component Analysis for Spectral Transformation.Photogrammetric Engineering &Remote Sensng,2000,66(7):841-847.
[17]赵英时.遥感应用分析原理与方法.北京:科学出版社,2003.
注:这篇文章以基于专家知识而不是简单地基于图像光谱角度,提出作者自己的遥感信息提取方法。随着高光谱遥感的发展,其信息提取方法有巨大的发展空间,更需要专家知识的结合。
免责声明:以上内容源自网络,版权归原作者所有,如有侵犯您的原创版权请告知,我们将尽快删除相关内容。