CN114612387B - 基于特征阈值的遥感影像融合方法、系统、设备及介质 - Google Patents
基于特征阈值的遥感影像融合方法、系统、设备及介质 Download PDFInfo
- Publication number
- CN114612387B CN114612387B CN202210143107.6A CN202210143107A CN114612387B CN 114612387 B CN114612387 B CN 114612387B CN 202210143107 A CN202210143107 A CN 202210143107A CN 114612387 B CN114612387 B CN 114612387B
- Authority
- CN
- China
- Prior art keywords
- band
- characteristic
- function
- vegetation
- enhancement
- 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
- 238000007500 overflow downdraw method Methods 0.000 title claims description 46
- 230000004927 fusion Effects 0.000 claims abstract description 74
- 230000002194 synthesizing effect Effects 0.000 claims abstract description 17
- 238000000034 method Methods 0.000 claims abstract description 16
- 230000006870 function Effects 0.000 claims description 161
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 34
- 238000004364 calculation method Methods 0.000 claims description 23
- 230000014509 gene expression Effects 0.000 claims description 20
- 238000009499 grossing Methods 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 12
- 239000002131 composite material Substances 0.000 claims description 11
- 238000010606 normalization Methods 0.000 claims description 9
- 230000002708 enhancing effect Effects 0.000 claims description 8
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 238000003786 synthesis reaction Methods 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 abstract description 40
- 230000000694 effects Effects 0.000 abstract description 28
- 230000000007 visual effect Effects 0.000 abstract description 15
- 238000010586 diagram Methods 0.000 description 17
- 238000001228 spectrum Methods 0.000 description 14
- 239000011435 rock Substances 0.000 description 9
- 238000004590 computer program Methods 0.000 description 7
- 239000002689 soil Substances 0.000 description 7
- 230000008859 change Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 4
- 239000003086 colorant Substances 0.000 description 4
- 230000007547 defect Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000011156 evaluation Methods 0.000 description 4
- 238000013507 mapping Methods 0.000 description 4
- 230000003287 optical effect Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 238000000513 principal component analysis Methods 0.000 description 3
- 238000011158 quantitative evaluation Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000036541 health Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- 101150029019 ATP6 gene Proteins 0.000 description 1
- 101100026598 Arabidopsis thaliana NOP2A gene Proteins 0.000 description 1
- 238000004737 colorimetric analysis Methods 0.000 description 1
- 150000001875 compounds Chemical group 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- -1 desert Substances 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000003344 environmental pollutant Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 231100000719 pollutant Toxicity 0.000 description 1
- 238000007639 printing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001308 synthesis method Methods 0.000 description 1
- 238000010189 synthetic method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 230000016776 visual perception Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/90—Determination of colour characteristics
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Quality & Reliability (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于特征阈值的遥感影像融合方法、系统、设备及介质，所述方法包括：获取具有原始波段的卫星遥感影像；将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；根据特征比值指数及其特征阈值，构建植被特征增强函数；根据植被特征增强函数，对红波段和绿波段进行植被特征增强；由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。本发明可以提高真彩色遥感影像目视分辨力和计算机解析力，挖掘历史遥感数据真彩色影像模式的应用潜力，改善现状和未来遥感数据真彩色影像模式的应用效果。
Description
技术领域
本发明涉及一种基于特征阈值的遥感影像融合方法、系统、设备及介质，属于卫星遥感影像融合领域。
背景技术
可见光多光谱卫星遥感数据的红、绿、蓝波段分别置于红、绿、蓝通道组合生成的图像为真彩色影像，其显著特征是水体和岩石、土壤、建筑物等裸露陆地物与地面目视的自然色一致，植被色彩就像在万米高空飞机上看到的地面植被色彩。
早期的卫星影像的多波段模式中不少都没有设置蓝色波段，难以合成真彩色影像。随着高分辨率遥感影像的发展，一般都设置包含蓝色波段的多个波段。这一发展引领了遥感卫星波段设置模式，使设置蓝色波段成为一时之风尚。目前绝大多数可见光遥感卫星如资源卫星、海洋水色卫星、气象卫星等均设置了包含红、绿、蓝波段的多个波段。合成真彩色影像的不仅成为可能，而且成为当前遥感影像制图、分类等应用研究的主流模式，如气象遥感卫星和对地观测卫星均把真彩色影像长期作为一项重要的业务化产品，广泛应用于国土、海洋、环境、林业、农业、气象、水利等。
尽管如此，由于真彩色卫星影像数据具有波段间相关性强的特点以及数据中的蓝色波段受水汽等大气干扰的特点，易使地物被笼罩在瑞利散射形成的蓝色“云雾”下，造成真彩色组合图像上植被地物暗淡不清，植被色彩偏蓝灰色调，与地面真实色彩不一致的情况；影像上植被与其它地物的混合像元无法目视识别植被特征，更无法判断植被的类别、健康状况和覆盖度。对于植被广泛分布的区域，造成真彩色影像整体图像质量显著降低，严重影响了真彩色影像的制图应用效果。以地面植被色为参照，卫星遥感真彩色影像与无人机、航空遥感等低空平台获取的真彩色遥感影像相比，存在显著的质量差距。而含有植被波段的标准假彩色图像上植被呈层次多样的红色，目视特征显著，使其成为遥感植被分析技术人员青睐的图种。两相比较，卫星遥感真彩色影像在植被分析与应用方面存在明显的短板。
可见，卫星遥感真彩色影像植被色彩不自然不真实是显而易见的缺陷。保持卫星遥感真彩色影像上的水及岩石、土壤、建筑物等裸露地物的色彩特征，改善卫星遥感真彩色影像的植被色彩特征是卫星遥感真彩色图像处理的关键。直接采用传统直方图调节等现有增强技术可以在一定程度上克服其外在缺点，但难度大，效率低，且不能来克服其固有缺陷。如何应用植被波段的优点来提高真彩色影像的应用潜力成为卫星遥感真彩色图像处理技术研究的难点。针对这一难点，不少技术人员从不同角度对其进行了深入研究，以改善真色彩影像的整体质量。这些研究大致可以分为以下四种思路：
(1)对遥感数据进行瑞利散射纠正，使得色彩信号图像接近于地面真彩色影像(陈春等，遥感信源色彩信号的提取与复现测绘科学，2006年1月，第31卷第1期；韩秀珍等，风云三号D星真彩色影像合成方法研究及应用，海洋气象学报，2019年5月，第39卷第2期)。
(2)利用白平衡法和基于色度学的色彩纠正改善真彩色影像的植被特征，获得更为真实的真彩色图像(游晶等，一种处理彩色多光谱图像的白平衡法，大气与环境光学学报，2012年7月第7卷第4期；黄红莲等，基于人工靶标的多光谱遥感图像真彩色合成，2016年11月，第45卷第11期)。
(3)通过对绿色波段增强处理获得相对较好的真彩色影像。早期主要多采用绿色波段与近红外波段整体加权组合运算方案获得新的绿波段(樊旭艳等，基于主成分分析的遥感图像模拟真彩色融合法，测绘科学技术学报，2006年8月，第23卷第4期；王海燕等，ALOS自然色影像变换及融合方法探讨，测绘技术装备，2012年第14卷第 1期；史园莉等，高分二号卫星遥感影像制图适用性分析，测绘通报，2017年第12期)；后来逐步发展为以归一化植被指数作为分类函数对影像植被像元进行分类加权处理得到新的绿波段(张伟等，基于植被指数的多光谱影像真彩色合成方法，测绘与空间地理信息，2010年12月33卷第6期)。
(4)利用归一化的植被指数分段对绿波段与近红外波段进行Contourlet融合，得到新的绿波段图像(丁慧梅，利用近红外提高多光谱遥感图像颜色自然性研究，硕士论文，2016年)。
发明内容
针对卫星遥感真彩色影像植被等地物暗淡，色彩不自然等缺点，本发明提供了一种基于特征阈值的遥感影像融合方法、系统、设备及介质，其可以提高真彩色遥感影像的目视分辨力和计算机解析力，挖掘历史遥感数据真彩色影像模式的应用潜力，改善现状和未来遥感数据真彩色影像模式的应用效果，为卫星遥感真彩色影像制图在普通用户手机和电脑等网络平台的导航应用、户外和室内等实体空间的观察、分析应用提供基础。
本发明的第一个目的在于提供一种基于特征阈值的遥感影像融合方法。
本发明的第二个目的在于提供一种基于特征阈值的遥感影像融合系统。
本发明的第三个目的在于提供一种计算机设备。
本发明的第四个目的在于提供一种存储介质。
本发明的第一个目的可以通过采取如下技术方案达到：
一种基于特征阈值的遥感影像融合方法，所述方法包括：
获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
根据特征比值指数及其特征阈值，构建植被特征增强函数；
根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
进一步的，所述将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数，统一表达如下式：
其中，n＝1,…,m，m为参与融合的真彩色波段数。
进一步的，所述将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，融合方法为准Brovey融合方法、Gram-Schmidt融合方法和Wavelet融合方法的其中一种；
所述准Brovey融合方法，具体公式如下：
其中，m＝1,2,3；i＝1,…,m；bi、bj表示参与融合的真彩色波段；b′i表示融合结果；
所述Gram-Schmidt或者Wavelet融合方法，具体公式如下：
b′i＝Fusion(NIR,b1,…,bm)
其中，m＝1,2,3，i＝1,…,m，b1,…,bm表示参与融合的真彩色波段；b′i表示融合结果；Fusion表示Gram-Schmidt或Wavelet融合方法。
进一步的，所述根据特征比值指数及其特征阈值，构建植被特征增强函数，具体如下：
根据特征比值指数，确定特征阈值；
利用特征比值指数及其特征阈值，构建阶跃函数；
根据平滑增强要求，构建平滑约束函数；
根据阶跃函数和平滑约束函数，得到植被特征增强函数。
进一步的，所述根据特征比值指数，确定特征阈值，具体如下：
若遥感影像无植被特征，则无需选取特征阈值；
若遥感影像存在植被特征，则特征阈值的取值范围为：
c∈[FIw,FIv]
其中，FIw表示纯水体的特征比值指数的阈值，FIv表示纯植被的特征比值指数的阈值。
进一步的，所述利用特征比值指数及其特征阈值，构建阶跃函数，如下式：
其中，δ(FI)表示阶跃函数，δ(FI)在阈值c处不连续。
进一步的，所述根据平滑增强要求，构建平滑约束函数，具体如下：
设有界连续函数f(FI)；
将阶跃函数δ(FI)与有界连续函数f(FI)相乘，得到函数F(FI)＝δ(FI)f(FI)；
当函数F(FI)满足以下条件时：
即：limFI→cf(FI)＝0满足平滑增强的要求，将所述有界连续函数f(FI)界定为平滑约束函数。
进一步的，所述平滑约束函数的统一表达，如下式：
其中，g(FI)表示平滑约束基函数的复合函数，g(FImax)表示平滑约束基函数的复合函数的最大值，k表示绿度调节系数；
所述平滑约束基函数包括特征比值指数平滑约束基函数和特征比值指数的同源归一化指数平滑约束基函数；
所述特征比值指数平滑约束基函数的表达式，如下式：
u(FI)＝FI-c
所述特征比值指数的同源归一化指数平滑约束基函数的表达式，如下式：
v(FI)＝NDFI-NDFIc
其中，NDFI表示FI的同源归一化指数，NDFIc表示特征阈值c对应的归一化指数。
进一步的，所述根据植被特征增强函数，对红波段和绿波段进行植被特征增强，如下式：
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
所述δ(FI)被替换成δ(NDFI)，如下式：
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
其中，R表示原始红波段，G表示原始绿波段；R′表示植被特征增强后的红波段，G′表示植被特征增强后的绿波段。
本发明的第二个目的可以通过采取如下技术方案达到：
一种基于特征阈值的遥感影像融合系统，所述系统包括：
获取单元，用于获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
计算单元，用于将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
构建单元，用于根据特征比值指数及其特征阈值，构建植被特征增强函数；
增强单元，用于根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
合成单元，用于由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
本发明的第三个目的可以通过采取如下技术方案达到：
一种计算机设备，包括处理器以及用于存储处理器可执行程序的存储器，所述处理器执行存储器存储的程序时，实现上述的遥感影像融合方法。
本发明的第四个目的可以通过采取如下技术方案达到：
一种存储介质，存储有程序，所述程序被处理器执行时，实现上述的遥感影像融合方法。
本发明相对于现有技术具有如下的有益效果：
1、本发明在增强真彩色影像上植被特征的同时，仍保持了真彩色影像上水体及土壤、岩石、建筑物等裸露地物的特征，大幅改善了真彩色影像的整体目视特征和效果，同时真彩色影像各波段间的相关性降低，影像色彩、纹理、层次更加丰富。
2、本发明有效改善了真彩色影像的目视分辨能力和计算机解析能力，有利于挖掘历史遥感数据真彩色影像模式的应用潜力，改善现状和未来遥感数据真彩色影像模式的应用效果，为卫星遥感真彩色影像制图在普通用户手机和电脑等网络平台的导航应用、户外和室内等实体空间的观察、分析应用提供基础。特别是在AI时代，遥感影像目视分类是AI自动分类训练样本选择的关键环节，更加凸显遥感影像目视分类能力在 AI自动分类中的基础作用。
3、本发明依赖于可见光卫星影像的数据内在特征，受主观因素影响小，具有普适性：适用于所有具有近红外、红、绿、蓝光谱数据的多波段卫星影像，也适用于全色波段与多光谱波段通过几何增强融合方法(如主成份分析(PCA)、傅立叶变换(FFT)、小波(Wavelet)变换、Gram-Schimdt变换、Curvelet变换)得到的近红外、红、绿、蓝的多波段卫星影像数据。
4、本发明的实施方案1-6中，增强算法各参数均具有明确的物理意义且表达式简洁明了，有助于提高卫星影像的处理效率：特征比值指数、特征阈值、基于特征阈值的阶跃函数、平滑约束函数和绿度调节系数等均具有明确的物理意义，并且表达式简洁明了，计算快捷高效；融合算法基于叠加原理构建，结构形式简洁明了，可基于图像处理语言快速编制相应的融合增强程序，与真彩色直方图调节等传统图像处理方法相比，处理难度大幅降低，处理效率大幅提高。
5、本发明的实施方案1-6中，增强算法具有多指数、多参数、多模式的特点，可灵活进行参数设置与模式选择，实现不同应用对象个性化需求的量身订制：增强算法其特征比值指数、特征阈值、绿度调节系数等参数均可以根据实际需要灵活选取和设置；增强方案既可以基于特征比值指数平滑约束函数的增强方案，也可以是基于特征比值指数的同源归一化指数平滑约束函数的增强方案，还可以是二者复合约束函数的增强方案，其增强模式及增强成果丰富多样，可满足不同应用人员、不同应用平台、不同用户的个性化需求。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
图1为本发明实施例1的基于特征比值指数的遥感影像融合方法的简易流程图。
图2为本发明实施例1的基于特征比值指数的遥感影像融合方法的具体流程图。
图3为本发明实施例1的增强前的真彩色影像图。
图4a为本发明实施例1的TM典型地物光谱曲线图。
图4b为本发明实施例1的OLI2典型地物光谱曲线图。
图5为本发明实施例1的比值植被指数分类原理示意图。
图6为本发明实施例1的基于特征比值指数的真彩色影像图。
图7为本发明实施例1的基于特征比值指数的归一化指数的真彩色影像图。
图8为本发明实施例1的平滑约束函数族的缺省函数曲线对比图。
图9为本发明实施例1的增强后的红波段图。
图10为本发明实施例1的增强后的绿波段图。
图11为本发明实施例1的增强后的真彩色影像合成图。
图12为本发明实施例2的增强后的红波段图。
图13为本发明实施例2的增强后的绿波段图。
图14为本发明实施例2的增强后的真彩色影像合成图。
图15为本发明实施例3的增强后的红波段图。
图16为本发明实施例3的增强后的绿波段图。
图17为本发明实施例3的增强后的真彩色影像合成图。
图18为本发明实施例4的增强后的红波段图。
图19为本发明实施例4的增强后的绿波段图。
图20为本发明实施例4的增强后的真彩色影像合成图。
图21为本发明实施例5的增强后的红波段图。
图22为本发明实施例5的增强后的绿波段图。
图23为本发明实施例5的增强后的真彩色影像合成图。
图24为本发明实施例6的增强后的红波段图。
图25为本发明实施例6的增强后的绿波段图。
图26为本发明实施例6的增强后的真彩色影像合成图。
图27为本发明实施例7的基于特征阈值的遥感影像融合系统的流程图。
图28为本发明实施例8的计算机设备的结构框图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1：
真彩色影像植被增强效果评价技术包括目视评价和定量评价，具体说明如下：
1)目视评价：正如人们所知，原始的遥感真彩色影像虽然在水体、裸露陆地等地物上与地面地物色彩一致，但在植被区域色彩暗淡、层次不清。一般真彩色影像需要通过对植被的增强获得在上述类地物均与地面色彩一致的真彩色影像。即植被以绿色为基调，不同类别、不同覆盖度的植被呈现深浅不一、浓淡各异的各种绿色；水体以蓝色为基调，除开因水面植被覆盖、高浓度含沙、高污染等组分不同呈现绿、黄、黑等外，主体为深浅不一、浓淡各异的各种蓝色；其它如岩石、裸露土壤、道路、居民地等裸露陆地则与其地面丰富的色彩一致，呈现灰、黑、白、赤、橙、黄、绿、青、蓝、紫等五颜六色。目视选取典型水域、裸露陆地、植被等地物类别，定性比较原始真彩色影像同增强的真彩色影像的色彩变化，可以目视评价真彩色影像植被增强的效果。
2)定量评价：定量评价在一定程度上是目视评价的指标的定量化。对于真彩色影像植被特征增强而言，既要改善真彩色影像植被的色彩，也要保障增强后的真彩色影像的层次、细节等的丰富度和清晰度。总体上可以从两个方面对增强后的影像进行定量评价：其一为植被色彩增强效果的定量描述与比较。在RGB、CMYK、IHS、CIE Lab 等描述色彩空间的模型中，一般认为RGB三原色色彩模型适合于计算机等屏幕显示， CMYK等印刷模型适用于彩色图像打印输出，而IHS、CIE Lab等色彩空间模型色彩描述方面符合人眼视觉感知模式。基于这种认识，一般在定量评价真彩色影像的效果时采用的方法是：将RGB三原色色彩空间描述的遥感影像转换为IHS或CIE Lab色彩空间描述的影像，读取植被地物增强前后在这些色彩空间中的色度、饱和度、强度等，分析其变化趋势和特征。其二，增强后的真彩色图像的质量指标统计与比较。一般来说图像处理的质量可以从三个方面进行评价：第一,增强后图像整体及植被区域的信息丰富程度，可用熵、联合熵度量；第二，增强后图像整体及植被区域的色彩丰富、明亮程度，可用波段统计特征——最大值、最小值、均值、方差及波段之间的相关性指标——相关系数、协方差等衡量；第三，增强后图像整体及植被区域的层次(边缘)、细节(纹理)和影像的清晰程度，可以用梯度、平均梯度等衡量。比较增强前后图像整体及植被区域各指标的差异，就能对光谱(灰度、色调)信息、边缘(层次、差异) 信息、纹理(细节)信息的变化方向进行分析。
本实施例提供了一种基于特征阈值的遥感影像融合方法，该方法基于特征比值指数具有指示水(含泥沙、污染物等各种水体组分的水体)、裸露地物(含土壤、岩石、人工建筑物等)、植被(含各种健康状态、各种类型)的功能，且沿着水-裸露地物- 植被等大类地物的分类轴其值从小到大变化的特征，设计了一个特征比值指数，并根据影像特点确定其一个概略分类的特征阈值；以特征比值指数及其特征阈值构建了一个阶跃函数作为影像概略分类函数；基于影像分类增强的平滑原理，设计了两个平滑约束基函数及其复合函数族；以概略分类函数与平滑约束函数之积建构了一个植被特征增强复合函数，其与真彩色影像红、绿波段的积作为植被特征增强的叠加分量；基于影像融合的叠加增强原理，在真彩色影像的原始红波段、绿波段数据上叠加上述特征增强分量，实现真彩色影像的红、绿两个波段上植被特征的增强；用植被特征增强后的红波段、绿波段与原始的蓝波段合成增强后的真彩色图像。
如图1和图2所示，本实施例提供了一种基于特征阈值的遥感影像融合方法，该方法包括以下步骤：
S101、获取具有原始波段的卫星遥感影像。
打开一幅具有原始波段的卫星遥感影像，其中原始波段包括近红外波段(NIR)、红光波段(R)、绿光波段(G)和蓝光波段(B)；图3为本实施例磨刀门2015年1 月23日植被特征增强前的GF-2真彩色组合彩色影像图(按照envi缺省设置0.1％拉伸的效果图)。
S102、将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数。
本实施例中的特征比值指数分为第一特征比值指数和第二特征比值指数，并以FI作为第一特征比值指数和第二特征比值指数的表达。其中第一特征比值指数是利用准Brovey融合方法得到；第二特征比值指数是利用Gram-Schmidt融合或Wavelet融合方法得到。
(1)利用准Brovey融合方法计算得到第一特征比值指数，具体为：将近红外波段与红、绿、蓝波段中的一个或多个波段进行乘积融合的结果作为分子，以及将融合前的波段或波段组合作为分母，计算得到第一特征比值指数。
第一特征比值指数是指可以反映水、裸露地物、植被等大类地物分布特征的比值指数，可以为比值植被指数(VI)、比值水体指数的倒数(IWI)、比值岩石指数的倒数(IRI)、两波段Brovey指数(DBI)和三波段Brovey指数(TBI)中的一种，具体表达如下：
A、比值植被指数VI为：
对于卫星遥感波段数据而言，VI≥0。
如图4a～图4b所示，遥感影像上水、裸露地物、植被等大类地物的比值植被指数的值从小到大变化，若以该指数作为增强系数由简单的乘积融合方法即能有效增强影像各波段上的植被特征。同时，该指数从水到陆单调递增的性质使之具有地物分类的作用，即总可以找到一个特征阈值把影像分为两类。如找到一个水体阈值就可以将影像分为水与非水两类；找到一个裸露地物值域就可以将影像分为裸露地物与非裸露地物两类；找到一个植被阈值就可以将影像分为植被与非植被两类。当然，由于影像“同谱异物、同物异谱”像元的广泛存在，相较于地物的精确分类而言，这种分类总会有有错分和漏分，只能称之为概略分类。可见，这一指数具有地物分类和植被特征增强的作用。
本实施例设VI的最小值、最大值、水体阈值、纯裸露地物阈值、纯植被阈值、平均值分别为VImin、VImax、VIw、VIb、VIv、VIm，其分类原理示意图如图5所示。
在实际应用中，当以近红外波段为分子，以真彩色的三个波段的组合为分母时，构建的第一特征比值指数均具有比值植被指数相似的功能，因而以下五个第一特征比值指数也是合适的概略分类函数，具体如下：
B、比值水体指数的倒数(IWI)：
C、比值岩石指数的倒数(IRI)：
D、两波段Brovey指数(DBI)：
E、三波段的Brovey指数(TBI)：
(2)利用Gram-Schmidt融合或Wavelet融合方法计算得到第二特征比值指数，具体为：将红外波段与红、绿、蓝波段中的一个或多个波段进行谱间衍生融合的结果作为分子，以及将融合前的波段或波段组合作为分母，计算得到第二特征比值指数；其中，谱间衍生融合采用Gram-Schmidt融合方法或Wavelet融合方法。
本实施例基于上述(1)的融合思路，根据Gram-Schmidt谱间融合方法或Wavelet谱间融合方法计算得到第二特征比值指数，进而扩展概略分类函数的家族。
本实施例中的第二特征比值指数可以为以下五种，具体如下：
F、若采用近红外波段与红波段、近红外波段与绿波段、近红外波段与蓝波段进行谱间融合的方案，则第二特征比值指数的计算公式如下：
其中，SFR、SFG、SFB分别为近红外波段与红波段、近红外波段与绿波段、近红外波段与蓝波段融合后的结果；
G、若采用近红外波段与红波段、绿波段这两个波段进行谱间融合的方案，则第二特征比值指数的计算公式如下：
其中，DFR、DFG为近红外波段与红波段、绿波段两个波段融合后的结果；
H、若采用近红外波段与红波段、蓝波段这两个波段进行谱间融合的方案，则第二特征比值指数的计算公式如下：
其中，DFR、DFB为近红外波段与红波段、蓝波段两个波段融合后的结果；
J、若采用近红外波段与绿波段、蓝波段这两个波段进行谱间融合的方案，则第二特征比值指数的计算公式如下：
其中，DFG、DFB为近红外波段与绿波段、蓝波段两个波段融合后的结果；
K、若采用近红外波段与红波段、绿波段、蓝波段这三个波段进行谱间融合的方案，则第二特征比值指数的计算公式如下：
其中，TFR、TFG、TFB为近红外波段与红波段、绿波段、蓝波段三波段融合后的结果。
本实施例中的特征比值指数VI、IWI、IRI、DBI、RBI、SFI、DFI、TFI都可作为概略分类的特征比值指数；缺省特征指数取Gram-Schmid融合衍生的第二特征比值指数 TF7。
(3)上述(1)和(2)的特征比值指数可按照统一的数学描述进行表达，具体如下：
本实施例设bi为真彩色波段中的一个波段，参与融合的波段为NIR及m个真彩色波段；融合后的波段成果数量为m个，记融合后的成果为b′i。
1、对于准Brovey融合，融合结果记为：
其中，m＝1,2,3，i＝1,m，bi、bj表示参与融合的真彩色波段；
2、对于Gram-Schmidt或者Wavelet融合，融合结果记为：
b′i＝Fusion(NIR,b1,…,bm)
其中，m＝1,2,3，i＝1,…,m，b1,…,bm表示参与融合的真彩色波段；Fusion为 GS(Gram-Schmidt)或者WL(Wavelet)融合方法。
3、特征比值指数(FI)可以统一表达为：
其中，n＝1,…,m，m为参与融合的真彩色波段数。
则特征比值指数(FI)的归一化指数为：
特征比值指数中的缺省值取：
b′i＝GS(NIR,b1,b2,b3)
本实施例以Gram-Schmidt融合衍生的第二特征比值指数为例进行后续运算，具体为：应用ENVI软件中的图像增加工具，以NIR作为增强源，采用Gram-Schmidt融合方法，与真彩色影像的红、绿、蓝三个波段进行融合运算，得到三个新的波段 TFR、TFG、TFB。
根据上述公式，建立波段运算表达式(1.0*b1+b2+b3)/(1.0*b4+b5+b6)计算第二特征比值指数FI，计算成果如图6所示(按照envi缺省设置0.1％拉伸的效果图)。其中， b1、b2、b3分别表示融合后的红、绿、蓝波段，b4、b5、b6分别表示原始红、绿、蓝波段。
进一步地，建立波段运算表达式(1.0*b1-1)/(1.0*b1+1)计算得到NDFI，计算成果如图7所示(按照envi缺省设置0.1％拉伸的效果图)。其中，b1表示第二特征比值指数 FI。
S103、根据特征比值指数及其特征阈值，构建植被特征增强函数。
在本实施例中，步骤S103具体包括：
S1031、根据特征比值指数，确定特征阈值。
从分类增强的角度出发，特征阈值可以是特征比值指数最大值与最小值中的任何值。但正如前述，待增强目标既要保持真彩色影像上水体及岩石、土壤、建筑物等裸露地物的特征，又要增强植被特征；对于一幅具有植被的影像而言，植被特征除了纯植被像元的植被特征之外，还包含大量的含有植被的混合像元的植被特征。理论上讲，特征阈值取得越小，包含植被特征的像元数越多；反之，越少。
设纯植被的特征比值指数的阈值为FIv，特征比值指数的平均值为FIm，纯裸露地物的阈值为FIb，纯水体的特征比值指数的阈值为FIw。按大于阈值分类，取阈值为FIm较之阈值FIv将包含更多的含有植被的混合像元；取阈值FIb将包含所有含植被的混合像元；如果影像上存在水体，对于纯水体而言，不包含任何植被特征，取阈值为FIw，既包含了纯裸露地物、纯植被，也包含所有的含植被的混合像元。可见无论阈值取纯水体阈值、取平均值、纯裸露地物阈值还是取纯植被阈值，都能起到植被特征增强的作用。
综上可知，取FIw或者FIb作为特征阈值，可以对影像上所有像元的植被特征进行增强，而又保持了纯水体及陆地上纯裸露地物的特征；取FIv作为特征阈值，只对纯植被特征进行增强，而遗漏了含植被的混合像元；取平均值FIm作为特征阈值，植被增强的像元数介于上述二者之间。故一般而言，特征阈值的取值范围为：
c∈[FIw,FIv]
鉴于影像地物组成的复杂性，无水体的影像无所谓FIw，有水体的影像，FIw、FIb、FIv等特征值需要通过人机交互确定，相对困难；而FIm可以直接通过统计计算得到，相对容易。一般取FIm为缺省值，也可以根据影像实际特点或用户增强需要可以取其它值，如FIw、FIb或者FIv，或者其它值。
如果影像上无植被，而仅是水体、沙漠、土壤、岩石、建筑物等裸露地物，这样的真彩色影像事实上无需增强处理，也就不必取特征阈值。
基于上述论述，本实施例应用ENVI统计工具计算FI的平均值FIm，并将平均值FIm作为本实施例的特征阈值c。其中，c＝FIm＝1.027063，NDFIc＝0.013351。
S1032、基于特征比值指数及其特征阈值，得到平滑约束函数族。
(1)基于特征阈值的阶跃函数的构建。
若以FI统一表达上述特征比值指数家族，设c为特征阈值，则对于一幅具体的卫星遥感影像而言，FI存在相应的统计最大值和统计最小值，分别记为：FImin，FImax， c∈(FImin,FImax)，则二分影像上地物的概略分类函数可以用一个阶跃函数表达如下：
显然，阶跃函数δ(FI)在阈值c处不连续。
设NDFI为FI的归一化指数，阈值c处的归一化指数为NDFIc，根据二者之间的等价性，上述阶跃函数可以用下列等价形式表述：
(2)基于特征阈值的平滑约束函数的构建
设有f(FI)为有界连续函数，δ(FI)为上述阶跃函数，当用这两个函数之积 F(FI)＝δ(FI)f(FI)对影像波段b的植被特征进行增强时，按照叠加融合的原理可构造如下融合表达式：
b′＝[δ(FI)f(FI)+1]×b
由于阶跃函数在阈值c处为不连续函数，容易使得融合后的影像出现锯齿、斑点等噪声，使得影像上的地物之间色彩过渡不平顺、不平滑。这种现象从数学层面上分析是因为增强后的影像破坏了原始影像的连续性，要使增强后的影像b′依然保持平滑的充分必要条件是函数F(FI)＝δ(FI)f(FI)为连续函数。即：
上述条件等价于：
满足上述约束条件的f(FI)称之为平滑约束函数。
显而易见，这样的函数无穷多，其中两个简洁的平滑约束基函数为：
u(FI)＝FI-c
或者
v(FI)＝NDFI-NDFIc
其中，u(FI)表示特征比值指数平滑约束基函数，v(FI)表示特征比值指数的同源归一化指数平滑约束基函数；NDFI为FI的同源归一化指数，NDFIc为阈值c对应的归一化指数，表达式分别为：
进一步地，由上述两个平滑约束基函数可以构建系列平滑约束函数，形成平滑约束函数族，其中六个典型的平滑约束函数如下：
(1)特征比值指数平滑约束函数：
其中，k为大于0的常系数，缺省值k＝0.5。
(2)归一化特征比值指数平滑约束函数：
其中，k为大于0的常系数，缺省值取k＝0.5。
(3)特征比值指数及其归一化特征比值指数乘积平滑约束函数：
其中，c1∈[0,c]，缺省值取c1＝c，k为大于0的常系数，缺省值取k＝0.5。
或者写为其等价表达式：
其中，c1∈[0,c]，缺省值取c1＝c，k为大于0的常系数，缺省值取k＝0.5。
(4)特征比值指数自乘平滑约束函数：
其中，c1∈[0,c]，缺省值取c1＝c，k为大于0的常系数，缺省值取k＝0.5。
(5)归一化特征比值指数自乘平滑约束函数：
其中，c1∈[0,c]，缺省值取c1＝c，k为大于0的常系数，缺省值取k＝0.5。
(6)特征比值指数及其归一化特征比值指数线性组合平滑约束函数：
其中，k为大于0的常系数，ck∈[0,1/(FImax-c)]，缺省值取ck＝0，缺省值取 k＝0.5。
上述六个平滑约束函数可统一表达为：
其中，g(FI)表示平滑约束基函数的复合函数，g(FImax)表示平滑约束复合函数的最大值，k表示常系数，具有调节植被绿度的作用，称之为绿度调节系数。
进一步地，当k从小到大变换时，植被绿度从原始真彩色影像的蓝灰色到青绿色到黄绿色变化，缺省值取k＝0.5。依据实际情况，可取更大或更小的k值，可供参考的几个特征值有：(NDFImax-NDFIc)、(NDFImax-NDFIc)2、(FImax-c)(NDFImax- NDFIc)、(FImax-c)、2(FImax-c)-(c+1)2(NDFImax-NDFIc)和(FImax-c)2。
本实施中的平滑约束函数族的缺省函数曲线如图8所示。
本实施例应用ENVI统计工具得到FI的最大值FImax，FImax＝2.279753， NDFImax＝0.390198。
S104、根据植被特征增强函数，对红波段和绿波段进行植被特征增强。
基于步骤S103所得的平滑约束函数族，本实施例选取特征比值指数平滑约束函数来构建植被特征增强波段，具体如下：
S1041、设真彩色影像的红、绿、蓝波段分别为R、G、B，按照叠加融合原理对红、绿波段进行植被特征增强，得到植被特征增强波段，具体公式如下：
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
或者：
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
其中，δ(FI)也可以用等价的δ(NDFI)替换。
S1042、利用特征比值指数平滑约束函数计算植被特征增强波段，具体如下：
g1(FI)＝FI-c＝FI-1.027063
g1(FImax)＝FImax-c＝1.25269
k＝0.5
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
根据上述公式，在ENVI中输入计算表达式：uint(((b1 gt 1.027063)*(b1-1.027063) *0.5/1.25269+1)*b2)；其中b1为特征比值指数FI，b2为真彩色影像的红波段或绿波段。植被特征增强的计算结果如图9和10所示(按照envi缺省设置0.1％拉伸的效果图)。
S105、由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
将增强后的红波段、绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影像，如图11所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见表1和表2；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表3。
表1植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表2植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表3植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
实施例2：
本实施例中的步骤与实施例1的步骤基本一致，但是在具体实现上略有不同，不同点在于步骤S4022。
进一步地，本实施例利用归一化特征比值指数平滑约束函数计算植被特征增强波段，具体如下：
g2(FI)＝NDFI-NDFIc＝NDFI-0.013351
g2(FImax)＝NDFImax-NDFIc＝0.376847
k＝0.5
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
根据上述公式输入计算表达式：uint(((b1gt0.013351)*(b1-0.013351)* 0.5/0.376847+1)*b2)；其中b1为特征比值指数的归一化指数NDFI，b2为真彩色影像的红波段或绿波段。植被特征增强的计算结果如下图12和13(按照envi缺省设置0. 1％拉伸的效果图)。
将上述增强后的红波段、增强后的绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影像，如图14所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见表4～表5；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表6。
表4植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表5植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表6植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
实施例3：
本实施例中的步骤与实施例1的步骤基本一致，但是在具体实现上略有不同，不同点在于步骤S4022。
进一步地，本实施例利用特征比值指数及其归一化特征比值指数乘积平滑约束函数计算植被特征增强波段，具体如下：
g3(FI)＝(FI-c)(NDFI-NDFIc)＝(FI-1.027063)(NDFI-0.013351)
g3(FImax)＝(FImax-c)(NDFImax-NDFIc)＝0.472073
k＝0.5
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
根据上述公式，在ENVI中输入计算表达式：uint(((b1 gt 1.027063)*(b1-1.027063) *(b2-0.013351)*0.5/0.472073+1)*b3)；其中b1为特征比值指数，b2为特征比值指数的归一化指数NDFI，b3为真彩色影像的红波段或绿波段；植被特征增强的计算结果如图 15和16所示(按照envi缺省设置0.1％拉伸的效果图)。
将上述增强后的红波段、增强后的绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影像，如图17所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见表7～表8；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表9。
表7植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表8植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表9植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
实施例4：
本实施例中的步骤与实施例1的步骤基本一致，但是在具体实现上略有不同，不同点在于步骤S4022。
进一步地，本实施例利用特征比值指数自乘平滑约束函数计算植被特征增强波段，具体如下：
g4(FI)＝(FI-c)2＝(FI-1.027063)2
g4(FImax)＝(FImax-c)2＝1.569232
k＝0.5
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
根据上述公式，在ENVI中输入计算表达式：uint(((b1 gt 1.027063)*(b1-1.027063) ^2*0.5/1.569232+1)*b2)；其中b1为特征比值指数FI，b2为真彩色影像的红波段或绿波段；植被特征增强的计算结果如图18和19所示(按照envi缺省设置0.1％拉伸的效果图)。
将上述增强后的红波段、增强后的绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影像，如图20所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见表10～表11；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表12。
表10植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表11植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表12植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
实施例5：
本实施例中的步骤与实施例1的步骤基本一致，但是在具体实现上略有不同，不同点在于步骤S4022。
进一步地，本实施例利用归一化特征比值指数自乘平滑约束函数计算植被特征增强波段，具体如下：
g5(FI)＝(NDFI-NDFIc)2＝(NDFI-0.013351)2
g5(FImax)＝(NDFImax-NDFIc)2＝0.142014
k＝0.5
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
根据上述公式输入计算表达式：uint(((b1 gt 0.013351)*(b1-0.013351)^2*0.5/0.1420 14+1)*b2)；其中b1为特征比值指数的归一化指数NDFI，b2为真彩色影像的红波段或绿波段；植被特征增强的计算结果如图21和22所示(按照envi缺省设置0.1％拉伸的效果图)。
将上述增强后的红波段、增强后的绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影像，如图23所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见见表13～表14；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表15。
表13植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表14植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表15植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
实施例6：
本实施例中的步骤与实施例1的步骤基本一致，但是在具体实现上略有不同，不同点在于步骤S4022。
进一步地，本实施例利用特征比值指数及其归一化特征比值指数线性组合平滑约束函数计算植被特征增强波段，具体如下：
当
ck＝0.4∈[0,0.798282]时
g6(FI)＝1.442769(FI-1.027063)-2.14236(NDFI-0.376847)
g6(FImax)＝1
k＝0.5
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
根据上述公式输入计算表达式：uint(((b1 gt 1.027063)*(1.442769*(b1-1.027063)-2. 14236*(b2-0.013351))*0.5+1)*b3)；其中b1为特征比值指数，b2为特征比值指数的归一化指数NDFI，b3为真彩色影像的红波段或绿波段；植被特征增强的计算结果如图24 和25所示(按照envi缺省设置0.1％拉伸的效果图)。
将上述增强后的红波段、增强后的绿波段与原始蓝波段按照红、绿、蓝通道合成彩色影，如图26所示(按照envi缺省设置0.1％拉伸的效果图)；存储增强后的图像文件即为所得。
本实施例植被特征增强后的彩色影像成果与原始真彩色影像RGB色彩模式下的全图及增强区域的统计特征对比分析表见表16～表17；HLS色彩模式下的全图及增强区域的统计特征对比分析表见表18。
表16植被特征增强后的真彩色影像与原始真彩色影像RGB模式统计特征对比分析表
表17植被特征增强后的真彩色影像与原始真彩色影像RGB模式植被区域统计特征对比分析表
表18植被特征增强后的真彩色影像与原始真彩色影像HLS模式统计特征对比分析表
根据上述实施例1-6中的表1-18可知，原始真彩色影像的植被信息存在不突出、不自然和不真实情况，如：纯植被像元色调偏灰，容易与其它地物混淆；混合植被像元植被色彩淹没在其它地物色彩中，植被信息被掩盖。这两种情况会造成真彩色影像植被信息的目视分辨率低，直接影响了真彩色影像模式应用效果；而经本发明增强后的真彩色影像，其各植被像元点都被增强，其植被色彩、纹理和层次也得到全面改善，有效提高了真彩色影像植被信息的目视分辨力和计算机解析能力，提高了真彩色影像模式植被分析能力和效果。
本领域技术人员可以理解，实现上述实施例方法中的全部或部分步骤可以通过程序来指令相关的硬件来完成，相应的程序可以存储于计算机可读取存储介质中。
应当注意，尽管在附图中以特定顺序描述了上述实施例的方法操作，但是这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。相反，描绘的步骤可以改变执行顺序。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
实施例7：
如图27所示，本实施例提供了一种基于特征比值指数的遥感影像融合系统，该系统包括获取单元2701、计算单元2702、构建单元2703、增强单元2704和合成单元2705，各个单元的具体功能如下：
获取单元2701，用于获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
计算单元2702，用于将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
构建单元2703，用于根据特征比值指数及其特征阈值，构建植被特征增强函数；
增强单元2704，用于根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
合成单元2705，用于由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
本实施例中各个单元的具体实现可以参见上述实施例1，在此不再一一赘述；需要说明的是，本实施例提供的系统仅以上述各功能单元的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配给不同的功能单元完成，即将内部结构划分成不同的功能单元，以完成以上描述的全部或者部分功能。
实施例8：
如图28所示，本实施例提供了一种计算机设备，其包括通过系统总线2801连接的处理器2802、存储器、输入装置2803、显示装置2804和网络接口2805。其中，处理器2802用于提供计算和控制能力，存储器包括非易失性存储介质2806和内存储器 2807，该非易失性存储介质2806存储有操作系统、计算机程序和数据库，该内存储器 2807为非易失性存储介质2806中的操作系统和计算机程序的运行提供环境，计算机程序被处理器2802执行时，实现上述实施例1的遥感影像融合方法，如下：
获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
根据特征比值指数及其特征阈值，构建植被特征增强函数；
根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
实施例9：
本实施例提供一种存储介质，该存储介质为计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现上述实施例1的遥感影像融合方法，如下：
获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
根据特征比值指数及其特征阈值，构建植被特征增强函数；
根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像。
需要说明的是，本实施例的计算机可读存储介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是但不限于电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器 (CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。
在本实施例中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本实施例中，计算机可读信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读信号介质还可以是计算机可读存储介质以外的任何计算机可读存储介质，该计算机可读信号介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读存储介质上包含的计算机程序可以用任何适当的介质传输，包括但不限于：电线、光缆、RF(射频)等等，或者上述的任意合适的组合。
上述计算机可读存储介质可以以一种或多种程序设计语言或其组合来编写用于执行本实施例的计算机程序，上述程序设计语言包括面向对象的程序设计语言—诸如Java、Python、C++，还包括常规的过程式程序设计语言—诸如C语言或类似的程序设计语言。程序可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络，包括局域网(LAN)或广域网(WAN)连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
综上所述，本发明的理论基础扎实、物理意义明确，应用对象广泛。增强后的影像色彩鲜明、信息丰富、易于目视和自动分类，尤其在当前高分辨率卫星遥感快速发展背景下，对于推动高分辨率影像在国内外各行各业的推广应用具有巨大的促进作用。
以上所述，仅为本发明专利较佳的实施例，但本发明专利的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明专利所公开的范围内，根据本发明专利的技术方案及其发明构思加以等同替换或改变，都属于本发明专利的保护范围。
Claims (6)
1.一种基于特征阈值的遥感影像融合方法，其特征在于，所述方法包括：
获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
根据特征比值指数及其特征阈值，构建植被特征增强函数；
根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像；
所述将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数，统一表达如下式：
其中，n＝1,…,m，m为参与融合的真彩色波段数，bi为真彩色波段中的一个波段，参与融合的波段为NIR及m个真彩色波段，b′i为融合后的成果；
所述根据特征比值指数及其特征阈值，构建植被特征增强函数，具体如下：
根据特征比值指数，确定特征阈值；
利用特征比值指数及其特征阈值，构建阶跃函数；
根据平滑增强要求，构建平滑约束函数；
根据阶跃函数和平滑约束函数，得到植被特征增强函数；
所述根据特征比值指数，确定特征阈值，具体如下：
若遥感影像无植被特征，则无需选取特征阈值；
若遥感影像存在植被特征，则特征阈值的取值范围为：
c∈[FIw,FIv]
其中，FI表示特征比值指数，FIw表示纯水体的特征比值指数的阈值，FIv表示纯植被的特征比值指数的阈值；
所述利用特征比值指数及其特征阈值，构建阶跃函数，如下式：
其中，δ(FI)表示阶跃函数，δ(FI)在阈值c处不连续；
或用等价形式表示，如下式：
其中，NDFI为FI的同源归一化指数，NDFIc为阈值c对应的归一化指数值；
所述根据平滑增强要求，构建平滑约束函数，具体如下：
设有界连续函数f(FI)；
将阶跃函数δ(FI)与有界连续函数f(FI)相乘，得到函数F(FI)＝δ(FI)f(FI)；
当函数F(FI)满足以下条件时：
即：limFI→cf(FI)＝0满足平滑增强的要求，将所述有界连续函数f(FI)界定为平滑约束函数；
所述根据植被特征增强函数，对红波段和绿波段进行植被特征增强，如下式：
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
或如下式：
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
其中，R表示原始红波段，G表示原始绿波段；R′表示植被特征增强后的红波段，G′表示植被特征增强后的绿波段。
2.根据权利要求1所述的遥感影像融合方法，其特征在于，所述将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，融合方法为准Brovey融合方法、Gram-Schmidt融合方法和Wavelet融合方法的其中一种；
所述准Brovey融合方法，具体公式如下：
其中，m＝1,2,3；i＝1,…,m；NIR表示近红外波段；bi、bj表示参与融合的真彩色波段；b′i表示融合结果；
所述Gram-Schmidt或者Wavelet融合方法，具体公式如下：
b′i＝Fusion(NIR,b1,…,bm)
其中，m＝1,2,3，i＝1,…,m，b1,…,bm表示参与融合的真彩色波段；b′i表示融合结果；Fusion表示Gram-Schmidt或Wavelet融合方法。
3.根据权利要求1所述的遥感影像融合方法，其特征在于，所述平滑约束函数的统一表达，如下式：
其中，g(FI)表示平滑约束基函数的复合函数，g(FImax)表示平滑约束基函数的复合函数的最大值，k表示绿度调节系数；
所述平滑约束基函数包括特征比值指数平滑约束基函数和特征比值指数的同源归一化指数平滑约束基函数；
所述特征比值指数平滑约束基函数的表达式，如下式：
u(FI)＝FI-c
所述特征比值指数的同源归一化指数平滑约束基函数的表达式，如下式：
v(FI)＝NDFI-NDFIc
其中，NDFI表示FI的同源归一化指数，NDFIc表示特征阈值c对应的归一化指数。
4.一种基于特征阈值的遥感影像融合系统，其特征在于，所述系统包括：
获取单元，用于获取具有原始波段的卫星遥感影像，所述原始波段包括近红外波段、红波段、绿波段和蓝波段；
计算单元，用于将近红外波段与红波段、绿波段和蓝波段中的一个或多个波段进行融合，将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数；
构建单元，用于根据特征比值指数及其特征阈值，构建植被特征增强函数；
增强单元，用于根据植被特征增强函数，对红波段和绿波段进行植被特征增强；
合成单元，用于由增强后的红波段、增强后的绿波段以及原始蓝波段，合成得到最终的彩色影像；
所述将融合后的波段或波段组合作为分子，以及将融合前的波段或波段组合作为分母，计算得到特征比值指数，统一表达如下式：
其中，n＝1,…,m，m为参与融合的真彩色波段数，bi为真彩色波段中的一个波段，参与融合的波段为NIR及m个真彩色波段，b′i为融合后的成果；
所述根据特征比值指数及其特征阈值，构建植被特征增强函数，具体如下：
根据特征比值指数，确定特征阈值；
利用特征比值指数及其特征阈值，构建阶跃函数；
根据平滑增强要求，构建平滑约束函数；
根据阶跃函数和平滑约束函数，得到植被特征增强函数；
所述根据特征比值指数，确定特征阈值，具体如下：
若遥感影像无植被特征，则无需选取特征阈值；
若遥感影像存在植被特征，则特征阈值的取值范围为：
c∈[FIw,FIv]
其中，FI表示特征比值指数，FIw表示纯水体的特征比值指数的阈值，FIv表示纯植被的特征比值指数的阈值；
所述利用特征比值指数及其特征阈值，构建阶跃函数，如下式：
其中，δ(FI)表示阶跃函数，δ(FI)在阈值c处不连续；
或用等价形式表示，如下式：
其中，NDFI表示FI的同源归一化指数，NDFIc为阈值c对应的归一化指数值；
所述根据平滑增强要求，构建平滑约束函数，具体如下：
设有界连续函数f(FI)；
将阶跃函数δ(FI)与有界连续函数f(FI)相乘，得到函数F(FI)＝δ(FI)f(FI)；
当函数F(FI)满足以下条件时：
即：limFI→cf(FI)＝0满足平滑增强的要求，将所述有界连续函数f(FI)界定为平滑约束函数；
所述根据植被特征增强函数，对红波段和绿波段进行植被特征增强，如下式：
R′＝[δ(FI)f(FI)+1]×R
G′＝[δ(FI)f(FI)+1]×G
或如下式：
R′＝[δ(NDFI)f(FI)+1]×R
G′＝[δ(NDFI)f(FI)+1]×G
其中，R表示原始红波段，G表示原始绿波段；R′表示植被特征增强后的红波段，G′表示植被特征增强后的绿波段。
5.一种计算机设备，包括处理器以及用于存储处理器可执行程序的存储器，其特征在于，所述处理器执行存储器存储的程序时，实现权利要求1-3任一项所述的遥感影像融合方法。
6.一种存储介质，存储有程序，其特征在于，所述程序被处理器执行时，实现权利要求1-3任一项所述的遥感影像融合方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210143107.6A CN114612387B (zh) | 2022-02-16 | 2022-02-16 | 基于特征阈值的遥感影像融合方法、系统、设备及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210143107.6A CN114612387B (zh) | 2022-02-16 | 2022-02-16 | 基于特征阈值的遥感影像融合方法、系统、设备及介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114612387A CN114612387A (zh) | 2022-06-10 |
| CN114612387B true CN114612387B (zh) | 2023-02-10 |
Family
ID=81859984
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210143107.6A Active CN114612387B (zh) | 2022-02-16 | 2022-02-16 | 基于特征阈值的遥感影像融合方法、系统、设备及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114612387B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116935222B (zh) * | 2023-07-25 | 2024-02-13 | 珠江水利委员会珠江水利科学研究院 | 一种遥感影像蓝波段模拟方法、系统及可读存储介质 |
| CN117237770A (zh) * | 2023-08-29 | 2023-12-15 | 珠江水利委员会珠江水利科学研究院 | 基于指数幂函数的遥感影像融合方法、系统、设备及介质 |
| CN117437164A (zh) * | 2023-11-08 | 2024-01-23 | 武汉大势智慧科技有限公司 | 一种三维模型纹理增强方法、装置、电子设备和介质 |
| CN117333402B (zh) * | 2023-12-01 | 2024-03-01 | 北京飞渡科技股份有限公司 | 倾斜摄影图像的局部更新处理方法及电子设备 |
| CN117830119A (zh) * | 2023-12-27 | 2024-04-05 | 珠江水利委员会珠江水利科学研究院 | 基于指数缩放变换的遥感影像融合方法、系统、设备及介质 |
| CN117475272B (zh) * | 2023-12-27 | 2024-04-05 | 珠江水利委员会珠江水利科学研究院 | 基于指数函数空间的特征曲线的遥感影像融合方法及系统 |
| CN117935079A (zh) * | 2024-01-29 | 2024-04-26 | 珠江水利委员会珠江水利科学研究院 | 一种遥感影像融合方法、系统及可读存储介质 |
| CN118015420A (zh) * | 2024-03-28 | 2024-05-10 | 珠江水利委员会珠江水利科学研究院 | 基于植被盖度特征幂次函数的遥感影像融合方法、系统及介质 |
| CN117975224A (zh) * | 2024-03-29 | 2024-05-03 | 珠江水利委员会珠江水利科学研究院 | 基于指数特征幂次函数的遥感影像融合方法、系统及介质 |
| CN117994668B (zh) * | 2024-04-07 | 2024-06-11 | 国家海洋环境监测中心 | 一种翅碱蓬覆盖面积的遥感图像计算方法、系统、设备及介质 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10586105B2 (en) * | 2016-12-30 | 2020-03-10 | International Business Machines Corporation | Method and system for crop type identification using satellite observation and weather data |
| US10997707B1 (en) * | 2018-02-27 | 2021-05-04 | Orbital Sidekick, Inc. | Aerial and space-based hyperspectral imaging system and method for hydrocarbon and chemical industry regulatory compliance, leak detection and product speciation |
| CN109472244B (zh) * | 2018-11-15 | 2021-10-01 | 珠江水利委员会珠江水利科学研究院 | 一种基于分离指数的土壤、岩石和植被识别方法 |
| CN112287871B (zh) * | 2020-11-12 | 2023-01-17 | 广东海洋大学 | 基于多特征与光谱融合的近岸水产养殖区遥感图像提取方法 |
| CN113240758B (zh) * | 2021-05-28 | 2022-03-08 | 珠江水利委员会珠江水利科学研究院 | 基于融合衍生指数的遥感影像融合方法、系统、设备及介质 |
-
2022
- 2022-02-16 CN CN202210143107.6A patent/CN114612387B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN114612387A (zh) | 2022-06-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114612387B (zh) | 基于特征阈值的遥感影像融合方法、系统、设备及介质 | |
| CN113240758B (zh) | 基于融合衍生指数的遥感影像融合方法、系统、设备及介质 | |
| CN113222054B (zh) | 基于特征比值指数的遥感影像融合方法、系统、设备及介质 | |
| Herold et al. | Object-oriented mapping and analysis of urban land use/cover using IKONOS data | |
| CN117237770A (zh) | 基于指数幂函数的遥感影像融合方法、系统、设备及介质 | |
| CN101840581A (zh) | 一种从卫星遥感影像中提取建筑物轮廓的方法 | |
| CN112733596A (zh) | 一种基于中高空间分辨率遥感影像融合的森林资源变化监测方法及应用 | |
| Hovi et al. | Quantifying the missing link between forest albedo and productivity in the boreal zone | |
| CN117475272B (zh) | 基于指数函数空间的特征曲线的遥感影像融合方法及系统 | |
| CN110779876A (zh) | 一种疫木识别方法、装置、设备及计算机可读存储介质 | |
| CN112052757B (zh) | 火烧迹地信息提取方法、装置、设备和存储介质 | |
| CN110879992A (zh) | 基于迁移学习的草原地表覆盖物分类方法和系统 | |
| CN114821349A (zh) | 顾及谐波模型系数和物候参数的森林生物量估算方法 | |
| Tompalski et al. | LiDAR and VHRS data for assessing living quality in cities–An approach based on 3D spatial indices | |
| Smith et al. | A new approach to quantifying abundances of materials in multispectral images | |
| Cilek et al. | The use of regression tree method for Sentinel-2 satellite data to mapping percent tree cover in different forest types | |
| CN107239886B (zh) | 基于高分卫星遥感数据的gdp密度分析系统 | |
| Zabala et al. | Segmentation and thematic classification of color orthophotos over non-compressed and JPEG 2000 compressed images | |
| Danoedoro et al. | Integrating spectral indices and geostatistics based on Landsat-8 imagery for surface clay content mapping in Gunung Kidul area, Yogyakarta, Indonesia | |
| CN116994071A (zh) | 一种基于自适应光谱残差的多光谱激光雷达点云分类方法 | |
| Matabishi et al. | Multiple endmember spectral mixture analysis of desis image to identify rooftops in Kigali | |
| CN111257854B (zh) | 一种通用的基于遥感影像分割单元的地形校正优化方法 | |
| CN112686861A (zh) | 一种土地利用变化遥感监测分析方法、装置及智能终端 | |
| CN112800993A (zh) | 一种基于离散小波变换的区域尺度工业厂房信息增强方法 | |
| Colaninno et al. | Defining Densities for Urban Residential Texture, through Land use Classification, from Landsat TM Imagery: Case Study of Spanish Mediterranean Coast |
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 |