CN108427964B - 一种遥感图像与地球化学的融合方法及系统 - Google Patents
一种遥感图像与地球化学的融合方法及系统 Download PDFInfo
- Publication number
- CN108427964B CN108427964B CN201810177878.0A CN201810177878A CN108427964B CN 108427964 B CN108427964 B CN 108427964B CN 201810177878 A CN201810177878 A CN 201810177878A CN 108427964 B CN108427964 B CN 108427964B
- Authority
- CN
- China
- Prior art keywords
- image data
- remote sensing
- geochemical
- image
- covariance matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
- G06F18/251—Fusion techniques of input or preprocessed data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Astronomy & Astrophysics (AREA)
- Remote Sensing (AREA)
- Multimedia (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种遥感图像与地球化学的融合方法及系统。方法包括:首先,对所述遥感图像进行最小噪声变换,获得多维图像数据;对所述地球化学图像进行点面转化,获得地球化学栅格图像数据,即使地球化学图像变成栅格数据;然后将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;并利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据,最后对替换后的多维图像数据进行最小噪声反变换,得到融合后的遥感化学图像,实现了地球化学和遥感数据的融合。
Description
技术领域
本发明涉及图像处理领域,特别涉及一种遥感图像与地球化学的融合方法及系统。
背景技术
遥感图像是栅格数据,而地球化学是点阵数据,通常遥感图像和地球化学的融合不能够很好的实现,仅仅是形式的融合,比如常规的主分类融合方法、彩色变换融合方法,融合的图像可识别程度很低。
遥感地球化学是利用遥感的手段反演元素丰度特征,因矿物的吸收特征主要是分子键振动形成的吸收特征或者是电子跃迁形成的吸收特征。主要的吸收特征(主要指吸收谷的波长)与分子振动和电子跃迁有关,这与分子晶体结构有关,如羟基、碳酸盐基团的吸收特征,和铁离子、锰离子的吸收特征,而与元素丰度关系不大。但是元素含量对分子吸收特征是有影响的,如高铝白云母和低铝白云母,其吸收特征有一定的差异性,铝含量的增加吸收谷向左偏移 (Swayze,1997)。在月球上的铁、钛元素含量与月壤光谱特征有关系,Lucey 给出了经验公式(Lucey et al.,1998)。这就说明元素含量对光谱有一定的影响。
由于可见光-短波红外的光谱都是矿物(或岩石)的分子振动和电子跃迁产生的。阳离子的大小和含量对矿物波谱特征有影响,随着矿物(或岩石)成分从酸性向基性过渡,它们的主要的吸收谷向在短波红外区域向波长较长的区域偏移。例如钠-白云母、钾-白云母,镁硅白云母、镁铁-多硅白云母,吸收峰依次在2115nm;2205nm;2220nm和2230nm。
参见附图1:为了分析元素对光谱特征的影响,从USGS最新谱库中获得白云母的光谱特征和元素含量相对特征,采用白云母AL含量的不同来进行分析,单一的蚀变矿物元素含量不同对光谱特征影响的区别如表1所示。参见附图2中可以看出,在热红外区,AL含量虽对波谱吸收特征有影响,但许多因素使得吸收特征变得复杂,编号为G108和G111白云母吸收特征就不明显。因此主要是利用短波红外区域进行区分。
表1
白云母编号 | AL2O3元素含量 | 吸收特征 |
G120 | 35.69 | 1.412191,2.200405 |
G119 | 33.34 | 1.411807,2.201339 |
G118 | 31.31 | 1.413346,2.207902 |
G117 | 35.54 | 1.4112305,2.1999385 |
G116 | 31.06 | 1.414503,2.213559 |
G114 | 33.09 | 1.412961,2.205085 |
G113 | 35.64 | 1.411807,2.197607 |
G111 | 28.03 | 1.416435,2.2197195 |
G108 | 35.0 | 1.410462,2.198539 |
从表1和图3中都可以看出白云母Al的含量与白云母光谱吸收谷位置特征是有一定的关系。为了分析白云母Al的含量与白云母光谱吸收谷位置特征, 首先分别对1.4μm左右和2.2μm左右这两处的位置特征做和Al含量进行做一个二维数据表,然后利用最小二乘法进行拟合,可以得出铝含量与位置特征的关系公式是典型的线性关系,铝含量与2.2μm左右吸收谷位置相关性为 0.9274,铝含量与1.4μm左右吸收谷位置相关性为0.8443,从相关性可以看出,铝含量与2.2μm左右吸收谷位置更加相关,那么铝含量与光谱吸收谷位置关系函数如下:
为了更好的表示元素含量与光谱吸收特征的关系,从USGS最新谱库上获得橄榄石的光谱特征与元素含量如表2所示,分析Fe的含量与橄榄石光谱吸收谷位置的关系,主要是2.2μm左右这个位置的特征,做一个二维数据表,然后利用最小二乘法进行拟合,得出橄榄石铁含量与吸收谷位置特征的关系公司,是典型的线性关系,铁含量与2.2μm左右吸收谷位置相关性为0.975,相关性非常好(图4)。铁含量与光谱吸收位置关系函数如下:
ConFeO=0.0003pos2.2+2.1919
式中,ConFeO是FeO的含量,pos2.2是橄榄石在2.2μm左右处的吸收位置。
表2
橄榄石编号 | FeO的含量 | 吸收特征 |
KI3005 | 62.82 | 2.214504 |
KI3054 | 30.59 | 2.202275 |
KI3188 | 41.43 | 2.207432 |
KI3189 | 34.63 | 2.204147 |
KI3291 | 53.65 | 2.209784 |
KI3377 | 59.75 | 2.212614 |
KI4143 | 47.65 | 2.207902 |
对于岩石来说,由于岩石是矿物的组合,是硅酸盐、碳酸盐、硫酸岩等一系列矿物的组合,吸收特征与元素含量的关系非常复杂,混合矿物光谱是线性和非线性的单矿物组合,而且这些组合矿物的吸收特征,含量较少的矿物的吸收特征可能被含量较多的那个矿物的吸收特征所掩盖。因此需要建立的参数特征也非常复杂,这需要来进行回归分析。
从上面可以看出遥感光谱与地球化学具有一定正相关的联系,可以通过遥感图像来反演地球化学。
现有技术主要是针对遥感来反演地球化学或者是利用高分辨率的波段对低分辨率的波段进行融合,不能真实反映地球化学数据的本身,是遥感所揭示的地球化学,而遥感数据本身的融合也不能反映地球化学,两者处于脱节。主要表现在(1)地球化学没有遥感数据本身所反映的信息量大;(2)遥感数据没有地球化学数据的精确。如何实现地球化学和遥感数据的融合,使融合后的图像既能够反映地球化学的准确性,也具有遥感数据的图像细腻性,成为一个亟待解决的技术问题。
发明内容
本发明的目的是提供一种遥感图像与地球化学的融合方法及系统,以实现地球化学和遥感数据的融合,使融合后的图像既能够反映地球化学的准确性,也具有遥感数据的图像细腻性。
为实现上述目的,本发明提供了如下方案:
一种遥感图像与地球化学的融合方法,所述融合方法包括如下步骤:
获取遥感图像和所述遥感图像对应的地球化学图像;
对所述遥感图像进行最小噪声变换,获得多维图像数据;
对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;
将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;
利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;
对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像。
可选的,所述对所述遥感图像进行最小噪声变换,获得多维图像数据,具体包括:
根据所述遥感图像,构建所述遥感图像的总协方差矩阵;
采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵;
根据所述协方差矩阵构建变换矩阵;
利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据。
可选的,所述采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵,具体包括:
利用高通滤波器模块对所述遥感图像进行处理,获得噪声方差估计;
可选的,所述根据所述协方差矩阵构建变换矩阵,具体包括:
将所述协方差矩阵对角化,获得对角化矩阵;
根据所述协方差矩阵和所述对角化矩阵,获得变换矩阵。
可选额,所述利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据,具体包括:
利用变换矩阵对所述总协方差矩阵进行变换,获得降噪后的总协方差矩阵;
计算所述降噪后的总协方差矩阵的特征向量矩阵,获得多维图像数据。
可选的,所述对所述地球化学图像进行点面转化,获得地球化学栅格图像数据,具体包括:
建立点面转化坐标系;
列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值,形成点矢量;
计算所述点矢量的极差值;
判断所述极差值否大于预设阈值,若是,则采用多重分形法对所述点矢量进行插值,获得地球化学栅格图像数据,若否,则采用克里格法对所述点矢量进行插值,获得地球化学栅格图像数据。
可选的,所述将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据,具体包括:
利用公式计算所述多维图像数据和所述地球化学栅格图像数据的相似性系数;其中,SSIM(X,Yi)表示地球化学栅格化图像数据X,表示多维图像数据Y中的第i维图像数据Yi的相似系数;μX表示地球化学栅格化图像数据X的均值,μYi表示第i维图像数据 Yi的均值,σX表示地球化学栅格化图像数据X的方差,σYi表示第i维图像数据 Yi的方差,σXYi是地球化学栅格化图像数据X和第i维图像数据Yi的协方差,
m表示行数,n表示列数,C1,C2,C3为常数;
从所述多维图像数据中选取与所述地球化学栅格图像数据的相似性系数最大的图像数据作为相似图像数据。
一种遥感图像与地球化学的融合系统,所述融合系统包括:
图像获取模块,用于获取遥感图像和所述遥感图像对应的地球化学图像;
最小噪声变换模块,用于对所述遥感图像进行最小噪声变换,获得多维图像数据;
点面转化模块,用于对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;
匹配模块,用于将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;
替换模块,用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;
最小噪声反变换模块,用于对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像。
可选的,所述最小噪声变换模块,具体包括:
总协方差矩阵构建子模块,用于根据所述遥感图像,构建所述遥感图像的总协方差矩阵;
噪声协方差矩阵构建子模块,用于采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵;
变换矩阵构建子模块,用于根据所述协方差矩阵构建变换矩阵;
变换子模块,用于利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据。
可选的,所述噪声协方差矩阵构建子模块,具体包括:
滤波单元,用于利用高通滤波器模块对所述遥感图像进行处理,获得噪声方差估计;
根据本发明提供的具体实施例,本发明公开了以下技术效果:
本发明公开了一种遥感图像与地球化学的融合方法及系统。方法包括:首先,对所述遥感图像进行最小噪声变换,获得多维图像数据;对所述地球化学图像进行点面转化,获得地球化学栅格图像数据,即使地球化学图像变成栅格数据;然后将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;并利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据,最后对替换后的多维图像数据进行最小噪声反变换,得到融合后的遥感化学图像,实现了地球化学和遥感数据的融合,使融合后的图像既能够反映地球化学的准确性,也具有遥感数据的图像细腻性,对于地质找矿、土壤质量调查、农作物调查等具有重要意义。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为白云母基性及酸性成分与主要吸收峰的偏移图;
图2为白云母吸收谷与铝含量的关系图;
图3为白云母吸收谷位置与Al含量相关关系图;
图4为橄榄石吸收谷位置与Fe含量相关关系图;
图5为本发明提供的一种遥感图像与地球化学的融合方法的流程图;
图6为本发明选取的遥感图像发热频域直方图;
图7为本发明获得的遥感图像;
图8为本发明获取的地球化学图像;
图9为本发明的地球化学栅格化图像;
图10为采用本发明的方法融合后的图像;
图11为本发明提供的一种遥感图像与地球化学的融合系统的组成框图;
图12为本发明提供的一种遥感图像与地球化学的融合系统的结构图。
具体实施方式
本发明的目的是提供一种遥感图像与地球化学的融合方法及系统,以实现地球化学和遥感数据的融合,使融合后的图像既能够反映地球化学的准确性,也具有遥感数据的图像细腻性。
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对发明作进一步详细的说明。
如图5所示,本发明提供了一种遥感图像与地球化学的融合方法,所述融合方法包括如下步骤:
步骤501,获取遥感图像和所述遥感图像对应的地球化学图像;具体包括:
1.遥感图像数据的选择和预处理
遥感数据需要选择无云少云的秋冬季,尽量避免植被的影响,同时要求不是雨后的数据,由于水对电磁波的吸收,使得信息反射极微弱,冬季数据要避免选择冰雪覆盖数据,冰雪对电磁波的强反射,使得信息透射极微弱。另外在波段选择上根据获取隐伏目标的不同,选择不同的波段组合,如岩石需要选择 TM(Thematic Mapper)遥感图像数据的第五和第七波段加上第一和第四波段的组合。虽然波段组合不同,但一般须有红外数据的参与。为了相对精确的获得辅助尖端光谱,需要对图像进行去除高亮、极暗地物,所谓高亮地物指的是诸如冰雪、云等集中高反射的地物,所谓极暗地物指的是水体、阴影等地物,处理后的图像在直方图上近似正态分布。
其中,TM遥感图像数据是指美国陆地卫星4~5号专题制图仪(thematic mapper)所获取的多波段扫描影像。有7个波段,其波谱范围:TM-1为0.45~ 0.52微米,TM-2为0.52~0.60微米,TM-3为0.63~0.69微米,以上为可见光波段;TM-4为0.76~0.90微米,为近红外波段;TM-5为1.55~1.75微米,TM-7为2.08~2.35微米,为中红外波段;TM-6为10.40~12.50微米,为热红外波段。影像空间分辨率除热红外波段为120米外,其余均为30米,像幅185×185公里2。每波段像元数达61662个(TM-6为15422个)。一景 TM影像总信息量为230兆字节),约相当于MSS影像的7倍。因TM影像具较高空间分辨率、波谱分辨率、极为丰富的信息量和较高定位精度,成为 20世纪80年代中后期得到世界各国广泛应用的重要的地球资源与环境遥感数据源。能满足有关农、林、水、土、地质、地理、测绘、区域规划、环境监测等专题分析和编制1∶10万或更大比例尺专题图,修测中小比例尺地图的要求。 Landsat-7,星上携带专题制图仪ETM,ETM具有8个波段,其中1-5波段和 7波段是多光谱波段,空间分辨率是30米,第六波段是热红外波段,空间分辨率是120米,第8波段为全色波段,分辨率为15米。景宽185公里,景面积为34225平方公里。此星一直在正常运转,2000年遥感卫星地面站开始接收数据。
在一个M×N大小的遥感图像中,遥感图像某一个波段像元xj,k (j=1,m;k=1,n)区间为[x0,xn],统计这个窗口内的直方图,公式如下:
其中i∈[0,n],xj,k=xi为逻辑操作。
偏度系数满足公式:
其中ε1为给定一个很小的正数。
峰度系数满足公式:
其中ε2为给定一个很小的正数。
如图6所示本发明选取的遥感图像的频域直方图,其偏度系数和峰度系数满足上述要求,那么数据选择较好。然后对选择的遥感图像进行处理获得反射率或者辐射率。
直接获取的遥感图像存在几何变形,遥感器增益和偏移参数等,经过预处理得到带有坐标信息的行星反射率图像。具体公式如下:
其中:R反为行星反射率;pi=3.14;d表示影像当天的日地距离;A是太阳高度角;Esun值为大气层外相应波长的太阳光谱辐照度;L是辐亮度,可通过下式求出:
L=gain×DN+bias
其中:gain为增益,bias为偏移。
对于热红外数据需要温度和比发射率的分离。公式如下:
其中:T为温度,λ为波长,ε为比发射率,c1,c2为常数, c1=3.74818×10-4Wμm2,c2=1.43878×104Kμm。R为光谱辐射亮度,可以通过下面公式计算出:
QCAL是数据的实际辐射,LMINλ是QCAL=0时的光谱辐射值,LMAXλ是在QCAL=QCALMAX的光谱辐射值,QCALMAX是数据的图像辐射值。R 的单位是W/(m2×sr×μm)。
处理后的数据反射率数据R反和比发射率数据ε暂定为”选定数据”。
2.基础数据的生成
基础数据是噪声变换后通过去除干扰信息形成的数据,避免干扰信息影响结果,比如一些游离在空间数据云集边缘的云、阴影等干扰,这会或多或少影响处理的精度。常见的去除干扰模式有比值法、高端或低端切割法、Q值法、光谱角法等。
(1)高端或低端切割法
主要是利用干扰地物在遥感图像上某个波段有特征的高反射或强吸收,即某波段干扰地物有高值或低值,比如水体在TM/ETM的第7波段有低值,采用低端切割方法处理,而云在TM/ETM的第1波段有高值,采用高端切割方法处理,白泥地在TM/ETM的第3波段有高值,采用高端切割方法处理等。公式如下:
其中,i=0,…,n,n指所使用的遥感图像波段总数,xi和yi指分别指i波段去除干扰信息前和后的的波段值,b∈[1,…,n],Cb是常数,xb是原始b波段对应的值。这个公式的目的就是:给定一个约束的条件,使得这个条件大于或小于某一个数值的图像保留下来,其他的全部被赋值为零。
(2)比值法
常用比值法去除阴影、水体、冰雪、白泥地等多种干扰。首先判断干扰地物的各个波段的波谱特征,比如TM/ETM图像的阴影区第1波段明显大于第 7波段,因此采用第7波段比第1波段的方法,设定一个阈值进行去除,植被采用第5波段比第4波段或者第3波段比第4波段的方法等。公式如下:
其中,i=0,…,n,n指所使用的遥感图像波段总数,xi和yi指分别指i波段去除“尖锐”信息前和后的的波段值,a∈[1,…,n],Ca是常数,xa,xb是原始a、 b波段对应的值。这个公式的目的就是:给定一个约束的条件,使得这个条件大于或小于某一个数值的图像保留下来,其他的全部被赋值为零。
(3)Q值法
主要解决了雪边或湖边湿地、干河道、冲积区、薄云等干扰。定义Q值如下:
Q=(xa×ka-xb×kb)/xc×kc
其中,xa,xb,xc为参与主成分分析的波段,ka,kb,kc为参与主成分变化的xa,xb,xc对应本征向量的值。
去除干扰后的数据形成基础数据。
3.地球化学点数据整理
把地球化学数据进行处理,坐标与遥感图像的坐标相对应,按照表3 的方式把地球化学点数据进行输入,分别为ID,X,Y,地球化学元素值 Cu,其中ID为地球化学点数据的索引,x,y为对应于遥感图像坐标的坐标值,Cu为地球化学数据值。
表3
步骤502,对所述遥感图像进行最小噪声变换,获得多维图像数据;具体包括:根据所述遥感图像,构建所述遥感图像的总协方差矩阵Cz;采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵 Cn;根据所述协方差矩阵构建变换矩阵P;利用所述变换矩阵P对所述总协方差矩阵Cz进行变换,获得多维图像数据Y。
所述利用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵Cn,具体包括:
利用高通滤波器模块对所述遥感图像进行处理,获得噪声方差估计;根据所述噪声方差估计利用公式构建噪声协方差矩阵;其中,Cn为噪声协方差矩阵,是噪声方差估计Bi的倒数。具体为,首先通过对所述遥感图像的协方差Kn×n进行分解:
Kn×n=DkEkDk
上式中,ρkl为Kn×n在(k,l)处的相关系数(k≠l)。
所述根据所述协方差矩阵构建变换矩阵,具体包括:
将所述协方差矩阵对角化,获得对角化矩阵;根据所述协方差矩阵和所述对角化矩阵,获得变换矩阵。具体为:将Cn对角化获得对角化矩阵Dn:
Dn=UTCnU
PTCnP=I
其中I为单位矩阵。
所述利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据,具体包括:
利用变换矩阵对所述总协方差矩阵进行变换,获得降噪后的总协方差矩阵;计算所述降噪后的总协方差矩阵的特征向量矩阵,获得多维图像数据。具体的,对选定数据A可以通过Y=PA变换到新空间。
利用矩阵P对遥感数据总协方差矩阵Cz进行变换,得到噪声调整后的总协方差矩阵Cz-a;
Cz-a=PTCzP
计算协方差矩阵Cz-a的特征向量矩阵V,使得VTCz-aV=Dz-a,Dz-a为特征向量矩阵V所对应的特征值按降序排列得对角矩阵,且VTV=I;
其中,I为单位矩阵。
其中,Tmnf=PV构成最小噪声变换矩阵-MNF(Minimum Noise Fraction) 变换矩阵。所述MNF变换矩阵是一个最小噪声变换矩阵,旨在降噪.。最小噪声分离变换(MinimumNoise Fraction Rotation,MNF Rotation)工具用于判定图像数据内在的维数(即波段数),分离数据中的噪声,减少随后处理中的计算需求量。MNF本质上是两次层叠的主成分变换。第一次变换(基于估计的噪声协方差矩阵)用于分离和重新调节数据中的噪声,这步操作使变换后的噪声数据只有最小的方差且没有波段间的相关。第二步是对噪声白化数据(Noise-whitened)的标准主成分变换。为了进一步进行波谱处理,通过检查最终特征值和相关图像来判定数据的内在维数。数据空间可被分为两部分:一部分与较大特征值和相对应的特征图像相关,其余部分与近似相同的特征值以及噪声占主导地位的图像相关。
步骤503,对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;具体包括:建立点面转化坐标系;列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值,形成点矢量;计算所述点矢量的极差值;判断所述极差值否大于预设阈值,若是,则采用多重分形法对所述点矢量进行插值,获得地球化学栅格图像数据,若否,则采用克里格法对所述点矢量进行插值,获得地球化学栅格图像数据。具体的,按照坐标对点数据进行点面转化,也就是进行点数据的栅格化,根据坐标就需要对点数据进行插值,主要插值方法有多重分形和克里格方法:
①列出地球化学图像(温度图像)上每个点的坐标和化学参数值(温度值),形成矢量图层,即为点矢量;P(T)i,j=(xi,j,yi,j,Ti,j);
其中,xi,j为图像第i行,第j列的坐标值,Ti,j为图像第i行,第j列的计算温度值。1<i<m,1<j<n,m和n为地球化学图像行列最大值。
②计算矢量点的极差值|R(T)|;
如果|R(T)|>ε,采用多重分形方法对点矢量进行插值,如果|R(T)|<ε,采用克里格方法对点矢量进行插值。
a、多重分形
把温度图像每一行(或每一列)看作为一个空间序列(设m行n列),求出每个序列的Ti个极差和标准差,建立R(Ti)/S(Ti)、Ti的数据对,若的散点分布在一条直线附近,则该直线的斜率为赫斯特指数ai,也就是分形维度。
(1)先计算每一行(列)的维度,步骤如下:
Tj为第i行第j个数据的温度值,n为最大的列数。
Tj为第i行第j个数据的温度值,n为最大的列数。
ai为行分维,Tj为第i行第j个数据的温度值,n为最大的列数。
同理可以计算列分维bj。
(2)建立分维矩阵,元素为di,j=(ai,bj)i=1,2,...m;j=1,2,...n。
(4)计算每个单位的分维向量距离:
对每一点Ti,j的距离进行计算Δdij,利用Δdij对点数据进行插值。
b、克里格处理方法
克里格方法是在有限区域内对区域化变量进行无偏最优估计的一种方法。
根据前面计算得出的温度值图像T,计算某一点Ti,j附近n个温度值表示为 Ti,j(xk),k=1,…,n。
(1)计算期望与协方差
其中,(Ti,j(xk)-m(x))T为(Ti,j(xk)-m(x))的转置向量。
(3)计算最优性
(4)计算克里格方程;
最终获得地球化学栅格化图像数据。
步骤504,将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;具体包括:对地球化学处理地球化学栅格化图像数据和最小噪声变换生成的多通道(多维)图像数据(矩阵)进行匹配处理。在这里主要对比结构的相似性,因此必须参考图像的亮度,对比度和结构三个方面进行度量。假如经过p个波段经过MNF变换后的多维图像数据为同样地球化学栅格化图像数据为X;利用公式计算所述多维图像数据和所述地球化学栅格图像数据的相似性系数;其中, SSIM(X,Yi)表示地球化学栅格化图像数据X,表示多维图像数据Y中的第i维图像数据Yi的相似系数;μX表示地球化学栅格化图像数据X的均值,μYi表示第i维图像数据Yi的均值,σX表示地球化学栅格化图像数据X的方差,σYi表示第i维图像数据Yi的方差,σXYi是地球化学栅格化图像数据X和第i维图像数据Yi的协方差,
m表示行数,n表示列数,C1,C2,C3为常数,具体的,通常取 C1=(K1×L)2,C2=(K2×L)2,k1可以取值0.01,K2取值0.03,L可以取值 255。从所述多维图像数据中选取与所述地球化学栅格图像数据的相似性系数最大的图像数据作为相似图像数据。
步骤505,利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;
步骤506,对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像。
步骤506所述对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像,之后还包括:
步骤507,对所述融合后的遥感化学图像进行图像优化分析:
定义y为随机变量(地球化学点数据X),x1,x2,…,xn为n个自变量(对应的变化后的数据YZ),共观测m次,首先假设y与n个自变量之间存在线性关系:
y=a0+a1×x1+a2x2+…+anxn+ε
式中,a0,a1,a2,…,an为回归系数,是常数,表示在其它自变量不变情况下,xj(j=1,2,…n)增加或减少一个单位时的平均变化量,ε为去除n个自变量对y 的影响后的随机误差,上式称为多元线性回归模型。多元线性回归其条件是(1) y与x1,x2,…,xn之间具有线性关系;(2)每个观测值yj(j=1,2,…,m)相互独立; (3)ε服从正态分布。
首先用a0+a1×x1+a2x2+…+anxn来估计y的均值E(y),假定ε服从均值为0,方差为σ2的正态分布,即ε~N(0,σ2),则y服从均值为E(y),方差为的σ2的正态分布,即y~N[E(y),σ2],那么m组样本观测数据:
x11,x12,…,x1n,y1
x21,x22,…,x2n,y2
…………………
xm1,xm2,…,xmn,ym
式中,xij表示xj在第i次的观测值。有如下公式:
上式为n元线性回归的数学模型,式中,a0,a1,a2,…,an为n+1个待定参数,ε1,ε2,…,εm为m个相互独立并且服从同一正态分布的随机变量。为了简化表示,利用矩阵形式:
则n元线性回归的数学模型为
Y=AX+E
根据公式来进行最小二乘估计,首先假设b0,b1,b2,…,bn分别为n+1个回归系数a0,a1,a2,…,an的最小二乘估计值,于是观测值如下表示:
yj=b0xj1+b1xj2+…+bnxjn+ej
达到最小,根据高等数学原理,极值在函数微分值为0处,确立方程,
从上述公式中得到正规方程,
根据矩阵X,系数两边等式用C和D表示,则,
那么,正规方程的矩阵形式为:CB=(X′X)B=X′y=D。
B为未知向量,如果矩阵系数C满秩,其逆矩阵存在,可以反解出未知向量B,B=A-1D=(X′X)X′y。
步骤508,图像数据检验;
回归方程的假设检验和评价一般采用离差分析法。定义总变异,
总变异自由度为m-1。
如果观测值给定,总变异确定,可以用SS和MS来衡量回归效果,SS越大回归效果越显著,MS越大回归效果不好。
为了检验总的回归效果,定义了无量纲指标-决定系数R2来表示,
R2反映回归离差对总变异的贡献比例。R=R1/2称为复相关系数,反映全部自变量与因变量的相关程度。R2和R值越大,回归效果越好。
上述位总体回归效果检验,不能说明每个自变量x1,x2,…,xn对因变量y 都重要,有些自变量对因变量可能不起作用,或者作用为其他自变量所代替,这就需要把这些自变量从回归方程中剔除,建议每个自变量xi是否显著,假设H0:ai=0,i=1,2,…n。
(1)F值检验
在H0:ai=0假设下,
对给定置信度α,从F值分布表中查与β对应的临界值Fβ,如果|Fi|〉F β,拒绝假设H0,认为n个自变量总体回归效果显著,反之,总体回归效果不显著。
(2)t检验
对给定检验水平β,从t值分布表中查与β对应的临界值tβ,如果|ti|〉t β,拒绝假设H0,认为ai与0值有显著差异,不应剔除,反之,应剔除。
(3)p值检验
假设H0:ai=0,服从自由度分别为1与m-n-1的p分布统计量,
对给定检验水平β,从p值分布表中可以查到临界值pβ(1,m-n-1),如果pi〉pβ(1,m-n-1),拒绝假设H0,认为xi对y值有重要作用,不应剔除,反之,应剔除。
步骤509,图像输出;
经过以上基本处理后,能够得出的融合后的遥感图像。选择最优层,进行组合式中,Si为为第i个波段的标准差,Ri,j为第i,j波段的相关系数。OIF越大,表明包含的信息量越大,因此,最大的OIF波段组合为最佳波段组合,采用最佳波段组合合成的图像作为底图。组合后的图像利用RGB合成一幅假彩色图像,这样做出来的图像会反映更多的去除背景噪声后的精细信息,更加适合人的视觉习惯,简单、易用。
输出一幅适合人眼观察的叠加图像。通过软件输出成为JPG或者TIF格式的最终图像。
具体的图像如图7-10所示,其中图7为遥感图像,图8为地球化学图像,图9为地球化学栅格化图像,图10为融合后的图像。
如图11和12所示本发明还提供了一种遥感图像与地球化学的融合系统,所述融合系统包括:
图像获取模块1101,用于获取遥感图像和所述遥感图像对应的地球化学图像。
最小噪声变换模块1102,用于对所述遥感图像进行最小噪声变换,获得多维图像数据;所述最小噪声变换模块1102,具体包括:总协方差矩阵构建子模块,用于根据所述遥感图像,构建所述遥感图像的总协方差矩阵;噪声协方差矩阵构建子模块,用于采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵;变换矩阵构建子模块,用于根据所述协方差矩阵构建变换矩阵;变换子模块,用于利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据。
所述噪声协方差矩阵构建子模块,具体包括:
滤波单元,用于利用高通滤波器模块对所述遥感图像进行处理,获得噪声方差估计;噪声协方差矩阵构建单元,用于根据所述噪声方差估计利用公式构建噪声协方差矩阵;其中,Cn为噪声协方差矩阵,是噪声方差估计Bi的倒数。
具体的,最小噪声变换模块1102中设置有最小噪声变换处理器。所述最小噪声变换处理器用于对遥感图像或遥感图像经预处理得到的基础数据。
点面转化模块1103,用于对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;具体的,所述点面转换模块1103中设置有点面转换装置,所述点面转换装置用于将所述地球化学图像转换成栅格数据,即地球化学栅格图像数据。
匹配模块1104,用于将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;所述匹配模块1104中设置有多通道匹配处理器,所述多通道匹配处理器用于进行多通道匹配处理。
替换模块1105,用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;所述替换模块1105中设置有替换合成装置,所述替换合成装置用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据。
最小噪声反变换模块1106,用于对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像。所述最小噪声反变换模块1106中设置有最小噪声反变换处理器,所述最小噪声反变换处理器用于对所述替换后的多维图像数据进行最小噪声反变换,获得融合后的遥感化学图像。
所述融合系统还包括图像优化分析装置、图像输出装置,所述图像优化分析装置用于对所述融合后的遥感化学图像进行参数优化;所述图像输出装置内设置有精度控制器和图像合成装置,所述精度控制器用于对融合化学图像进行检验,提高融合化学图像的精度,所述图像合成装置用于从融合化学图像中选择最优层并进行合成,得到RGB合成的假彩色图像,使图像反映更多的去除背景噪声后的精细信息,更加适合人的视觉习惯、简单、易用,并将最终图像输出。
本发明主要针对遥感图像与地球化学融合难的问题,采用空间变换方式,把水系或者土壤的特征表征出来,采用水系或者土壤地球化学数据点阵变换栅格后替换,从而实现遥感图像与地球化学高度融合,是两者之间的一种重要的补充,能够大幅度提升两种数据的使用性和效率。对于地质找矿、土壤质量调查、农作物调查等具有重要的意义。
本说明书中各个实施例采用递进的方式描述,每个实施例重点说明的都是与其他实施例的不同之处,各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言,由于其与实施例公开的方法相对应,所以描述的比较简单,相关之处参见方法部分说明即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想,所描述的实施例仅仅是本发明的一部分实施例,而不是全部的实施例,基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
Claims (9)
1.一种遥感图像与地球化学的融合方法,其特征在于,所述融合方法包括如下步骤:
获取遥感图像和所述遥感图像对应的地球化学图像;
对所述遥感图像进行最小噪声分离变换,获得多维图像数据;
对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;具体包括:建立点面转化坐标系;列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值,形成点矢量;计算所述点矢量的极差值;判断所述极差值否大于预设阈值,若是,则采用多重分形法对所述点矢量进行插值,获得地球化学栅格图像数据,若否,则采用克里格法对所述点矢量进行插值,获得地球化学栅格图像数据;
将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;
利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;
对所述替换后的多维图像数据进行最小噪声分离反变换,获得融合后的遥感化学图像。
2.根据权利要求1所述的一种遥感图像与地球化学的融合方法,其特征在于,所述对所述遥感图像进行最小噪声分离变换,获得多维图像数据,具体包括:
根据所述遥感图像,构建所述遥感图像的总协方差矩阵;
采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵;
根据所述协方差矩阵构建变换矩阵;
利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据。
4.根据权利要求2所述的一种遥感图像与地球化学的融合方法,其特征在于,所述根据所述协方差矩阵构建变换矩阵,具体包括:
将所述协方差矩阵对角化,获得对角化矩阵;
根据所述协方差矩阵和所述对角化矩阵,获得所述变换矩阵。
5.根据权利要求2所述的一种遥感图像与地球化学的融合方法,其特征在于,所述利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据,具体包括:
利用所述变换矩阵对所述总协方差矩阵进行变换,获得降噪后的总协方差矩阵;
计算所述降噪后的总协方差矩阵的特征向量矩阵,获得多维图像数据。
6.根据权利要求1所述的一种遥感图像与地球化学的融合方法,其特征在于,所述将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据,具体包括:
利用公式计算所述多维图像数据和所述地球化学栅格图像数据的相似性系数;其中,SSIM(X,Yi)表示地球化学栅格化图像数据X与多维图像数据Y中的第i维图像数据Yi的相似系数;μX表示地球化学栅格化图像数据X的均值,μYi表示第i维图像数据Yi的均值,σX表示地球化学栅格化图像数据X的方差,σYi表示第i维图像数据Yi的方差,σXYi是地球化学栅格化图像数据X和第i维图像数据Yi的协方差,
m表示行数,n表示列数,C1,C2,C3为常数;
从所述多维图像数据中选取与所述地球化学栅格图像数据的相似性系数最大的图像数据作为相似图像数据。
7.一种遥感图像与地球化学的融合系统,其特征在于,所述融合系统包括:
图像获取模块,用于获取遥感图像和所述遥感图像对应的地球化学图像;
最小噪声变换模块,用于对所述遥感图像进行最小噪声分离变换,获得多维图像数据;
点面转化模块,用于对所述地球化学图像进行点面转化,获得地球化学栅格图像数据;具体包括:建立点面转化坐标系;列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值,形成点矢量;计算所述点矢量的极差值;判断所述极差值否大于预设阈值,若是,则采用多重分形法对所述点矢量进行插值,获得地球化学栅格图像数据,若否,则采用克里格法对所述点矢量进行插值,获得地球化学栅格图像数据;
匹配模块,用于将所述多维图像数据和所述地球化学栅格图像数据进行匹配,获得相似图像数据;
替换模块,用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据,获得替换后的多维图像数据;
最小噪声反变换模块,用于对所述替换后的多维图像数据进行最小噪声分离反变换,获得融合后的遥感化学图像。
8.根据权利要求7所述的一种遥感图像与地球化学的融合系统,其特征在于,所述最小噪声变换模块,具体包括:
总协方差矩阵构建子模块,用于根据所述遥感图像,构建所述遥感图像的总协方差矩阵;
噪声协方差矩阵构建子模块,用于采用高通滤波器模板对所述遥感图像进行处理,获得所述遥感图像的噪声协方差矩阵;
变换矩阵构建子模块,用于根据所述协方差矩阵构建变换矩阵;
变换子模块,用于利用所述变换矩阵对所述总协方差矩阵进行变换,获得多维图像数据。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810177878.0A CN108427964B (zh) | 2018-03-05 | 2018-03-05 | 一种遥感图像与地球化学的融合方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810177878.0A CN108427964B (zh) | 2018-03-05 | 2018-03-05 | 一种遥感图像与地球化学的融合方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108427964A CN108427964A (zh) | 2018-08-21 |
CN108427964B true CN108427964B (zh) | 2020-05-12 |
Family
ID=63157659
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810177878.0A Active CN108427964B (zh) | 2018-03-05 | 2018-03-05 | 一种遥感图像与地球化学的融合方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108427964B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109145881B (zh) * | 2018-10-09 | 2021-06-22 | 中国地质科学院矿产资源研究所 | 一种遥感图像膏盐信息提取方法及装置 |
CN109345459B (zh) * | 2018-10-30 | 2020-08-28 | 中国科学院遥感与数字地球研究所 | 一种提升地球化学元素图层分辨率的方法及系统 |
CN109902557A (zh) * | 2019-01-14 | 2019-06-18 | 云南大学 | 蚀变特征矿物的蚀变异常信息的提取方法 |
CN109740576B (zh) * | 2019-01-31 | 2023-01-10 | 交通运输部天津水运工程科学研究所 | 一种监管公路临时占地利用与恢复的卫星遥感方法 |
CN110263872B (zh) * | 2019-06-26 | 2022-05-17 | 上海鹰瞳医疗科技有限公司 | 训练数据处理方法及装置 |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP1875041A2 (en) * | 2005-04-05 | 2008-01-09 | Raymond P. Murphy | Well bore fluid redistribution and fluid disposal |
CN102721650B (zh) * | 2012-06-13 | 2014-05-14 | 中国地质科学院矿产资源研究所 | 基于特征指标的矿物成分遥感信息提取方法及装置 |
CN103383348B (zh) * | 2013-05-28 | 2015-09-30 | 吉林大学 | 植被覆盖区高光谱遥感蚀变矿物提取方法 |
CN104123559B (zh) * | 2014-06-30 | 2017-07-07 | 中国地质科学院矿产资源研究所 | 盐湖区地下含钾卤水资源的多源遥感判别方法和系统 |
-
2018
- 2018-03-05 CN CN201810177878.0A patent/CN108427964B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN108427964A (zh) | 2018-08-21 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108427964B (zh) | 一种遥感图像与地球化学的融合方法及系统 | |
Chen et al. | SAR and multispectral image fusion using generalized IHS transform based on à trous wavelet and EMD decompositions | |
Bourdin et al. | Temperature structure of the intergalactic medium within seven nearby and bright clusters of galaxies observed with XMM-Newton | |
CN108765359A (zh) | 一种基于jskf模型和nsct技术的高光谱遥感影像与全色图像的融合方法 | |
CN111337434A (zh) | 一种矿区复垦植被生物量估算方法及系统 | |
CN111986162B (zh) | 基于粗定位和协同表示的高光谱异常点快速检测方法 | |
CN109829872B (zh) | 一种用于内陆水体遥感的多时相多源遥感影像融合方法 | |
CN113970376A (zh) | 一种基于海洋区域再分析资料的卫星红外载荷定标方法 | |
CN111008664A (zh) | 一种基于空谱联合特征的高光谱海冰检测方法 | |
Al-Wassai et al. | Multisensor images fusion based on feature-level | |
CN116309136A (zh) | 一种基于sar先验知识引导的遥感影像云区重建方法 | |
CN104951800A (zh) | 一种面向资源开采型地区的遥感影像融合方法 | |
CN117115669B (zh) | 双条件质量约束的对象级地物样本自适应生成方法及系统 | |
CN113744134A (zh) | 基于光谱解混卷积神经网络的高光谱图像超分辨方法 | |
Meng et al. | Improving the spatial resolution of hyperspectral image using panchromatic and multispectral images: An integrated method | |
Muhammad et al. | Evaluation of wavelet transform algorithms for multi-resolution image fusion | |
CN109696406B (zh) | 一种基于复合端元的月表高光谱图像阴影区域解混方法 | |
CN111383203B (zh) | 基于分区域拟合的全色与多光谱遥感图像融合方法 | |
Gross et al. | Application of spatial resolution enhancement and spectral mixture analysis to hyperspectral images | |
Krawczyk et al. | Investigation of interpretation possibilities of spectral high-dimensional measurements by means of principal component analysis: a concept for physical interpretation of those measurements | |
Zhang et al. | A spatial-temporal-spectral blending model using satellite images | |
Smara et al. | Multisource ERS-1 and optical data for vegetal cover assessment and monitoring in a semi-arid region of Algeria | |
Yang et al. | Practical image fusion method based on spectral mixture analysis | |
Kim et al. | Spatiotemporal fusion of high resolution land surface temperature using thermal sharpened images from regression-based urban indices | |
He | Quantitative hyperspectral imaging pipeline to recover surface images from CRISM radiance data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |