一种遥感图像与地球化学的融合方法及系统

摘要

本发明公开了一种遥感图像与地球化学的融合方法及系统。方法包括：首先，对所述遥感图像进行最小噪声变换，获得多维图像数据；对所述地球化学图像进行点面转化，获得地球化学栅格图像数据，即使地球化学图像变成栅格数据；然后将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；并利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据，最后对替换后的多维图像数据进行最小噪声反变换，得到融合后的遥感化学图像，实现了地球化学和遥感数据的融合。

说明书

技术领域

本发明涉及图像处理领域，特别涉及一种遥感图像与地球化学的融合方法及系统。

背景技术

遥感图像是栅格数据，而地球化学是点阵数据，通常遥感图像和地球化学的融合不能够很好的实现，仅仅是形式的融合，比如常规的主分类融合方法、彩色变换融合方法，融合的图像可识别程度很低。

遥感地球化学是利用遥感的手段反演元素丰度特征，因矿物的吸收特征主要是分子键振动形成的吸收特征或者是电子跃迁形成的吸收特征。主要的吸收特征(主要指吸收谷的波长)与分子振动和电子跃迁有关，这与分子晶体结构有关，如羟基、碳酸盐基团的吸收特征，和铁离子、锰离子的吸收特征，而与元素丰度关系不大。但是元素含量对分子吸收特征是有影响的，如高铝白云母和低铝白云母，其吸收特征有一定的差异性，铝含量的增加吸收谷向左偏移 (Swayze,1997)。在月球上的铁、钛元素含量与月壤光谱特征有关系，Lucey 给出了经验公式(Lucey et al.,1998)。这就说明元素含量对光谱有一定的影响。

由于可见光-短波红外的光谱都是矿物(或岩石)的分子振动和电子跃迁产生的。阳离子的大小和含量对矿物波谱特征有影响，随着矿物(或岩石)成分从酸性向基性过渡，它们的主要的吸收谷向在短波红外区域向波长较长的区域偏移。例如钠-白云母、钾-白云母，镁硅白云母、镁铁-多硅白云母，吸收峰依次在2115nm；2205nm；2220nm和2230nm。

参见附图1：为了分析元素对光谱特征的影响，从USGS最新谱库中获得白云母的光谱特征和元素含量相对特征，采用白云母AL含量的不同来进行分析，单一的蚀变矿物元素含量不同对光谱特征影响的区别如表1所示。参见附图2中可以看出，在热红外区，AL含量虽对波谱吸收特征有影响，但许多因素使得吸收特征变得复杂，编号为G108和G111白云母吸收特征就不明显。因此主要是利用短波红外区域进行区分。

表1

白云母编号AL2O3元素含量吸收特征G12035.691.412191，2.200405G11933.341.411807，2.201339G11831.311.413346，2.207902G11735.541.4112305，2.1999385G11631.061.414503，2.213559G11433.091.412961，2.205085G11335.641.411807，2.197607G11128.031.416435，2.2197195G10835.01.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左右吸收谷位置更加相关，那么铝含量与光谱吸收谷位置关系函数如下：

式中，是Al2O3的含量，pos_2.2是白云母在2.2μm左右处的吸收位置。

为了更好的表示元素含量与光谱吸收特征的关系，从USGS最新谱库上获得橄榄石的光谱特征与元素含量如表2所示，分析Fe的含量与橄榄石光谱吸收谷位置的关系，主要是2.2μm左右这个位置的特征，做一个二维数据表，然后利用最小二乘法进行拟合，得出橄榄石铁含量与吸收谷位置特征的关系公司，是典型的线性关系，铁含量与2.2μm左右吸收谷位置相关性为0.975，相关性非常好(图4)。铁含量与光谱吸收位置关系函数如下：

Con_FeO＝0.0003pos_2.2+2.1919

式中，Con_FeO是FeO的含量，pos_2.2是橄榄石在2.2μm左右处的吸收位置。

表2

橄榄石编号FeO的含量吸收特征KI300562.822.214504KI305430.592.202275KI318841.432.207432KI318934.632.204147KI329153.652.209784KI337759.752.212614KI414347.652.207902

对于岩石来说，由于岩石是矿物的组合，是硅酸盐、碳酸盐、硫酸岩等一系列矿物的组合，吸收特征与元素含量的关系非常复杂，混合矿物光谱是线性和非线性的单矿物组合，而且这些组合矿物的吸收特征，含量较少的矿物的吸收特征可能被含量较多的那个矿物的吸收特征所掩盖。因此需要建立的参数特征也非常复杂，这需要来进行回归分析。

从上面可以看出遥感光谱与地球化学具有一定正相关的联系，可以通过遥感图像来反演地球化学。

现有技术主要是针对遥感来反演地球化学或者是利用高分辨率的波段对低分辨率的波段进行融合，不能真实反映地球化学数据的本身，是遥感所揭示的地球化学，而遥感数据本身的融合也不能反映地球化学，两者处于脱节。主要表现在(1)地球化学没有遥感数据本身所反映的信息量大；(2)遥感数据没有地球化学数据的精确。如何实现地球化学和遥感数据的融合，使融合后的图像既能够反映地球化学的准确性，也具有遥感数据的图像细腻性，成为一个亟待解决的技术问题。

发明内容

本发明的目的是提供一种遥感图像与地球化学的融合方法及系统，以实现地球化学和遥感数据的融合，使融合后的图像既能够反映地球化学的准确性，也具有遥感数据的图像细腻性。

为实现上述目的，本发明提供了如下方案：

一种遥感图像与地球化学的融合方法，所述融合方法包括如下步骤：

获取遥感图像和所述遥感图像对应的地球化学图像；

对所述遥感图像进行最小噪声变换，获得多维图像数据；

对所述地球化学图像进行点面转化，获得地球化学栅格图像数据；

将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；

利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据；

对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像。

可选的，所述对所述遥感图像进行最小噪声变换，获得多维图像数据，具体包括：

根据所述遥感图像，构建所述遥感图像的总协方差矩阵；

采用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵；

根据所述协方差矩阵构建变换矩阵；

利用所述变换矩阵对所述总协方差矩阵进行变换，获得多维图像数据。

可选的，所述采用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵，具体包括：

利用高通滤波器模块对所述遥感图像进行处理，获得噪声方差估计；

根据所述噪声方差估计利用公式构建噪声协方差矩阵；其中，C_n为噪声协方差矩阵，是噪声方差估计B_i的倒数。

可选的，所述根据所述协方差矩阵构建变换矩阵，具体包括：

将所述协方差矩阵对角化，获得对角化矩阵；

根据所述协方差矩阵和所述对角化矩阵，获得变换矩阵。

可选额，所述利用所述变换矩阵对所述总协方差矩阵进行变换，获得多维图像数据，具体包括：

利用变换矩阵对所述总协方差矩阵进行变换，获得降噪后的总协方差矩阵；

计算所述降噪后的总协方差矩阵的特征向量矩阵，获得多维图像数据。

可选的，所述对所述地球化学图像进行点面转化，获得地球化学栅格图像数据，具体包括：

建立点面转化坐标系；

列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值，形成点矢量；

计算所述点矢量的极差值；

判断所述极差值否大于预设阈值，若是，则采用多重分形法对所述点矢量进行插值，获得地球化学栅格图像数据，若否，则采用克里格法对所述点矢量进行插值，获得地球化学栅格图像数据。

可选的，所述将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据，具体包括：

利用公式计算所述多维图像数据和所述地球化学栅格图像数据的相似性系数；其中，SSIM(X，Y_i)表示地球化学栅格化图像数据X，表示多维图像数据Y中的第i维图像数据Y_i的相似系数；μ_X表示地球化学栅格化图像数据X的均值，μ_Yi表示第i维图像数据>i的均值，σ_X表示地球化学栅格化图像数据X的方差，σ_Yi表示第i维图像数据>i的方差，σ_XYi是地球化学栅格化图像数据X和第i维图像数据Y_i的协方差,

m表示行数，n表示列数，C₁，C₂，C₃为常数；

从所述多维图像数据中选取与所述地球化学栅格图像数据的相似性系数最大的图像数据作为相似图像数据。

一种遥感图像与地球化学的融合系统，所述融合系统包括：

图像获取模块，用于获取遥感图像和所述遥感图像对应的地球化学图像；

最小噪声变换模块，用于对所述遥感图像进行最小噪声变换，获得多维图像数据；

点面转化模块，用于对所述地球化学图像进行点面转化，获得地球化学栅格图像数据；

匹配模块，用于将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；

替换模块，用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据；

最小噪声反变换模块，用于对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像。

可选的，所述最小噪声变换模块，具体包括：

总协方差矩阵构建子模块，用于根据所述遥感图像，构建所述遥感图像的总协方差矩阵；

噪声协方差矩阵构建子模块，用于采用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵；

变换矩阵构建子模块，用于根据所述协方差矩阵构建变换矩阵；

变换子模块，用于利用所述变换矩阵对所述总协方差矩阵进行变换，获得多维图像数据。

可选的，所述噪声协方差矩阵构建子模块，具体包括：

滤波单元，用于利用高通滤波器模块对所述遥感图像进行处理，获得噪声方差估计；

噪声协方差矩阵构建单元，用于根据所述噪声方差估计利用公式构建噪声协方差矩阵；其中，C_n为噪声协方差矩阵，是噪声方差估计B_i的倒数。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种遥感图像与地球化学的融合方法及系统。方法包括：首先，对所述遥感图像进行最小噪声变换，获得多维图像数据；对所述地球化学图像进行点面转化，获得地球化学栅格图像数据，即使地球化学图像变成栅格数据；然后将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；并利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据，最后对替换后的多维图像数据进行最小噪声反变换，得到融合后的遥感化学图像，实现了地球化学和遥感数据的融合，使融合后的图像既能够反映地球化学的准确性，也具有遥感数据的图像细腻性，对于地质找矿、土壤质量调查、农作物调查等具有重要意义。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图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大小的遥感图像中，遥感图像某一个波段像元x_j,k>0，x_n]，统计这个窗口内的直方图，公式如下：

其中i∈[0,n]，x_j,k＝x_i为逻辑操作。

这个遥感图像某一个波段像元值为x_j,k(j＝1,m；k＝1,n)，像元均值为标准差为σ，利用偏度系数和峰度系数进行判断直方图是否近似正态分布。

偏度系数满足公式：

其中ε₁为给定一个很小的正数。

峰度系数满足公式：

其中ε₂为给定一个很小的正数。

如图6所示本发明选取的遥感图像的频域直方图，其偏度系数和峰度系数满足上述要求，那么数据选择较好。然后对选择的遥感图像进行处理获得反射率或者辐射率。

直接获取的遥感图像存在几何变形，遥感器增益和偏移参数等，经过预处理得到带有坐标信息的行星反射率图像。具体公式如下：

其中：R_反为行星反射率；pi＝3.14；d表示影像当天的日地距离；A是太阳高度角；E_sun值为大气层外相应波长的太阳光谱辐照度；L是辐亮度，可通过下式求出：

L＝gain×DN+bias

其中：gain为增益，bias为偏移。

对于热红外数据需要温度和比发射率的分离。公式如下：

其中：T为温度，λ为波长，ε为比发射率，c₁，c₂为常数，>1＝3.74818×10^-4Wμm²，c₂＝1.43878×10⁴Kμm。R为光谱辐射亮度，可以通过下面公式计算出：

QCAL是数据的实际辐射，LMINλ是QCAL＝0时的光谱辐射值，LMAXλ是在QCAL＝QCALMAX的光谱辐射值，QCALMAX是数据的图像辐射值。R 的单位是W/(m²×sr×μm)。

处理后的数据反射率数据R_反和比发射率数据ε暂定为”选定数据”。

2.基础数据的生成

基础数据是噪声变换后通过去除干扰信息形成的数据，避免干扰信息影响结果，比如一些游离在空间数据云集边缘的云、阴影等干扰，这会或多或少影响处理的精度。常见的去除干扰模式有比值法、高端或低端切割法、Q值法、光谱角法等。

(1)高端或低端切割法

主要是利用干扰地物在遥感图像上某个波段有特征的高反射或强吸收，即某波段干扰地物有高值或低值，比如水体在TM/ETM的第7波段有低值，采用低端切割方法处理，而云在TM/ETM的第1波段有高值，采用高端切割方法处理，白泥地在TM/ETM的第3波段有高值，采用高端切割方法处理等。公式如下：

其中，i＝0,…,n，n指所使用的遥感图像波段总数，x_i和y_i指分别指i波段去除干扰信息前和后的的波段值，b∈[1,…,n]，C_b是常数，x_b是原始b波段对应的值。这个公式的目的就是：给定一个约束的条件，使得这个条件大于或小于某一个数值的图像保留下来，其他的全部被赋值为零。

(2)比值法

常用比值法去除阴影、水体、冰雪、白泥地等多种干扰。首先判断干扰地物的各个波段的波谱特征，比如TM/ETM图像的阴影区第1波段明显大于第 7波段，因此采用第7波段比第1波段的方法，设定一个阈值进行去除，植被采用第5波段比第4波段或者第3波段比第4波段的方法等。公式如下：

其中，i＝0,…,n，n指所使用的遥感图像波段总数，x_i和y_i指分别指i波段去除“尖锐”信息前和后的的波段值，a∈[1,…,n]，C_a是常数，x_a，x_b是原始a、>

(3)Q值法

主要解决了雪边或湖边湿地、干河道、冲积区、薄云等干扰。定义Q值如下：

Q＝(x_a×k^a-x_b×k^b)/x_c×k^c

其中，x_a，x_b，x_c为参与主成分分析的波段，k^a，k^b，k^c为参与主成分变化的x_a，x_b，x_c对应本征向量的值。

去除干扰后的数据形成基础数据。

3.地球化学点数据整理

把地球化学数据进行处理，坐标与遥感图像的坐标相对应，按照表3 的方式把地球化学点数据进行输入，分别为ID，X，Y，地球化学元素值 Cu，其中ID为地球化学点数据的索引，x，y为对应于遥感图像坐标的坐标值，Cu为地球化学数据值。

表3

步骤502，对所述遥感图像进行最小噪声变换，获得多维图像数据；具体包括：根据所述遥感图像，构建所述遥感图像的总协方差矩阵C_z；采用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵>n；根据所述协方差矩阵构建变换矩阵P；利用所述变换矩阵P对所述总协方差矩阵C_z进行变换，获得多维图像数据Y。

所述利用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵C_n，具体包括：

利用高通滤波器模块对所述遥感图像进行处理，获得噪声方差估计；根据所述噪声方差估计利用公式构建噪声协方差矩阵；其中，C_n为噪声协方差矩阵，是噪声方差估计B_i的倒数。具体为，首先通过对所述遥感图像的协方差K_n×n进行分解：

K_n×n＝D_kE_kD_k

D_k＝diag{σ₁，…，σ_n}为一个对角矩阵，n为波段数，K_n×n的对角元素为第i波段样本变量B_i的方差，且，

上式中，ρ_kl为K_n×n在(k,l)处的相关系数(k≠l)。

也可以分解成：

其中，为一个对角矩阵，为的对角元素，并且，

式中，ξ_kl为在(k,l)处的相关系数(k≠l),ξ_i和σ_i存在如下关系：

其中，为B_l与共n-1个波段的多元相关系数，可以通过多元回归理论获得，是B_i噪声方差估计的倒数。因此，噪声协方差矩阵C_n可以通过对角阵来估计。

所述根据所述协方差矩阵构建变换矩阵，具体包括：

将所述协方差矩阵对角化，获得对角化矩阵；根据所述协方差矩阵和所述对角化矩阵，获得变换矩阵。具体为：将C_n对角化获得对角化矩阵D_n：

D_n＝U^TC_nU

式中，D_n为C_n的特征值按照降序排列的对角矩阵，U为由特征向量组成的正交矩阵，在此基础上，获得变换矩阵且，

P^TC_nP＝I

其中I为单位矩阵。

所述利用所述变换矩阵对所述总协方差矩阵进行变换，获得多维图像数据，具体包括：

利用变换矩阵对所述总协方差矩阵进行变换，获得降噪后的总协方差矩阵；计算所述降噪后的总协方差矩阵的特征向量矩阵，获得多维图像数据。具体的，对选定数据A可以通过Y＝PA变换到新空间。

利用矩阵P对遥感数据总协方差矩阵C_z进行变换，得到噪声调整后的总协方差矩阵C_z-a；

C_z-a＝P^TC_zP

计算协方差矩阵C_z-a的特征向量矩阵V，使得V^TC_z-aV＝D_z-a,D_z-a为特征向量矩阵V所对应的特征值按降序排列得对角矩阵，且V^TV＝I；

其中，I为单位矩阵。

其中，T_mnf＝PV构成最小噪声变换矩阵-MNF(Minimum>

步骤503，对所述地球化学图像进行点面转化，获得地球化学栅格图像数据；具体包括：建立点面转化坐标系；列出所述地球化学图像上每个点在所述点面转化坐标系中的坐标和化学参数值，形成点矢量；计算所述点矢量的极差值；判断所述极差值否大于预设阈值，若是，则采用多重分形法对所述点矢量进行插值，获得地球化学栅格图像数据，若否，则采用克里格法对所述点矢量进行插值，获得地球化学栅格图像数据。具体的，按照坐标对点数据进行点面转化，也就是进行点数据的栅格化，根据坐标就需要对点数据进行插值，主要插值方法有多重分形和克里格方法：

①列出地球化学图像(温度图像)上每个点的坐标和化学参数值(温度值)，形成矢量图层，即为点矢量；P(T)_i，j＝(x_i，j，y_i，j，T_i，j)；

其中，x_i，j为图像第i行，第j列的坐标值，T_i，j为图像第i行，第j列的计算温度值。1<i<m，1<j<n，m和n为地球化学图像行列最大值。

②计算矢量点的极差值|R(T)|；

如果|R(T)|>ε，采用多重分形方法对点矢量进行插值，如果|R(T)|<ε，采用克里格方法对点矢量进行插值。

a、多重分形

把温度图像每一行(或每一列)看作为一个空间序列(设m行n列)，求出每个序列的T_i个极差和标准差，建立R(T_i)/S(T_i)、T_i的数据对，若的散点分布在一条直线附近，则该直线的斜率为赫斯特指数a_i，也就是分形维度。

(1)先计算每一行(列)的维度，步骤如下：

每一行(行号i)的极差：

T_j为第i行第j个数据的温度值,n为最大的列数。

每一行(行号i)的标准差：

Tj为第i行第j个数据的温度值,n为最大的列数。

建立极差与标准差散点关系，

a_i为行分维，T_j为第i行第j个数据的温度值，n为最大的列数。

同理可以计算列分维b_j。

(2)建立分维矩阵，元素为d_i，j＝(a_i，b_j)i＝1，2，...m；j＝1，2，...n。

(3)利用方程计算行总维度A。

其中，i 为行号。采用同样的方法计算列总维度B。

(4)计算每个单位的分维向量距离：

对每一点T_i，j的距离进行计算Δd_ij，利用Δd_ij对点数据进行插值。

b、克里格处理方法

克里格方法是在有限区域内对区域化变量进行无偏最优估计的一种方法。

根据前面计算得出的温度值图像T，计算某一点T_i，j附近n个温度值表示为>i，j(x_k)，k＝1，…，n。

(1)计算期望与协方差

数学期望向量

计算协方差

其中，(T_i，j(x_k)-m(x))^T为(T_i，j(x_k)-m(x))的转置向量。

(2)计算无偏性

其中λ_k为待求值。

(3)计算最优性

在满足无偏性条件下，获得方差

构建拉格朗日方程

其中，μ和λ_k为待求值，分别对F求μ和λ_k的导数并令其为0。也就是可以得到克里格方程组

(4)计算克里格方程；

(5)进行插值；利用公式对每个点的领域T_i，j(x)都进行插值。

最终获得地球化学栅格化图像数据。

步骤504，将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；具体包括：对地球化学处理地球化学栅格化图像数据和最小噪声变换生成的多通道(多维)图像数据(矩阵)进行匹配处理。在这里主要对比结构的相似性，因此必须参考图像的亮度，对比度和结构三个方面进行度量。假如经过p个波段经过MNF变换后的多维图像数据为同样地球化学栅格化图像数据为X；利用公式计算所述多维图像数据和所述地球化学栅格图像数据的相似性系数；其中， SSIM(X，Y_i)表示地球化学栅格化图像数据X，表示多维图像数据Y中的第i维图像数据Y_i的相似系数；μ_X表示地球化学栅格化图像数据X的均值，μ_Yi表示第i维图像数据Y_i的均值，σ_X表示地球化学栅格化图像数据X的方差，σ_Yi表示第i维图像数据Y_i的方差，σ_XYi是地球化学栅格化图像数据X和第i维图像数据Y_i的协方差,

m表示行数，n表示列数，C₁，C₂，C₃为常数，具体的，通常取>1＝(K1×L)²，C₂＝(K2×L)²，k1可以取值0.01，K2取值0.03，L可以取值 255。从所述多维图像数据中选取与所述地球化学栅格图像数据的相似性系数最大的图像数据作为相似图像数据。

步骤505，利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据；

步骤506，对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像。

步骤506所述对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像，之后还包括：

步骤507，对所述融合后的遥感化学图像进行图像优化分析：

定义y为随机变量(地球化学点数据X)，x₁,x₂,…,x_n为n个自变量(对应的变化后的数据Y_Z)，共观测m次，首先假设y与n个自变量之间存在线性关系：

y＝a₀+a₁×x₁+a₂x₂+…+a_nx_n+ε

式中，a₀,a₁,a₂,…,a_n为回归系数，是常数，表示在其它自变量不变情况下，x_j(j＝1,2,…n)增加或减少一个单位时的平均变化量，ε为去除n个自变量对y>1,x₂,…,x_n之间具有线性关系；(2)每个观测值yj(j＝1,2,…,m)相互独立；>

首先用a₀+a₁×x₁+a₂x₂+…+a_nx_n来估计y的均值E(y)，假定ε服从均值为0，方差为σ²的正态分布，即ε～N(0,σ²),则y服从均值为E(y)，方差为的σ²的正态分布，即y～N[E(y),σ²]，那么m组样本观测数据：

x₁₁,x₁₂,…,x_1n,y₁

x₂₁,x₂₂,…,x_2n,y₂

…………………

x_m1,x_m2,…,x_mn,y_m

式中，x_ij表示x_j在第i次的观测值。有如下公式：

上式为n元线性回归的数学模型，式中，a₀,a₁,a₂,…,a_n为n+1个待定参数，ε1,ε2,…,εm为m个相互独立并且服从同一正态分布的随机变量。为了简化表示，利用矩阵形式：

则n元线性回归的数学模型为

Y＝AX+E

根据公式来进行最小二乘估计，首先假设b0,b1,b2,…,bn分别为n+1个回归系数a₀,a₁,a₂,…,a_n的最小二乘估计值，于是观测值如下表示：

y_j＝b₀x_j1+b₁x_j2+…+b_nx_jn+e_j

e_j为误差ε_j的估计值，称残差，假设为y_j的估计值，那么，

上述式中，j＝1,2,…,m。残差e_j表示实际值y_j与估计值偏离程度。为了使估计值与实际值y_j拟合的最好，必须使残差平方和：

达到最小，根据高等数学原理，极值在函数微分值为0处，确立方程，

从上述公式中得到正规方程，

根据矩阵X，系数两边等式用C和D表示，则，

那么，正规方程的矩阵形式为：CB＝(X′X)B＝X′y＝D。

B为未知向量，如果矩阵系数C满秩，其逆矩阵存在，可以反解出未知向量B，B＝A^-1D＝(X′X)X′y。

步骤508，图像数据检验；

回归方程的假设检验和评价一般采用离差分析法。定义总变异，

SS为回归平方和，是回归值与均值差值的平方和，反映了自变量X 的变化引起Y的波动，自由度df_回＝n(n为自变量个数)。

MS为残差平方和，是实测值y_j与回归值差值的平方和，是由试验误差及其它因素引起的，自由度df_残＝m-n-1。

总变异自由度为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检验

在H0：ai＝0假设下，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，图像输出；

经过以上基本处理后，能够得出的融合后的遥感图像。选择最优层，进行组合式中，S_i为为第i个波段的标准差，R_i，j为第i，j波段的相关系数。OIF越大，表明包含的信息量越大，因此，最大的OIF波段组合为最佳波段组合，采用最佳波段组合合成的图像作为底图。组合后的图像利用RGB合成一幅假彩色图像，这样做出来的图像会反映更多的去除背景噪声后的精细信息，更加适合人的视觉习惯，简单、易用。

输出一幅适合人眼观察的叠加图像。通过软件输出成为JPG或者TIF格式的最终图像。

具体的图像如图7-10所示，其中图7为遥感图像，图8为地球化学图像，图9为地球化学栅格化图像，图10为融合后的图像。

如图11和12所示本发明还提供了一种遥感图像与地球化学的融合系统，所述融合系统包括：

图像获取模块1101，用于获取遥感图像和所述遥感图像对应的地球化学图像。

最小噪声变换模块1102，用于对所述遥感图像进行最小噪声变换，获得多维图像数据；所述最小噪声变换模块1102，具体包括：总协方差矩阵构建子模块，用于根据所述遥感图像，构建所述遥感图像的总协方差矩阵；噪声协方差矩阵构建子模块，用于采用高通滤波器模板对所述遥感图像进行处理，获得所述遥感图像的噪声协方差矩阵；变换矩阵构建子模块，用于根据所述协方差矩阵构建变换矩阵；变换子模块，用于利用所述变换矩阵对所述总协方差矩阵进行变换，获得多维图像数据。

所述噪声协方差矩阵构建子模块，具体包括：

滤波单元，用于利用高通滤波器模块对所述遥感图像进行处理，获得噪声方差估计；噪声协方差矩阵构建单元，用于根据所述噪声方差估计利用公式构建噪声协方差矩阵；其中，C_n为噪声协方差矩阵，是噪声方差估计B_i的倒数。

具体的，最小噪声变换模块1102中设置有最小噪声变换处理器。所述最小噪声变换处理器用于对遥感图像或遥感图像经预处理得到的基础数据。

点面转化模块1103，用于对所述地球化学图像进行点面转化，获得地球化学栅格图像数据；具体的，所述点面转换模块1103中设置有点面转换装置，所述点面转换装置用于将所述地球化学图像转换成栅格数据，即地球化学栅格图像数据。

匹配模块1104，用于将所述多维图像数据和所述地球化学栅格图像数据进行匹配，获得相似图像数据；所述匹配模块1104中设置有多通道匹配处理器，所述多通道匹配处理器用于进行多通道匹配处理。

替换模块1105，用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据，获得替换后的多维图像数据；所述替换模块1105中设置有替换合成装置，所述替换合成装置用于利用所述地球化学栅格图像数据代替所述多维图像数据中的相似图像数据。

最小噪声反变换模块1106，用于对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像。所述最小噪声反变换模块1106中设置有最小噪声反变换处理器，所述最小噪声反变换处理器用于对所述替换后的多维图像数据进行最小噪声反变换，获得融合后的遥感化学图像。

所述融合系统还包括图像优化分析装置、图像输出装置，所述图像优化分析装置用于对所述融合后的遥感化学图像进行参数优化；所述图像输出装置内设置有精度控制器和图像合成装置，所述精度控制器用于对融合化学图像进行检验，提高融合化学图像的精度，所述图像合成装置用于从融合化学图像中选择最优层并进行合成，得到RGB合成的假彩色图像，使图像反映更多的去除背景噪声后的精细信息，更加适合人的视觉习惯、简单、易用，并将最终图像输出。

本发明主要针对遥感图像与地球化学融合难的问题，采用空间变换方式，把水系或者土壤的特征表征出来，采用水系或者土壤地球化学数据点阵变换栅格后替换，从而实现遥感图像与地球化学高度融合，是两者之间的一种重要的补充，能够大幅度提升两种数据的使用性和效率。对于地质找矿、土壤质量调查、农作物调查等具有重要的意义。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。