CN112233190A - 一种基于区域网平差的卫星遥感影像色彩均衡方法 - Google Patents

一种基于区域网平差的卫星遥感影像色彩均衡方法 Download PDF

Info

Publication number
CN112233190A
CN112233190A CN202010425601.2A CN202010425601A CN112233190A CN 112233190 A CN112233190 A CN 112233190A CN 202010425601 A CN202010425601 A CN 202010425601A CN 112233190 A CN112233190 A CN 112233190A
Authority
CN
China
Prior art keywords
image
grids
area
equation
color
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.)
Granted
Application number
CN202010425601.2A
Other languages
English (en)
Other versions
CN112233190B (zh
Inventor
刘世杰
童小华
许雄
卢其剑
高晓峰
金雁敏
谢欢
魏超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tongji University
Original Assignee
Tongji University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tongji University filed Critical Tongji University
Priority to CN202010425601.2A priority Critical patent/CN112233190B/zh
Publication of CN112233190A publication Critical patent/CN112233190A/zh
Application granted granted Critical
Publication of CN112233190B publication Critical patent/CN112233190B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/90Determination of colour characteristics
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/22Matching criteria, e.g. proximity measures
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/40Image enhancement or restoration using histogram techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/90Dynamic range modification of images or parts thereof
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/90Dynamic range modification of images or parts thereof
    • G06T5/94Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/20Analysis of motion
    • G06T7/269Analysis of motion using gradient-based methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Multimedia (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种基于区域网平差的卫星遥感影像色彩均衡方法,包括:S1、根据影像空间位置提取影像重叠区,对影像重叠区划分均匀格网,利用HOG特征提取算法提取格网特征,利用格网特征进行同名格网的优化筛选;S2、以多项式为平差模型,针对系数阵秩亏问题引入基准影像,针对像素值分布不均问题引入虚拟控制,针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束,对测区内各影像的色彩变换参数进行引入权值的整体求解;S3、利用S2得到的色彩变换参数对每景影像做色彩校正,得到色彩均衡的影像序列。与现有技术相比,本发明方法能够有效地处理相邻影像间的亮度和色彩差异,使整个测区内的颜色过渡自然,取得了较好的效果。

Description

一种基于区域网平差的卫星遥感影像色彩均衡方法
技术领域
本发明涉及卫星遥感影像处理技术,尤其是涉及一种基于区域网平差的卫星遥感影像色彩均衡方法。
背景技术
遥感影像是获取全球地理信息的重要数据来源。然而,由于光学遥感影像的成像时间,成像时的光照情况,以及各景影像的地物地形不同等各种因素的影响,会造成影像间的色彩不一致现象。这种现象是影响后续数字产品生产质量的直接因素,即使影像采取了高精度的几何定位和几何拼接方法,但由于影像间的色彩不一致,仍会出现明显的影像分界线,如若不加以处理,在视觉上很难做到测区“一张图”的效果。
目前对于卫星遥感影像的匀色处理研究较少,不同于近景影像和航空影像,卫星影像的重叠度远小于前者,如若将重叠区域的整体灰度均值作为平差变量,那么在求解时其所求参数远多于方程数,而导致无法求解结果;如若采用诸如sift等匹配方法提取同名点用于平差计算,在一景遥感影像中其像素数量可达千万甚至数亿数量级,在重叠区也会有数百万的像素数量,即使在其中只提取20%的点用于平差计算,那么整个测区也会有数百万个同名点,在后期平差计算时会造成系数阵过大,导致计算开销过大甚至出现无法解算的情况,且重叠区百分之二十的同名点可能在部分影像中无法体现两景影像整体的灰度差异,导致最后影像色彩均衡处理的失败。
发明内容
本发明的目的就是为了克服上述现有技术存在的卫星遥感影像间的色彩不一致的缺陷而提供一种基于区域网平差的卫星遥感影像色彩均衡方法。
本发明的目的可以通过以下技术方案来实现:
一种基于区域网平差的卫星遥感影像色彩均衡方法,包括:
S1、根据影像空间位置提取影像重叠区,对影像重叠区划分均匀格网,利用HOG特征提取算法提取格网特征,利用格网特征进行同名格网的优化筛选;
S2、以多项式为平差模型,针对系数阵秩亏问题引入基准影像,针对像素值分布不均问题引入虚拟控制,针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束,对测区内各影像的色彩变换参数进行引入权值的整体求解;
S3、利用S2得到的色彩变换参数对每景影像做色彩校正,得到色彩均衡的影像序列。
优选的,所述S1中利用格网特征进行同名格网的优化筛选的过程包括:
采用梯度距离GD来对格网内容进行相似性度量;
利用设定的筛选准则进行筛选,所述筛选准则包括:
1)将所有同名格网按梯度距离升序进行排列,取其前部一定比例的格网;
2)计算同名格网的平均梯度值
Figure BDA0002498599160000021
Figure BDA0002498599160000022
Figure BDA0002498599160000025
的同名格网,i、j表示一对拥有重叠区域的影像各自的序号,k为影像格网序号;
3)选取梯度距离GD<0.2的同名格网;
4)按同名格网灰度均值差值的三倍中误差剔除局部异常值。
优选的,所述S1中采用基于OpenMP的多线程加速编译处理方案。
优选的,所述平差模型包括:
对于测区内的每对具有重叠关系的影像列出多个如下的误差方程:
Figure BDA0002498599160000023
式中,xi,xj为影像i和影像j位于重叠区同名地物的灰度均值,ai、bi、ci、aj、bj、cj为灰度变换的多项式系数。
优选的,所述针对系数阵秩亏问题引入基准影像的过程包括:
在测区内选取一景或多景影像作为基准影像,引入误差方程形式的基准影像方程:
Figure BDA0002498599160000024
式中,xrefer为基准影像中与影像i中对应的同名格网的平均灰度值。
优选的,所述针对像素值分布不均问题引入虚拟控制的过程包括:
对所有除基准影像外的待变换影像引入虚拟控制点,引入误差方程形式的虚拟控制方程:
v2=x*(Mrefer/Mi)-(aix2+bix+ci)
式中,x取值为1或255,(Mrefer/Mi)为基准影像灰度均值与第i幅待变换影像灰度均值的比值。
优选的,所述基准影像灰度均值Mrefer为:
Figure BDA0002498599160000031
式中,
Figure BDA0002498599160000032
为第l幅基准影像的灰度均值,Pl为第l幅基准影像的权重:
Figure BDA0002498599160000033
式中,di为待变换影像与各基准影像间的距离,λ为指数值,n为基准影像数量。
优选的,所述针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束的过程包括:
建立误差方程形式的约束公式:
Figure BDA0002498599160000034
式中,mi和mj分别为具有弱连接重叠关系的第i景和第j景影像其重叠区域的灰度均值。
优选的,所述S2引入权值的整体求解包括:
建立权矩阵如下:
Figure BDA0002498599160000035
式中,
Figure BDA0002498599160000036
为具有重叠区域的影像i和影像j间方程的权值,其值为该重叠区同名格网数量与测区内所有同名格网数量的比值;kij为影像i和影像j的重叠区域中同名格网的数量;
Figure BDA0002498599160000037
为基准影像的权值,s为基准影像方程数量;
Figure BDA0002498599160000038
Figure BDA0002498599160000039
分别为虚拟控制方程和困难地区的约束方程权值,q和h分别为虚拟控制方程和约束方程的数量。
优选的,所述S2采用迭代方法求解平差模型,参数初始值设为0,当所有参数的改正数均小于10-15时停止迭代。
与现有技术相比,本发明具有以下优点:
1、本方法针对多时相卫星遥感影像间的色彩不均衡现象,实现了基于区域网平差的卫星遥感影像色彩均衡,应用于中国大陆区域Landsat5卫星影像匀色处理中表明,经本发明方法处理后的测区内没有出现明显的亮块或色块,且各影像间的色彩过渡较好,相邻影像间的色调和亮度趋于一致,充分说明了本发明方法的可行性和有效性。
2、本方法对测区内的每景影像都建立了色彩校正参数,且充分考虑了影像重叠区间的色彩差异,利用区域网平差整体求解各影像色彩变换参数,同时可通过在测区内选取不同位置、不同数量的基准影像的方式控制测区内影像的基本色调。
3、在求基准影像均值时给予每景影像一定的权重,提高了色彩均衡处理的成功率。
附图说明
图1为本发明方法的流程示意图;
图2为实施例中同名格网筛选结果图;
图3为实施例中全国原始影像缩略图;
图4为实施例中影像增强处理前后对比图;
图5为实施例中单景影像匀光前后对比图;
图6为实施例中经匀光和影像增强处理后的全国影像缩略图;
图7为实施例中本发明方法匀色后的全国影像缩略图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施,给出了详细的实施方式和具体的操作过程,但本发明的保护范围不限于下述的实施例。
对于测区内的遥感影像序列,如果每景影像在各通道内部不存在明显的色彩和亮度差异,那么影像间的色彩差异主要由不同通道之间亮度的整体偏移产生的。理论上,位于不同影像上的同一地物应该有一致的颜色信息,然而由于影像成像时间,成像时的环境不同等方面原因,导致影像重叠区的同一地物之间产生了色彩差异现象。因此,可对影像建立适当的函数模型对其进行色彩校正,缩小其灰度差异,使各影像色彩一致或接近。
如图1所示,本申请提出一种基于区域网平差的卫星遥感影像色彩均衡方法,包括:
S1、根据影像空间位置提取影像重叠区,对影像重叠区划分均匀格网,利用HOG特征提取算法提取格网特征,利用格网特征进行同名格网的优化筛选;
S2、以多项式为平差模型,针对系数阵秩亏问题引入基准影像,针对像素值分布不均问题引入虚拟控制,针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束,对测区内各影像的色彩变换参数进行引入权值的整体求解;
S3、利用S2得到的色彩变换参数对每景影像做色彩校正,得到色彩均衡的影像序列。
HOG是一种在计算机视觉和图像处理中应用较为普遍的特征描述子,该特征描述子是基于图像局部区域的,通过计算图像局部区域的梯度信息来表示物体边缘特征,通过计算图像局部区域的方向梯度直方图来构成特征,可用于描述图像的局部纹理特征。
对于一对拥有重叠区域的影像{Ii,Ij}其重叠区域为{Qi,Qj}。{Qi,Qj}基于HOG的格网特征提取处理步骤如下:
①计算影像梯度信息
分别计算影像Qi、Qj的梯度信息。梯度信息包括梯度幅值G(x,y)和梯度方向α(x,y),两者的计算计算公式如式(1)至式(4):
梯度幅值:
Gx(x,y)=I(x+1,y)-I(x-1,y) (1)
Gy(x,y)=I(x,y+1)-I(x,y-1) (2)
Figure BDA0002498599160000051
梯度方向:
Figure BDA0002498599160000052
对于处于影像边缘的像素,会在某一侧缺少进行运算的像素。以影像左侧边缘为例,可对影像左侧扩充一列,扩充值为影像第一列像素的灰度值再进行梯度信息统计,对于其余三个边缘可采用同样方式处理。
②分格网统计影像方向梯度直方图
将影像Qi、Qj分别均匀分割成互不重叠的64*64像素的方格网,这些格网可以称为细胞单元(cell)。对每个cell将其按角度范围分为36个bin,即每个bin包含10°,将cell内的每一个像素根据其梯度方向α在36个bin中进行加权累加,权值为该点的梯度幅值。对cell内所有像素点按上述方法处理后即可得到当前cell的方向梯度直方图。
对影像Qi、Qj内的所有格网按上述方法处理,得到影像所有格网的梯度方向直方图
Figure BDA0002498599160000061
其中k为影像格网序号。
步骤S1中利用格网特征进行同名格网的优化筛选的过程包括:
(1)采用梯度距离GD来对格网内容进行相似性度量:
同名格网内容的相似性需要一个度量,需采用指标进行定量描述,本方法采用梯度距离GD来对格网内容进行相似性度量,影像{Qi,Qj}的第k个同名格网的梯度距离计算公式如式(5):
Figure BDA0002498599160000062
式中,L为HOG算法中bin的数量,本实施例L取值为36;
Figure BDA0002498599160000063
为影像{Qi,Qj}第k个同名格网的第t个bin的值。一般情况下,同名格网间的梯度距离GD越小,两者间的相似程度越高。
(2)利用设定的筛选准则进行筛选:
由于云雾条件干扰、影像的成像质量等因素,很难保证测区内所有影像是同一时期的,该测区内的影像序列一般会由不同时相的影像甚至是异源影像组成。因此很难保证重叠区域内所有格网内的地物内容是一致或相近的,如图2所示,左右影像为一对影像的部分重叠区,其影像内容由于大块云体的影响而产生了明显的差异。本申请是基于不同影像中的同一个地物拥有一致色彩信息的,如果地物内容差异较大,统计后的灰度均值意义不大,引入平差计算后会对结果造成一定的影响,因此需要对格网按照一定的筛选准则进行筛选,确定用于后续平差计算的同名格网。筛选准则包括:
1)将所有同名格网按梯度距离升序进行排列,取其前部一定比例的格网,例如取其前20%的格网;
2)计算同名格网的平均梯度值
Figure BDA0002498599160000064
Figure BDA0002498599160000065
Figure BDA0002498599160000066
的同名格网,i、j表示一对拥有重叠区域的影像各自的序号;因为平均梯度值小于等于1,表明该区域的梯度变化较为缓慢,一般为云体、水体或其他梯度特征不明显的区域,这种地区可信度和稳定度低,在处理前应该予以剔除;
3)选取梯度距离GD<0.2的同名格网;
4)即使按上述步骤进行了筛选,结果仍然可能会有局部异常值,如部分包含云体边缘区域的格网仍未剔除,在统计格网灰度均值时会出现较大偏差,会对后续处理产生较大影响。为了保证所筛选格网内容的一致性,需要对其进行进一步的剔除处理,可按同名格网灰度均值差值的三倍中误差进一步剔除局部异常值。
经上述步骤筛选过后的格网,可以认为同名格网中的内容具有高度的相似性,其灰度均值可以进入下一步的平差处理。图2的右图为采用上述方法筛选后剩余同名格网分布示意图,图中白色格网所选区域为经上述方法筛选后的结果,从图中可以看出格网已避开有云遮挡区,各同名格网的地物内容存在较高的一致性,效果较好。
基于HOG的同名格网筛选,在实际运算中涉及到大量的统计计算,如需根据地理坐标提取重叠区,计算重叠区影像的梯度值和梯度方向,并根据梯度值和梯度方向统计各个格网的梯度方向直方图,然后根据梯度方向直方图计算各同名格网的梯度距离并统计各格网灰度均值,利用梯度距离、梯度值和灰度均值差值的3倍中误差进行格网筛选,当测区内影像较多,需要进行数次循环,耗费大量时间。考虑到各次循环间不存在关联性,可以对程序采用多线程并行方式来加速计算。本实施例中,步骤S1采用基于OpenMP的多线程加速编译处理方案。
OpenMP(Open Multi-Processing)是一种可用于共享内存并行系统的多线程编程方案,可支持C、C++和Fortran等编程语言。OpenMP提供了并行算法的高级抽象描述,适用于拥有多核CPU计算机的并行编程。编译器可自动识别程序中所添加的pragma指令,自动开辟多个线程将程序并行处理,利用OpenMP进行多线程并行运算能够降低并行编程的难度和复杂度,且当编译器不支持OpenMP时,程序会退化成普通(串行)程序,程序中已有的OpenMP指令不会影响程序的正常编译运行。
为了缩小影像间的灰度差异,使各影像色彩一致或接近,需要对影像建立适当的函数模型对其进行色彩校正。假设使用的灰度变换函数为F,那么影像i和影像j之间存在以下关系:
F(xi)=F(xj) (7)
其中,xi,xj为影像i和影像j位于重叠区同名地物的灰度均值,如若采用二次多项式作为灰度变换函数,则式(7)可改写为:
Figure BDA0002498599160000081
式中,ai、bi、ci、aj、bj、cj为灰度变换的多项式系数。由于影像中的地物存在色彩不一致现象,所以两者间存在一定的差值,可列出误差方程如式(9):
Figure BDA0002498599160000082
对于测区内的每对具有重叠关系的影像都可列出多个式(9)的误差方程作为平差模型。因此,所求解方程为超定方程组,可利用最小二乘求解出参数的最优解,用每景影像的参数对其进行灰度校正,即可得到色彩均衡的影像序列。
根据测量平差原理可知,式(9)所组成的误差方程组缺少必要的起算数据,系数阵秩亏,要得到唯一解便要引入一定的基准条件。因此,可在测区内选取一景或多景质量较好,色调、反差适中且与周围影像连接较强的影像作为基准影像。因此步骤S2中针对系数阵秩亏问题引入基准影像的过程包括:在测区内选取一景或多景影像作为基准影像,引入误差方程形式的基准影像方程。
对于作为基准的影像,其多项式系数为已知的,即其二次项系数a和常数项系数c为0,一次项系数b为1,当影像l为基准影像时,其基准影像方程如式(10)所示:
Figure BDA0002498599160000083
对于基准影像而言,其格网的灰度均值是作为已知值参与运算,因此与其相连的重叠区影像方程需在式(8)的基础上进行修改,结果如式(11)所示:
Figure BDA0002498599160000084
式中,xrefer为基准影像中与影像i中对应的同名格网的平均灰度值。将其改写为误差方程形式,得到基准影像方程:
Figure BDA0002498599160000085
当影像数量较少,引入一张质量较好,色调、反差适中且与周围影像连接较强的影像作为基准影像即可得到较好的效果。而对于测区内影像较多,覆盖范围较广时,只用一幅影像作为基准控制会产生误差传递,导致边缘地区影像出现较大的色彩差异,从而致使色彩均衡处理的失败。因此,对于影像较多,覆盖范围较广的测区需要选取多景基准影像,基准影像应均匀分布在整个测区内,如测区内无法找到符合条件的影像作为基准影像,可适当采用影像增强方法制作符合要求的基准影像。
步骤S2中针对像素值分布不均问题引入虚拟控制的过程包括:对所有除基准影像外的待变换影像引入虚拟控制点,引入误差方程形式的虚拟控制方程。
由于所提取格网的灰度均值大部分都在整景影像灰度均值上下波动,会造成类似于经典平差中点分布不均的情况,即在接近影像灰度均值的灰度值能够得到较好的校正,而大于或小于影像均值的灰度值校正后会有较大偏差。因此,需要引入一定的控制条件,考虑到正常情况下一景8位影像灰度值在[1,255]之间,可对所有除基准影像外的待变换影像引入虚拟控制点,其方程如下:
aix2+bix+ci=x*(Mrefer/Mi) (13)
式中,x取值为1或255,(Mrefer/Mi)为基准影像灰度均值与第i幅待变换影像灰度均值的比值。
将式(13)改为误差方程形式,得到虚拟控制方程:
v2=x*(Mrefer/Mi)-(aix2+bix+ci) (14)
在采用多幅影像作为基准时,式(14)中的基准影像灰度均值Mrefer如果直接采用多幅影像的灰度均值的平均值会导致基准影像与周围的影像存在明显的亮度差异,从而导致色彩均衡处理的失败。因此,可在求基准影像均值时给予每景影像一定的权重,考虑到各景影像在空间中存在一定的关系,可根据待变换影像与每景基准影像间的距离赋予权重,即距离越远权重越小,距离越近权重越大。在实际操作中,可取整景影像的中心点坐标作为整景影像的空间坐标,计算待变换影像与各基准影像间的距离di,计算公式如式(15):
Figure BDA0002498599160000091
式中(Xi,Yi)和(Xl,Yl)分别为影像i和第l景基准影像中心点坐标值。
在得到待变换影像与各基准影像的距离后,可根据反距离加权法,计算各个基准影像的权重,计算公式如式(16):
Figure BDA0002498599160000092
式中,Pl为第l幅基准影像的权重,di为待变换影像与各基准影像间的距离,λ为指数值,n为基准影像数量。
在计算基准影像平均灰度均值时,各基准影像所占权重大小受参数λ的影响,即随着待变换影像与基准影像间的距离的增大,基准影像对待变换影像的影响的权重按指数规律减少,本申请中选取的λ值为1。在计算中,各基准影像对待变换影像作用的权重是成比例的,且这些权重的值的总和为1。
将各基准影像的权重带入基准影像灰度均值的计算中即可得到Mrefer,计算公式如式(17):
Figure BDA0002498599160000101
式中,
Figure BDA0002498599160000102
为第l幅基准影像的灰度均值。
基于HOG的重叠区同名格网的选取方法在实际应用中会由于影像重叠区所存在的云体及其阴影,地物类型(如大面积沙漠、水体或雪等特征不明显,梯度变化较为缓慢的区域)等因素的干扰,导致在影像重叠区提取的格网较少甚至无法提取格网的现象出现,在后续处理中影像的色调和亮度会偏向连接较强即格网数量较多的相邻影像,最终在弱连接处出现明显的色调和亮度差异,导致最终色彩均衡的失败。步骤S2针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束,约束公式如式(18)所示:
Figure BDA0002498599160000103
式中,mi和mj分别为具有弱连接重叠关系的第i景和第j景影像其重叠区域的灰度均值。
将上式改为误差方程形式的约束公式:
Figure BDA0002498599160000104
在实际处理时,可将少于一定数量同名格网的影像重叠区定义为困难地区即弱连接区域,对于此区域需要利用式(19)加以约束。
在式(9)和式(12)中,x表示在影像重叠区经筛选后的各个格网的平均灰度值。每个重叠区的格网数量代表了两景影像的连接强度,决定了该重叠区域对整个测区的影响力。由于各重叠区所提取的同名格网数量不同,各景影像之间的连接强度也不同,因此需要在求解时引入一定的权值,每个重叠区的权值大小为该重叠区同名格网数量与测区内所有同名格网数量的比值。
对于式(10)中的基准影像方程,由于是将基准影像当作已知起算数据,因此需要给所有基准影像方程赋予的权值为1。
对于式(14)中的虚拟控制方程,为了保证匀色结果,其赋予的权值应小于基准影像,且大于式(9)和式(12)的权值。
对于式(19)中针对弱连接区的约束方程,其权值可与虚拟控制方程的权值相等。
根据上述条件,S2中建立权矩阵如下:
Figure BDA0002498599160000111
式中,
Figure BDA0002498599160000112
为具有重叠区域的影像i和影像j间方程的权值,其值为该重叠区同名格网数量与测区内所有同名格网数量的比值;kij为影像i和影像j的重叠区域中同名格网的数量;
Figure BDA0002498599160000113
为基准影像的权值,s为基准影像方程数量;
Figure BDA0002498599160000114
Figure BDA0002498599160000115
分别为虚拟控制方程和困难地区的约束方程权值,q和h分别为虚拟控制方程和约束方程的数量。
本实施例中,步骤S2采用迭代方法求解平差模型,对测区内所有影像按式(9)、式(12)、式(14)和式(19)列出误差方程组,并将其改写为矩阵形式:
V=BX-L (21)
在求解中还需将基准影像的条件方程加入其中,其形式如式(22)所示:
V=[B Cl]TX-[L Ll]T (22)
当影像l为基准影像时,其Cl为3×3的单位矩阵,Ll=[0 1 0]T
根据最小二乘原理,可得参数最优解为X=(BTPB)-1BTPL,其中P为式(20)中确定的权阵。在进行迭代求解时,参数初始值可设为0,当所有参数的改正数均小于10-15时可停止迭代,即可得到各景影像的校正参数a、b和c。最后利用所得参数对每景影像进行多项式校正即可得到匀色后的结果影像。
实施例
本实施例针对覆盖中国大陆的Landsat5卫星遥感影像(波段B1、B2和B3)进行基于区域网平差的影像间色彩一致性处理实验。由于部分影像存在严重的影像质量问题,我国部分地区的影像缺失,经筛选后剩余513景影像,每景影像的大小为7500×8500像素,其地面分辨率为30m,影像成像时间为2006年至2010年之间,所有影像均已做过辐射校正和几何校正。未作任何影像增强和匀光匀色处理的原始影像如图3所示,从图中可以看出大量影像显示出低对比度,且影像呈现整体偏暗的情况势必会影像后续的影像匀光匀色处理,视觉效果极差,影像间存在明显的亮度和色彩不一致现象,需要对其进行进一步的处理。
一、影像增强处理:
针对图3中大量影像呈现低对比度、亮度过暗的情况,需要对其进行影像增强处理,本实施例主要采用1%线性裁剪拉伸进行增强处理,部分影像增强前后对比如图4所示,经影像增强后的影像视觉效果有明显增强,影像整体亮度有了提升,更多细节信息得以显现,为下一步的匀光匀色处理提供保障。
二、基于Mask算法的单景影像匀光处理:
单景影像的色彩一致性处理,主要是为了解决单幅影像内部出现的亮度分布不均和色调差异,它与影像间的色彩均衡处理都是影像产品生产过程中的重要组成部分。本申请方法需要保证单景影像内部不存在明显的亮度和色彩不均衡现象,因此需要采用单景影像匀光算法对影像进行匀光处理。
Mask匀光算法是一种可行且有效的单景影像匀光算法,Mask匀光的主要思想是将光照不均匀的影像看作是光照均匀的影像叠加了一个光照不均匀的背景影像的结果,可利用高斯低通滤波估计背景影像,通过影像相减操作从原始影像中去除背景影像,并对结果做一定的拉伸处理即可得到光照均匀的结果影像。
因此,本实施例对所有影像做了基于Mask算法的单景影像匀光处理,部分影像处理效果如图5所示。从图中可以看出在经单景影像匀光处理前,各影像内部均存在较为明显的亮度和色彩不均衡现象,而经Mask匀光算法处理后影像在保持原有影像色调的基础上其亮度和色调分布更加均衡,效果较好。
三、基于区域网平差的色彩均衡处理:
基于HOG的同名格网筛选:
对影像进行增强和匀光后,即可开始进行本申请基于区域网平差的色彩均衡处理。首先,需要对影像重叠区进行同名格网的筛选,在提取格网时采用基于OpenMP的多线程加速,由于内存限制,总共开辟9个线程进行加速处理,各波段提取格网平均耗时277.264min。
重叠区格网筛选情况如表1所示,全国513景影像间共有1369个重叠关系,但是由于影像重叠区所存在的云体及其阴影,地物类型(如大面积沙漠、水体或雪等特征不明显,梯度变化较为缓慢的区域)等因素的干扰,导致各波段有大量重叠区未提取到格网,因此在后续处理中需要利用影像重叠区的整体灰度均值关系加以约束,同时在处理中将少于20个同名格网的重叠区定义为弱连接区,同样采取上述方法处理。在处理中由于部分地区影像重叠区相似程度较高,提取的格网数量远多于所需数量,可以对其进行一定限制,本实施例限制每个影像重叠区至多有500个同名格网。
表1各波段同名格网筛选情况统计表
Figure BDA0002498599160000131
基准影像的选取:
基准影像的选取在匀色过程中至关重要,基准影像的色调决定了它周围影像的色调,选取的基准影像应是成像质量较好,色调、反差适中且与周围影像连接较强的影像,整个测区中基准影像不应过多,且需要均匀分布于整个测区,选取的基准影像还应考虑到周围影像的地物类型。如在图3中西部地区的地物多为荒漠,植被覆盖较少的裸地,其色调以土黄色为主,而东南地区的地物植被覆盖较多,影像应以绿色基调为主,在上述地区都应挑选质量较好的一景影像作为基准影像,以保证匀色后荒漠地区不会出现绿色基调和植被较多地区出现黄色基调的现象。本实施例中,按上述方式在全国范围内选取了均匀分布的10景影像作为基准影像,并将其带入后续处理中。
匀色结果:
将上述得到的同名格网灰度均值、基准影像、弱连接重叠区的灰度均值按本申请方法列误差方程并进行迭代求解,在得到每景影像的灰度变换参数后,逐影像分通道进行灰度变换,得到最终匀色结果。
对图3中的所有原始影像进行影像增强处理和单景影像匀光处理,结果如图6所示,从图中可以看出,在各单景影像内部的色调和亮度较为均匀,但整个测区内有明显的亮度差异和色彩差异,影像呈现“马赛克”状,视觉效果极差,尤其是在我国东南地区和东北地区,每景影像都存在着巨大的色彩差异,在西部地区虽然每景影像色调类似,但也存在着亮度不一的情况,严重影响了整个测区的色彩质量。
针对图6中的情况,采用本申请方法处理过后的结果如图7所示,从图中可以看出,经处理后,整个测区没有出现明显的亮块或色块,且各影像间的色彩过渡较好,测区内相邻影像色调和亮度趋于一致,充分说明了本申请方法的可行性和有效性。

Claims (10)

1.一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,包括:
S1、根据影像空间位置提取影像重叠区,对影像重叠区划分均匀格网,利用HOG特征提取算法提取格网特征,利用格网特征进行同名格网的优化筛选;
S2、以多项式为平差模型,针对系数阵秩亏问题引入基准影像,针对像素值分布不均问题引入虚拟控制,针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束,对测区内各影像的色彩变换参数进行引入权值的整体求解;
S3、利用S2得到的色彩变换参数对每景影像做色彩校正,得到色彩均衡的影像序列。
2.根据权利要求1所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述S1中利用格网特征进行同名格网的优化筛选的过程包括:
采用梯度距离GD来对格网内容进行相似性度量;
利用设定的筛选准则进行筛选,所述筛选准则包括:
1)将所有同名格网按梯度距离升序进行排列,取其前部一定比例的格网;
2)计算同名格网的平均梯度值
Figure FDA0002498599150000012
Figure FDA0002498599150000014
Figure FDA0002498599150000013
的同名格网,i、j表示一对拥有重叠区域的影像各自的序号,k为影像格网序号;
3)选取梯度距离GD<0.2的同名格网;
4)按同名格网灰度均值差值的三倍中误差剔除局部异常值。
3.根据权利要求1所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述S1中采用基于OpenMP的多线程加速编译处理方案。
4.根据权利要求1所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述平差模型包括:
对于测区内的每对具有重叠关系的影像列出多个如下的误差方程:
Figure FDA0002498599150000011
式中,xi,xj为影像i和影像j位于重叠区同名地物的灰度均值,ai、bi、ci、aj、bj、cj为灰度变换的多项式系数。
5.根据权利要求4所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述针对系数阵秩亏问题引入基准影像的过程包括:
在测区内选取一景或多景影像作为基准影像,引入误差方程形式的基准影像方程:
Figure FDA0002498599150000021
式中,xrefer为基准影像中与影像i中对应的同名格网的平均灰度值。
6.根据权利要求5所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述针对像素值分布不均问题引入虚拟控制的过程包括:
对所有除基准影像外的待变换影像引入虚拟控制点,引入误差方程形式的虚拟控制方程:
v2=x*(Mrefer/Mi)-(aix2+bix+ci)
式中,x取值为1或255,(Mrefer/Mi)为基准影像灰度均值与第i幅待变换影像灰度均值的比值。
7.根据权利要求6所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述基准影像灰度均值Mrefer为:
Figure FDA0002498599150000022
式中,
Figure FDA0002498599150000023
为第l幅基准影像的灰度均值,Pl为第l幅基准影像的权重:
Figure FDA0002498599150000024
式中,di为待变换影像与各基准影像间的距离,λ为指数值,n为基准影像数量。
8.根据权利要求6所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述针对呈现弱连接的困难区域采用重叠区的整体灰度均值加以约束的过程包括:
建立误差方程形式的约束公式:
Figure FDA0002498599150000025
式中,mi和mj分别为具有弱连接重叠关系的第i景和第j景影像其重叠区域的灰度均值。
9.根据权利要求8所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述S2引入权值的整体求解包括:
建立权矩阵如下:
Figure FDA0002498599150000031
式中,
Figure FDA0002498599150000032
为具有重叠区域的影像i和影像j间方程的权值,其值为该重叠区同名格网数量与测区内所有同名格网数量的比值;kij为影像i和影像j的重叠区域中同名格网的数量;Ps refer为基准影像的权值,s为基准影像方程数量;
Figure FDA0002498599150000033
Figure FDA0002498599150000034
分别为虚拟控制方程和困难地区的约束方程权值,q和h分别为虚拟控制方程和约束方程的数量。
10.根据权利要求1所述的一种基于区域网平差的卫星遥感影像色彩均衡方法,其特征在于,所述S2采用迭代方法求解平差模型,参数初始值设为0,当所有参数的改正数均小于10-15时停止迭代。
CN202010425601.2A 2020-05-19 2020-05-19 一种基于区域网平差的卫星遥感影像色彩均衡方法 Active CN112233190B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010425601.2A CN112233190B (zh) 2020-05-19 2020-05-19 一种基于区域网平差的卫星遥感影像色彩均衡方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010425601.2A CN112233190B (zh) 2020-05-19 2020-05-19 一种基于区域网平差的卫星遥感影像色彩均衡方法

Publications (2)

Publication Number Publication Date
CN112233190A true CN112233190A (zh) 2021-01-15
CN112233190B CN112233190B (zh) 2023-04-07

Family

ID=74111677

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010425601.2A Active CN112233190B (zh) 2020-05-19 2020-05-19 一种基于区域网平差的卫星遥感影像色彩均衡方法

Country Status (1)

Country Link
CN (1) CN112233190B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112884676A (zh) * 2021-03-18 2021-06-01 国家海洋信息中心 一种基于分空间渐进控制的大范围航空遥感影像调色方法
CN113409366A (zh) * 2021-06-29 2021-09-17 中国科学院空天信息创新研究院 一种遥感卫星影像的辐射校正方法及装置
CN113421193A (zh) * 2021-05-05 2021-09-21 桂林理工大学 均值-方差成本函数最小模型的多幅影像镶嵌辐射均衡方法
CN113899386A (zh) * 2021-09-27 2022-01-07 武汉大学 基于立体基准网的多源光学卫星遥感影像协同区域网平差方法及系统
CN114841881A (zh) * 2022-04-29 2022-08-02 武汉大学 基于辐射云控制的全球卫星影像自然颜色恢复方法及系统
CN116597184A (zh) * 2023-07-11 2023-08-15 中国人民解放军63921部队 最小二乘影像匹配方法
CN117670732A (zh) * 2023-11-28 2024-03-08 武汉大学 基于伯恩斯坦多项式的大场景无人机影像匀光匀色方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009010636A (ja) * 2007-06-27 2009-01-15 Sony Deutsche Gmbh 適応ヒストグラム等化方法及び適応ヒストグラム等化装置
CN107480727A (zh) * 2017-08-28 2017-12-15 荆门程远电子科技有限公司 一种sift和orb相结合的无人机影像快速匹配方法
CN107563964A (zh) * 2017-08-22 2018-01-09 长光卫星技术有限公司 大面阵亚米级夜景遥感影像的快速拼接方法
CN110599424A (zh) * 2019-09-16 2019-12-20 北京航天宏图信息技术股份有限公司 影像自动匀色处理方法及装置、电子设备、存储介质

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009010636A (ja) * 2007-06-27 2009-01-15 Sony Deutsche Gmbh 適応ヒストグラム等化方法及び適応ヒストグラム等化装置
CN107563964A (zh) * 2017-08-22 2018-01-09 长光卫星技术有限公司 大面阵亚米级夜景遥感影像的快速拼接方法
CN107480727A (zh) * 2017-08-28 2017-12-15 荆门程远电子科技有限公司 一种sift和orb相结合的无人机影像快速匹配方法
CN110599424A (zh) * 2019-09-16 2019-12-20 北京航天宏图信息技术股份有限公司 影像自动匀色处理方法及装置、电子设备、存储介质

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
卢其剑 等: "基于区域网平差的遥感影像色彩均衡算法", 《东华理工大学学报(自然科学版)》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112884676A (zh) * 2021-03-18 2021-06-01 国家海洋信息中心 一种基于分空间渐进控制的大范围航空遥感影像调色方法
CN113421193A (zh) * 2021-05-05 2021-09-21 桂林理工大学 均值-方差成本函数最小模型的多幅影像镶嵌辐射均衡方法
CN113421193B (zh) * 2021-05-05 2022-09-30 桂林理工大学 均值-方差成本函数最小模型的多幅影像镶嵌辐射均衡方法
CN113409366A (zh) * 2021-06-29 2021-09-17 中国科学院空天信息创新研究院 一种遥感卫星影像的辐射校正方法及装置
CN113409366B (zh) * 2021-06-29 2022-11-22 中国科学院空天信息创新研究院 一种遥感卫星影像的辐射校正方法及装置
CN113899386A (zh) * 2021-09-27 2022-01-07 武汉大学 基于立体基准网的多源光学卫星遥感影像协同区域网平差方法及系统
CN113899386B (zh) * 2021-09-27 2023-11-21 武汉大学 基于立体基准网的多源光学卫星遥感影像协同区域网平差方法及系统
CN114841881A (zh) * 2022-04-29 2022-08-02 武汉大学 基于辐射云控制的全球卫星影像自然颜色恢复方法及系统
CN114841881B (zh) * 2022-04-29 2024-10-15 武汉大学 基于辐射云控制的全球卫星影像自然颜色恢复方法及系统
CN116597184A (zh) * 2023-07-11 2023-08-15 中国人民解放军63921部队 最小二乘影像匹配方法
CN116597184B (zh) * 2023-07-11 2023-09-22 中国人民解放军63921部队 最小二乘影像匹配方法
CN117670732A (zh) * 2023-11-28 2024-03-08 武汉大学 基于伯恩斯坦多项式的大场景无人机影像匀光匀色方法

Also Published As

Publication number Publication date
CN112233190B (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
CN112233190B (zh) 一种基于区域网平差的卫星遥感影像色彩均衡方法
Gruszczyński et al. Comparison of low-altitude UAV photogrammetry with terrestrial laser scanning as data-source methods for terrain covered in low vegetation
CN109934154B (zh) 一种遥感影像变化检测方法及检测装置
CN111091502A (zh) 遥感影像的匀色方法及系统、存储介质、电子设备
CN112307901B (zh) 一种面向滑坡检测的sar与光学影像融合方法及系统
US20050175253A1 (en) Method for producing cloud free and cloud-shadow free images
CN114881620B (zh) 一种基于卫星遥感的国土空间监测方法与系统
CN112733596A (zh) 一种基于中高空间分辨率遥感影像融合的森林资源变化监测方法及应用
CN117765051B (zh) 一种园林绿化养护监测预警系统及方法
CN116310800B (zh) 一种基于深度学习的梯田自动提取方法及装置
CN111192298B (zh) 一种夜光遥感影像相对辐射校正方法
KR20210096925A (ko) 대용량 항공정사영상의 효율적인 색상보정 방법
CN113935917B (zh) 一种基于云图运算和多尺度生成对抗网络的光学遥感影像薄云去除方法
CN113888416A (zh) 卫星遥感图像数据的处理方法
CN116309110A (zh) 一种基于轻量化深度神经网络的低光照图像去雾方法
CN117115669B (zh) 双条件质量约束的对象级地物样本自适应生成方法及系统
CN114663771A (zh) 一种基于分区分层理论的山区耕地智能化提取方法
CN114331937A (zh) 基于反馈迭代调节的低光照条件下多源图像融合方法
CN117437489A (zh) 基于决策树模型的城市绿地提取方法
CN116703744B (zh) 一种基于卷积神经网络的遥感影像匀光匀色方法和装置
KR102397148B1 (ko) 대용량 항공정사영상의 저해상도 컬러영상을 이용한 색상보정 방법
CN115689941A (zh) 一种跨域生成对抗的sar影像补偿方法及计算机可读介质
CN113781587B (zh) 基于最佳路径的遥感影像色彩一致性处理方法
CN113743373A (zh) 基于深度学习的高分遥感影像耕地变化检测装置及方法
CN110322454B (zh) 一种基于光谱差异最大化的高分辨遥感图像多尺度分割优化方法

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