CN109671111B - 基于可见光遥感图像的温度场调制方法 - Google Patents

基于可见光遥感图像的温度场调制方法 Download PDF

Info

Publication number
CN109671111B
CN109671111B CN201811333483.1A CN201811333483A CN109671111B CN 109671111 B CN109671111 B CN 109671111B CN 201811333483 A CN201811333483 A CN 201811333483A CN 109671111 B CN109671111 B CN 109671111B
Authority
CN
China
Prior art keywords
modulation
image
shadow
temperature
temperature field
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
Application number
CN201811333483.1A
Other languages
English (en)
Other versions
CN109671111A (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.)
Xidian University
Original Assignee
Xidian 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 Xidian University filed Critical Xidian University
Priority to CN201811333483.1A priority Critical patent/CN109671111B/zh
Publication of CN109671111A publication Critical patent/CN109671111A/zh
Application granted granted Critical
Publication of CN109671111B publication Critical patent/CN109671111B/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/40Analysis of texture
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • G06F18/232Non-hierarchical techniques
    • G06F18/2321Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/40Image enhancement or restoration by the use of histogram techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/20Image preprocessing
    • G06V10/25Determination of region of interest [ROI] or a volume of interest [VOI]
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明属于红外技术领域,具体涉及基于可见光遥感图像的温度场调制方法,包括:根据可见光遥感图像获取红外纹理调制模板;对所述红外纹理调制模板进行第一调制处理得到初级红外纹理温度场;对所述调制图像进行第二调制处理得到红外纹理温度场图像。本发明实施例,显著增强可见光遥感图像增质后的质量,使得调整生成的红外纹理温度场真实性大大提高;提高红外纹理温度场逼真度;使得调制更加符合实际情况,更加准确。

Description

基于可见光遥感图像的温度场调制方法
技术领域
本发明属于图像处理领域,具体涉及基于可见光遥感图像的温度场调制方法。
背景技术
红外成像系统具有精度高、抗干扰能力强、使用灵活等特点,越来越受到各国的关注和倡导。红外成像系统的研制往往要考虑各种性能指标,找到能够克服时间、环境、地域的限制,降低成本,缩短周期,同时生成各种环境下的高真实感目标仿真图像,成为一个迫切的要求。
现有的温度场调制生成方法中,虽然进行了固化阴影的去除,但未进行运动物体的去除。固化的运动物体遮盖所处位置的材质信息,导致场景仿真的灵活性和真实性受到限制,还会给后续的地物材质分割和红外纹理温度场的调制带来不便。
现有的温度场调制生成研究中,所有材质类型只使用一种调制模型,调制方式单一,未根据材质类型的特点选择适应的调制模型。导致生成的红外纹理温度场逼真度不够高,且在调制一些材质区域时会发生与实际不符的现象。
现有的红外纹理温度场调制方法,在处理大规模场景时,仅提取了高程信息中的海拔高度信息,进行红外纹理温度场的调制,未考虑太阳朝向带来的温度变化。本发明不仅提取了高程数据中的海拔高度信息,而且还提取了高程数据蕴含的太阳高度角、太阳方位角、经纬度等信息。首先利用原来没有考虑朝向的温度场分布,结合高程数据确定出不同时刻的太阳角度。然后确定出不同时刻下阳光照射比例。最后利用阳光照射比例信息,对红外纹理温度场进行进一步细节调制。
发明内容
为了解决现有技术中存在的上述问题,本发明提供了一种基于可见光遥感图像的温度场调制方法。本发明要解决的技术问题通过以下技术方案实现:
本发明实施例提供了基于可见光遥感图像的温度场调制方法,步骤如下:
根据可见光遥感图像获取红外纹理调制模板;
对所述红外纹理调制模板进行第一调制处理得到调制图像;
对所述调制图像进行第二调制处理得到红外纹理温度场图像。
在本发明的一个实施例中,对所述可见光遥感图像获取红外纹理调制模板,包括,
根据所述可见光遥感图像得到阴影区域;
对所述阴影区域进行补偿,得到补偿图像;
去除所述补偿图像中的运动物体,得到红外纹理调制模板。
在本发明的一个实施例中,对所述红外纹理调制模板进行第一调制处理得到调制图像,包括,
对红外文理调制模板进行分割处理与材质识别,得到分割图像;
根据分割图像建立热特征预测模型;
将所述热特征预测模型的宏观温度与材质区域分布图进行一一映射,得到基本红外纹理;
对所述宏观温度分布图进行分类调制得到调制图像。
在本发明的一个实施例中对所述宏观温度分布图进行分类调制得到调制图像,包括,
对所述基本红外纹理进行平坦表面调制,得到初级分类调制模型;
对所述初级分类调制模型进行起伏表面调制,得到分类调制模型;
对所述分类调制模型进行分析得到调制图像。
在本发明的一个实施例中,对所述调制图像进行第二调制处理得到红外纹理温度场图像,包括,
对调制图像进行海拔高度调制,得到海拔高度调制模型;
对所述海拔高度调制模型进行阳光照射比例调制得到红外纹理温度场。
在本发明的一个实施例中,对所述海拔高度调制模型进行阳光照射比例调制得到红外纹理温度场,包括,
对海拔高度调制模型进行山体阴影计算,得到初级阳光照射比调制模型;
对初级阳光照射比调制模型进行阳光照射比例调制,得到红外纹理温度场。
与现有技术相比,本发明的有益效果:
(1)本发明采用先对可见光遥感图像预处理获得高质量遥感图像,显著增强可见光遥感图像增质后的质量,使得调整生成的红外纹理温度场真实性大大提高;再将高质量遥感图像作为红外纹理的调制模板进行分类细节调制,提高红外纹理温度场逼真度;最后分析新图像中蕴含的环境因素对图像进行基于海拔高度和阳光照射比的进一步调制,使得调制更加符合实际情况,更加准确;
(2)本发明实施例中将图像块修复算法应用到运动物体去除中,并结合阴影剔除算法,防止遥感图像中阴影、运动的物体导致纹理固化;
(3)本发明实施例中采用地物材质分割将图像分割成多个感兴趣的区域,并建立热特征模型,通过在宏观均值温度上增加温度扰动,可以使物体表面辐射分布发生改变,使得红外纹理的细节增强,再对生成的红外纹理温度场进行分类调制,通过建立多种红外纹理温度场调制模型的方式,生成较为真实的红外纹理温度场;
(4)本发明实施例中采用分别建立平坦表面和起伏表面的调制模型,并对其结果进行分析,排除复杂地貌对调制准确度的影响;
(5)本发明实施例中同时对红外纹理温度场的海拔高度与阳光照射比例进行调制,提高了生成红外纹理温度场的真实性;
(6)本发明实施例中对图像中阴影进行了计算并基于阳光照射比例调制红外纹理温度场,提高了生成红外纹理温度场的真实性与准确性。
附图说明
图1为本发明实施例提供的基于可见光遥感图像的温度场调制方法的过程框架图;
图2为本发明实施例提供的基于可见光遥感图像的温度场调制方法的修复块O(p)的修复过程;
图3为本发明实施例提供的基于可见光遥感图像的温度场调制方法的可见光遥感图像去除运动物体对比图;
图4为本发明实施例提供的基于可见光遥感图像的温度场调制方法的红外纹理温度场调制过程框架图;
图5为本发明实施例提供的基于可见光遥感图像的温度场调制方法的材质分割结果图;
图6为本发明实施例提供的基于可见光遥感图像的温度场调制方法的材质宏观温度分布图;
图7为本发明实施例提供的基于可见光遥感图像的温度场调制方法的基于可见光遥感图像纹理的分类调制框架;
图8为本发明实施例提供的基于可见光遥感图像的温度场调制方法的基于可见光遥感图像的红外纹理温度场调制结果对比图;
图9为本发明实施例提供的基于可见光遥感图像的温度场调制方法的红外纹理温度场调制结果对比图;
图10为本发明实施例提供的基于可见光遥感图像的温度场调制方法的高程数据;
图11为本发明实施例提供的山体阴影计算公式中各参数的空间几何关系表示;
图12为本发明实施例提供的以一个像元为中心的3*3像素领域表示图;
图13为本发明实施例提供的可见光遥感图像与其各像素点的山体阴影值对比图;
图14为本发明提供的基于可见光遥感图像的温度场调制方法的流程框图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述,但本发明的实施方式不限于此。
如图1和图14所示,本发明实施例提供了基于可见光遥感图像的温度场调制方法,步骤如下:
根据可见光遥感图像获取红外纹理调制模板;
对所述红外纹理调制模板进行第一调制处理得到调制图像;
对所述调制图像进行第二调制处理得到红外纹理温度场图像。
在本发明的一个实施例中,对所述可见光遥感图像获取红外纹理调制模板,包括:
根据所述可见光遥感图像得到阴影区域:在HSV色彩空间,图像中的阴影区域与非阴影区域,在色度和饱和度通道表现出很大的差异性。利用这一点,将可见光遥感图像从RGB空间转换到HSV空间,可以更好的提取出阴影区域。
(1)转换空间:将RGB图像转换到HSV空间,把亮度信息与色调信息区分开来,由RGB空间转换到HSV空间的转换公式为:
Figure GDA0004230063140000031
Figure GDA0004230063140000032
Figure GDA0004230063140000033
其中,R、G、B为RGB空间的三通道,R为红通道,G为绿通道,B为蓝通道,H、S、V为HSV空间的三通道,H为色度,S为饱和度,V为颜色信息,即I,α为归一化因子。
(2)降低亮度通道值:读取色度分量H与亮度分量I,将(I+1)/(H+1)的值赋给亮度I;
(3)归一化直方图:统计0~i灰度级的像素个数占总像素个数的比例w0,统计这些灰度级像素的平均灰度μ0;统计i~255灰度级的像素个数占总像素个数的比例w1,统计这些灰度级像素的平均灰度μ1
(4)计算方差:计算0~i灰度级的像素和i~255灰度级的像素方差g,g=w0*w1*(μ01)2
(5)迭代结束:i从0开始迭代,i∈[0~255],当i>255时,迭代结束,否则,返回步骤(3)继续进行计算;
(6)分割图像:找出方差中的最大值gmax对应的i值作为全局阈值输出,大于阈值的区域为非阴影区域,小于阈值的区域为阴影区域。
(7)获得阴影区域:利用二值形态学算法,对分割的区域进行腐蚀或膨胀计算,最终得到阴影区域。
对所述阴影区域进行补偿,得到补偿图像:
对阴影检测后将对其进行补偿,该方案采用Finlayson和Trezzi提出的颜色恒常性计算方法,该方法要求图像中的阴影区域和非阴影区域满足Minkowski范式,
Figure GDA0004230063140000041
其中,e为光源颜色,f为图像三色通道的灰度值,k为比例系数,χ∈[1,∞)且为整数,属于权重因子,M为图像长度,N为图像宽度。阴影补偿基本步骤为:
(1)RGB通道处理:将RGB图像分为三个通道进行一一处理,将被作为Minkowski范式输入的各通道的灰度值;
(2)分割阴影和非阴影区域;利用阴影检测区域,将各通道的灰度图像分为阴影区域与非阴影区域;
(3)计算光源颜色:将各通道的阴影区域与非阴影区域利用Minkowski范式进行计算,得出各通道的阴影区域与非阴影区域的光源颜色;
(4)全部恢复到标准光照条件下:阴影区域恢复到标准光照条件下的表达式为:
Figure GDA0004230063140000042
其中,yb(i,j)为某通道阴影区域恢复到标准光照条件下的灰度值,y(i,j)为该通道阴影区域的灰度值,e1表示阴影区域的光照颜色。
非阴影区域恢复到标准光照条件下的表达式可表示为:
Figure GDA0004230063140000043
其中,gb(i,j)为某通道非阴影区域恢复到标准光照条件下的灰度值,g(i,j)为该通道非阴影区域的灰度值,e2表示非阴影区域的光照颜色。
在标准光照条件下,yb(i,j)和gb(i,j)的场景的光照条件一样,所以两者的光源颜色也应该相等。
(5)恢复到非阴影光照条件下:标准光照条件恢复到非阴影区域的表达式可表示为:g(i,j)=gb(i,j)*e2
则阴影区域恢复到非阴影区域的灰度值yg(i,j)应满足:
Figure GDA0004230063140000051
(6)得到补偿图像:将各通道的阴影区域恢复到非阴影区域的灰度值yg(i,j)保持非阴影区域g(i,j)值不变,赋给各通道对应的区域,即可得到各通道的灰度值。然后,将各通道按原来的R,G,B通道输出,即可得到补偿图像。
在实际计算中,只利用推导出的式对阴影区域进行补偿既可,由于式中只包含阴影区域的灰度值和两个光源颜色,可以避免破坏掉非阴影区域的色彩信息和纹理细节。
去除所述补偿图像中的运动物体,得到红外纹理调制模板:
基于图像块修复的去除算法思想为:1.利用一定的方法,将图像分成待修复区域ζ和原图像区域ψ,边界线为
Figure GDA0004230063140000052
边界线/>
Figure GDA0004230063140000053
上任意一点以p为中心的邻域为O(p),以O(p)为修复块进行修复。计算各点邻域的优先级,得到修复块修复的先后顺序;2.找出与待修复块O(p)相似的原图像区域块ψ(p'),进行邻域待修复块O(p)的像素填充,邻域O(p)里的修复区域也有了像素内容,如图2所示。3.不断进行上述步骤,待修复区域逐渐减少,最终完成全部的待修复区域的修复。具体的修复步骤如下:
(1)获得优先修复块
优先级直接决定图像的修补顺序,它与待修复区域内所含像素内容有关。根据Criminisi算法可知,边界线
Figure GDA0004230063140000054
上任意一点以p为中心的邻域O(p),其优先权P(p)可表述为:P(p)=D(p)C(p),其中,C(p)是表示置信度的物理量,描述邻域O(p)内所有像素点的置信度,满足/>
Figure GDA0004230063140000055
S(O(p))表示邻域O(p)的面积。
在原图像区域ψ的像素点置信度C(p)的初始值设为1,在待修复区域ζ的像素点置信度C(p)的初始值设为0。O(p)内的所有像素点的置信度C(p),也可以说是p点周围的可靠信息的多少。当邻域O(p)大部分像素处在原图像区域ψ,则p点周围的可靠信息越多,p点的像素内容越容易被确定,其置信度C(p)越高,反之,邻域O(p)在待修复区域占的比重大,则p点周围的可靠信息少,p的像素内容越不容易确定。
其中,D(p)是与邻域O(p)结构有关的物理量满足
Figure GDA0004230063140000061
其中,/>
Figure GDA0004230063140000062
表示的是沿p点的等照度线的方向向量,即垂直于亮度梯度方向,np为p点沿边界线/>
Figure GDA0004230063140000063
切线的单位法向量,α是归一化项。
修复块的先后顺序能直接影响修复的质量。通过上面确定的优先权P(p),虽然能够确定待修复块的先后顺序,但是,
Figure GDA0004230063140000064
为0时,也即/>
Figure GDA0004230063140000065
与np垂直时,数据项D(p)为0,P(p)直接就等于0了,导致该修复区域未被修复到,修复质量下降。且在后续的更新置信度的过程中,经过反复的迭代,置信度也会降低为0,导致无法确定待修复区域的修复顺序。通过分析优先权公式可得,优先修复结构特征很强的邻域,也即对应的D(p)比较大的邻域能够保证修复方向,确定待修复区域的先后顺序。
确定以D(p)的大小为准的先后修复顺序,可以在原有的Criminisi算法上增加D(p)的权重,改进的优先权公式为:P(p)=D(p)C(p)+w×D(p)利用上式对以边界线上的像素点为中心的待修补区域计算得到P(p),P(p)值最大的待修复块就是优先选择的修复块。
(2)建立自适应的修复块的尺寸模型
优先的修补块一般采用尺寸为9×9的固定尺寸进行修复,但是图像中各种材质的纹理结构类型非常丰富,有些材质的纹理很平滑,在这些材质区域如果采用过小的修复尺寸,会导致块效应的产生;有些材质的纹理结构很复杂,具有较强的颜色纹理,在这些材质区域如果采用过大的修复尺寸,会导致边界效应。为了能够避免这种因修复尺寸的大小带来的影响,本文采用自适应的根据材质纹理结构匹配尺寸模型,对修复尺寸进行适当调整。
图像中的颜色信息反映其纹理结构,纹理结构信息可以用其像素点颜色的梯度▽I(q)来表示。本文利用修复区域所处的材质区域的像素点的梯度来自适应的选择修复块的尺寸CC,具体如下:
Figure GDA0004230063140000066
(3)最佳的匹配块修复
边界线
Figure GDA0004230063140000068
上的任意一点p邻域O(p)通过优先权公式算出所有优先值。通过得到的最大优先值对该领域进行先行修复,设优先值最大的是以特定点p1为中心的待修复块O(p1),通过计算待修复块O(p1)处于的原图像区域的材质纹理像素点梯度/>
Figure GDA0004230063140000067
得到其对应的修复尺寸,再从原图像区域ψ搜索出与其纹理结构特征匹配最佳的区域块,对待修复块O(p1)进行填充修复。待修复块O(p1)与原图像区域块ψ(p1')的匹配关系如下:Ω=argmin(d(O(p1),ψ(p1')))其中,d(O(p1),ψ(p1')为待修复块O(p1)与原图像区域块ψ(p1')的像素差距,用其像素颜色的平方差和(SSD)计算获得:
Figure GDA0004230063140000071
其中,n为邻域O(p1)中的原图像区域的像素个数,/>
Figure GDA0004230063140000072
为邻域O(p1)中的p1点变换到Lab色彩空间的各通道的值,/>
Figure GDA0004230063140000073
为原图像区域块ψ(p1')中的p1'点变换到Lab色彩空间的各通道的值,其中变换的目的在于为了得到更好的视觉效果的图像。从原图像区域中搜索遍历到距离最小的匹配块,对待修复块O(p1)进行修复填充。
(4)更新置信度和优先级
通过最佳匹配块修复了待修复块O(p1),邻域中的待修复区域修复完成。待修复区域和原图像区域发生变化,这些新的完成修复的像素点的置信度也发生了变化。即:
Figure GDA0004230063140000074
新修复完成的像素点置信度等于p1点的置信度。在不断的迭代过程中,置信度不断减小,越靠近边界线/>
Figure GDA0004230063140000075
的像素点置信度越高,越远离边界线/>
Figure GDA0004230063140000076
的像素点置信度越低。对于处理遥感图像中的杂物去除运用中,杂物一般属于小物体,所以,在进行杂物去除时,不断的迭代过程造成的置信度减小问题,可以忽略。
填充修复了待修复块O(p1),图像的修复区域发生变化,边界线也发生了变化,在新的边界线上进行新的优先级的计算,距离较为远的点,置信度不变,只需要计算被新修复区域影响到的像素点,更新其待修复块的优先级即可。重复上述步骤,直到修复完成。
图3为可见光遥感图像去除运动物体的结果图,其中(a)(b)为两源遥感图像,(a’)(b’)分别为对应的遥感图像去除运动物体的结果。可以看出,该运动物体去除方法,很好的去除了图像中运动的汽车,得到更高质量的红外纹理调制模板,为后续的遥感图像的分割、温度场的调制和红外场景的仿真,奠定良好的基础。
在本发明的一个实施例中,如图4所示,对所述红外纹理调制模板进行第一调制处理得到调制图像,包括:
对红外文理调制模板进行分割处理与材质识别,得到分割图像:
图像分割通常处理的是像素边缘和物体的边界。可见光遥感图像中不同地物的温度因材质的差异表现出明显的不同,必须进行材质分割得到不同地物的温度映射。目前,分割遥感图像的算法有基于阈值、基于区域生长、基于边缘检测、基于聚类分析的算法等。前三种分割算法一般适用于一些简单的图像分割,精度不高,而均值聚类分析算法具有分割结果稳定、算法速度快且能适用于复杂图像的特点。但是在分割可见光遥感图像时,由于材质颜色、灰度、纹理复杂多变,只仅仅依靠算法,很难准确得到分割结果。经过综合考虑,本文采用基于K均值聚类分析算法结合人工的方法分割地物材质。基于K均值聚类分析算法分割地物材质的具体实现步骤如下:
(1)选取聚类中心
利用可见光遥感图像的灰度直方图,分别求出每个波峰对应的灰度值,并将它们按从大到小的顺序进行排列,选取前K个灰度级作为样本的初始值。K个初始聚类中心值为{C1(l),C2(l),…,Ci(l),…,CK(l)},其中,Ci(l)是经l次迭代后的第i个聚类中心值。
(2)特征聚类
通过比较像素灰度与K个初始聚类中心的相似程度,对图像的所有像素进行分类。根据归类公式||X(λ)-Cj(l)||<||X(λ)-Ci(l)||,将像素划分到各个聚类中,归类条件为X(λ)∈Sj(l)。
其中i,j=1,2,…,K,i≠j,Sj(l)代表第l次迭代时聚类j的样本,λ为向量特征维数,X(λ)为样本向量,||X(p)-Cj(l)||为每个样本向量与经l次迭代后的第j个聚类中心值的距离,||X(λ)-Ci(l)||为每个样本向量与经l次迭代后的第i个聚类中心值的距离,当||X(p)-Cj(l)||小于||X(λ)-Ci(l)||的欧氏距离时,表示每个样本向量与经l次迭代后第j个聚类中心值的距离是每个样本向量与这K个聚类中心值的距离的最小值,第l次迭代后,X(λ)样本向量就归类到Sj(l)中。
(3)更新聚类中心值
用步骤(2)得到的新聚类数据,更新K个聚类中心。当新聚类到其类别中的各向量的距离加权和Jj最小时,生成新聚类中心样本Cj(l+1),其中,
Figure GDA0004230063140000081
其中j=1,2,…,K,Nj是归类为Sj的样本向量的数量。
(4)判定收敛条件
当新聚类中心满足|Cj(l+1)-Cj(l)|≤δ时,认为聚类收敛,否则,返回步骤(2)继续迭代。其中,δ是聚类中心变化判断系数。
通过K均值聚类分析的算法步骤得到分割结果后,再用人工的方法对分割结果进行完善,完成遥感图像材质区域分割。图5给出了利用该方法得出的分割结果图,其中,图(a)为可见光遥感图像,图(b)是将图(a)分割为道路和森林两种材质,其中蓝色区域为道路,绿色为森林。可以看出,利用人机交互方法,分割的材质边界自然,且更符合客观事实。
根据分割图像建立热特征预测模型:
所有高于绝对零度的物体都能发出红外辐射,红外图像表现各种材质表面的辐射亮度。在特定的波长范围内,材质表面的辐射亮度分布主要与其表面温度有关,所以,获得材质的温度场分布能够为红外仿真提供关键数据。地物材质的温度与自身的材料、结构特征和导热性能等有关,也有大气与环境的参数有关。因此要获得不同材质的温度分布,需要对不同的地物材质建立相应的热特征预测模型。
大规模可见光遥感图像中,地物材质的长度和宽度,往往比其厚度大几个数量级。为简化计算,将这种大范围的地物材质看作是长度无限长、宽度无限宽,但厚度有限的一维导热材质,热量只在厚度方向进行传导。材质在厚度方向上满足一维导热微分方程:
Figure GDA0004230063140000091
其中,T为绝对温度;t为时间;λg为介质导热率;C为介质热容量;ρm为材料密度;Z为厚度方向的坐标。下边界条件Th可通过实测与经验获得,上边界条件T0可根据地物的温度分布特点,建立热平衡方程求得。上边界主要进行的热量交换有物体的自身辐射、太阳能辐射、大气辐射、地表与大气的显热交换、潜热交换,以及向下传递的能量。满足上边界条件的热平衡方程为:αsEsunlEsky-Mg-Hg-LEg-Gh=0,其中αsEsun项是吸收的太阳短波辐射能量,εlEsky项是吸收的大气长波辐射能量,Mg是地表热辐射能量,Hg是地表与大气的显热交换通量,LEg是地表与大气的潜热交换通量,Gh是地表对下层的热传导项,αs是地表短波吸收率,εl是地表长波发射率。
其中,地表热辐射Mg,通过斯蒂芬-玻尔兹曼公式求得:Mg=εgσT0 4,σ为斯蒂芬-玻尔兹曼常数;εg为物体表面全波段发射率,T0为上边界条件的温度。
其中,显热交换通量Hg,符合显热交换公式:
Figure GDA0004230063140000092
ρa为近地面空气密度,CP为空气定压比热,ra为空气动力学阻力;Ta为参考高度处大气温度。
其中,地表与大气的潜热交换通量LEg,符合潜热交换通量公式:
Figure GDA0004230063140000093
γ为干湿表常数;es(T0)为上边界条件温度为T0时的饱和水汽压;ea为近地面空气水汽压。
其中,地表对下层的热传导项Gh,符合地表的热传导公式:
Figure GDA0004230063140000094
综上,联立一维导热方程与热平衡方程,再通过设置大气与环境参数,即可求得不同地物材质的温度近似解。
将所述热特征预测模型的宏观温度与材质区域分布图进行一一映射,得到基本红外纹理:
在仿真所需要的时间内,同一材质的温度均值可以看作是不变的。将计算得到的各材质的宏观温度与材质区域分布图进行一一映射,得到各材质的宏观温度分布图。图6为地理位置经度120.5、纬度22.5,时间2011年5月23日上午10时某地的宏观温度分布图,进行0~255灰度量化。
对地物材质而言,白天的整体温度要比夜晚高,道路的温度日夜差别较大,植被的温度日夜差别不大。在上午10时,道路温度比植被温度高。从图6中可知,设置各材质温度后,同一材质无纹理细节,显然这样的结果不满足真实性的要求,必须增强温度场的纹理细节。通过在宏观均值温度上添加温度扰动,可以使物体表面辐射分布发生改变,进而实现红外纹理的细节的增强。
对所述宏观温度分布图进行分类调制得到调制图像:
目前,基于可见光遥感图像生成红外纹理温度场,一般采用灰度信息或纹理信息进行调制,调制方法均采用一种模型,调制模式单一。本方案在现有的调制方法基础上,首先根据不同材质的具体特点,对材质进行分类。然后选择合适的调制模型进行分类调制,旨在提高红外纹理温度场分布的真实性和可靠性。基于可见光遥感图像纹理的分类调制框架是以增质预处理后的遥感图像为红外纹理调制模板,具体框架如图7所示:
在本发明的一个实施例中,对所述宏观温度分布图进行第三调制处理得到调制图像,包括,
对所述基本红外纹理进行平坦表面调制,得到初级分类调制模型;平坦表面对应的材质例如道路、水面,其表面没有明显的高低起伏,能够被太阳光均匀的照射,该材质表面的反射率变化是造成材质表面颜色差异的主要原因。从可见光角度来看,白色反射率高,黑色反射率低,材质的反射率与灰度值成正比。从热力学的角度来看,材质的反射率越高,反射太阳光能量越多,吸收的太阳光能量越低,温度越低,材质的反射率与温度呈负相关。综上,平坦表面材质的温度与灰度值呈负相关。平坦表面材质的具体调制方法参考了现有的红外纹理温度场调制方法。
根据平坦表面材质灰度值与温度的关系,按以下步骤进行平坦表面的温度场调制:
(1)计算平坦表面材质的可见光吸收率αν
由于平坦表面材质的灰度值与其温度呈负相关,与其反射率呈正相关。根据材质的灰度值计算得到材质的反射率,即可求得可见光吸收率αν,步骤如下:
根据可见光遥感图像能量传输与传感器成像模型之间的关系,各波段传感器的像素灰度值Gn可表示成:
Figure GDA0004230063140000101
其中n=r,g,b,r,g,b分别为RGB图像的红通道、绿通道和蓝通道,an是传感器增益,bn是传感器偏置,ρn是相应波段的反射率,Esun_n是相应波段的太阳辐射,ar≈ag≈ab=a,br≈bg≈bb=0。
当反射率ρn=1时,灰度Gn=Gmax=255,得出平坦表面红绿蓝波段的反射率ρn为:
Figure GDA0004230063140000111
再根据平坦表面的反射能量与反射率成正比的性质,平坦表面的可见光反射率ρv为:
Figure GDA0004230063140000112
其中,系数cr,cg,cb分别是红、绿、蓝波段太阳辐射占可见光波段太阳辐射的比例,Esun_r为红波段对应的太阳辐射照度,Esun_g为绿波段对应的太阳辐射照度,Esun_b为红波段对应的太阳辐射照度,Esun是太阳辐射总照度。
综上,再根据基尔霍夫定律可得,可见光吸收率αν为:
Figure GDA0004230063140000113
(2)建立平坦表面调制模型
地物材质在接收太阳直接辐射时,地物的短波吸收率αs明显影响温度的高低,αs与温度的关系可以近似看作是线性关系。平坦表面能均匀接收太阳直接辐射,其短波吸收率αs与温度的关系如下:
Figure GDA0004230063140000114
其中,
Figure GDA0004230063140000115
是在t时刻平均短波吸收率为/>
Figure GDA0004230063140000116
的平坦表面温度预测值;T(αs,t)是在t时刻短波吸收率为αs的平坦表面温度调制值;k1(t)是t时刻平坦表面温度随短波吸收率变化的梯度参数;αs为平坦表面短波吸收率;αν为平坦表面吸收率;/>
Figure GDA0004230063140000117
为平坦表面平均吸收率。
材质的短波吸收率αs=k2αν=k2(1-ρν),其中,k2是将平坦表面可见光波段吸收率延拓为短波波段材质吸收率的比例因子。
综上分析可得:平坦表面的调制模型公式为:
Figure GDA0004230063140000118
(3)简化计算的复杂度后的平坦表面调制模型
可见光在0.4~0.76μm范围内太阳辐射能量占太阳辐射总能量的46%,因此可以通过可见光波段地物材质对太阳辐射的能量来近似地物材质对太阳辐射的总能量。为简化计算的复杂度,将t时刻平均短波吸收率为
Figure GDA0004230063140000119
的平坦表面温度预测值近似看作t时刻平坦表面的宏观温度Tg(t),将在t时刻短波吸收率为αs的平坦表面温度调制值T(αs,t)近似为t时刻平坦表面的调制温度T1(t),得到的平坦表面调制模型公式:/>
Figure GDA0004230063140000121
平坦表面温度调制的具体实现为:首先,读取可见光遥感图像中平坦表面的各像素点的RGB通道值,将其代入(3-17)式,计算出各像素点的可见光的吸收率αν;然后,根据分割结果,计算出各材质的平均可见光吸收率
Figure GDA0004230063140000122
最后,利用(3-34)式,将得到的可见光吸收率αν、平均可见光吸收率/>
Figure GDA0004230063140000123
和通过热特征预测得到的宏观温度代入简化的调制模型中,即可对平坦表面的各材质温度进行细节调制。
对所述初级分类调制模型进行起伏表面调制,得到分类调制模型:
起伏表面对应的材质例如森林、草地、植被等,表面有一定高低起伏或孔洞,会导致同一材质存在一定的遮挡问题,遮挡部分不能接收到太阳辐射,未遮挡部分能直接接收到太阳辐射。所以,在小尺度范围内,被遮挡的区域表现为暗区,未被遮挡的区域为亮区。当材质在小尺度范围内的暗区较少时,该区域绝大部分表面能够被太阳照射;而暗区较多时,该区域内绝大部分表面不能接受到太阳光的照射。这种与照射面积呈正相关的太阳照射能量差异,即起伏表面材质在遥感图像中表现出的颜色差异的重要原因。从可见光角度来看,该材质在小尺度范围内照射面积越大,接收到太阳光辐射越多,反射的可见光辐射越多,对应区域的遥感图像灰度值越高。再从热物理角度来看,照射面积越大的材质,吸收更多的短波辐射,表面温度也越高,在小尺度范围内,照射面积越大则高温面积越大,平均温度也越高。所以,起伏表面材质的温度与灰度值呈正相关。
根据起伏表面对应材质温度与灰度值呈正相关的性质,按以下步骤对起伏表面进行温度场调制:
(1)建立直射比例hν与灰度之间的关系
Figure GDA0004230063140000124
表示起伏表面的平均直射比例,hν为某一像素点处的实际直射比例,Δhν为hν偏离/>
Figure GDA0004230063140000125
的量,称为直射比例的扰动量。则灰度值与直射比例之间的关系可表示为:
Figure GDA0004230063140000126
其中,灰度值Gc为分割后的遥感图像中的起伏表面区域进行灰度化得到起伏表面地物类型每一个像素点灰度值;hν为每个像素点接收太阳辐射能量与太阳完全直射像素点的太阳辐射能量的直射比例;Gtotal_in为太阳完全直射像素点时hν为1时像素点的灰度值,它是所调制的起伏表面的最大灰度值;Gtotal_out为像素点全部处于阴影时hν为0时像素点的灰度值,它为该起伏表面的最小灰度值;式中/>
Figure GDA0004230063140000127
是将起伏表面的像素点的值压缩到0~1,式中/>
Figure GDA0004230063140000128
的范围将在[-1/2,1/2]内,则起伏表面地物类型的直射比例的扰动值在[-kg/2,kg/2]之间,可以通过调节kg来改变直射比例的扰动细节,kg为调制比例系数;
(2)建立起伏表面调制模型
根据直射比例与温度呈正相关的关系,它们的公式可以表示为:
Figure GDA0004230063140000131
其中,Tc(t)为t时刻起伏表面的宏观温度,T2(t)为t时刻起伏表面的调制后的温度,k3(t)为t时刻材质温度随直射比例变化的梯度参数;/>
Figure GDA0004230063140000132
表示起伏表面的平均直射比例,当起伏表面材质的直射比例为/>
Figure GDA0004230063140000133
时,调制后的温度仍旧等于宏观温度;
代入式可得,起伏表面温度调制模型公式为:
Figure GDA0004230063140000134
起伏表面温度调制的具体实现为:首先,对可见光遥感图像中的起伏表面进行灰度和模糊处理,得到各起伏表面区域的灰度图像,得到每一点的灰度值Gc,找到材质区域中的最大灰度值Gtotal_in和最小灰度值Gtotal_out;然后通过热特征预测模型算出的各起伏表面的宏观温度代入起伏表面调制模型(3-23)中,即可对起伏表面的各材质温度进行细节调制。
对所述分类调制模型进行分析得到调制图像:
将调制的平坦表面各种材质的温度场与起伏表面各种材质的温度场,逐一映射到各材质区域各像素点上,生成一张包含平坦表面与起伏表面的材质温度场,即红外纹理温度场。将红外纹理温度场映射到0~255灰度级的灰度图像进行显示,各像素点的温度与灰度之间应符合:
Figure GDA0004230063140000135
其中,T3(t)x为t时刻第x个像素点对应的调制后的温度值,T3(t)max、T3(t)min分别为t时刻整个温度场的最大、最小温度值,Grayx为第x个像素点对应的灰度值;
经过上述温度场调制的步骤,即可完成基于可见光遥感图像纹理的红外纹理温度场的分类调制。图8展示了可见光遥感图像的调制结果对比图,其中,(a)为源可见光遥感图像,(a’)增质后的可见光遥感图像,(b)为源遥感图像利用现有的调制方法得到的红外纹理温度场调制结果,(b’)为增质后利用现有的调制方法得到的红外纹理温度场调制结果,(c)为源遥感图像利用该调制方法得到的红外纹理温度场调制结果,(c’)为增质后利用该调制方法得到的红外纹理温度场调制结果。
平坦表面对应材质调制出来的温度场,应与其区域的灰度值成反比,与其吸收率成正比。平坦表面对应材质区域吸收率越高,灰度值越低,温度越高。起伏表面对应材质调制出来的温度场,应与其区域的灰度值成正比,与直射比例成正比。平坦表面对应材质区域直射比例越高,灰度值越高,温度越高。图8中,道路材质表面属于平坦表面,森林材质表面属于起伏表面。
从图8可以看出,(b)(c)(b’)(c’)红外纹理温度场调制结果中道路材质的白色斑马线区域较周围其他区域灰度值偏高,反射的太阳光能量多,吸收的太阳光能量低。从热物理角度来看,吸收的太阳光能量越低,对应的温度值越低。现有的调制方法和本论文的调制方法对道路材质的调制结果基本相同,符合客观物理规律。
但是,两种调制方法调制出的森林材质的温度场有很大的不同,如图(a)红色区域中森林材质的阴影区域,由于其接收太阳辐射的能量较低,该阴影区域应该显示较低的灰度值。(b)(b’)中该阴影区域调制出来的温度比周围的温度还要高,温度场分布显然不合理,其调制模型没有基于森林材质的具体特点,调制模式单一,物理性有待加强。而本文的分类调制方法调制出来的结果图(c)(c’)很好的做到了这一点,红色区域中的阴影区域具有较低的温度值。所以,根据材质的不同特点进行分类调制的方法,调制出来的地物温度场分布更接近于真实,将其作为数据源应用于地物场景的红外仿真中,能够大大提高仿真的逼真度。
从图8中还可以看出,图中的汽车根本不属于道路材质,但是在分割和后续的计算时依然将汽车看作道路材质进行处理,且在实际仿真中,汽车的位置应随着时间发生变化,而图(b)(c)的红外纹理温度场中汽车的位置固定不变。运动的物体会带来纹理固化,利用(a)源遥感图像的温度场调制结果,作为基础数据进行后续的仿真,显然违背了提高仿真灵活性的初衷。运动物体的存在,不仅在一定程度上影响了仿真的时间灵活性,而且对材质的分割也会造成一定的困扰。必须要通过图像增质技术消除了运动物体带来的影响,获得高质量的遥感图像(a’)。利用(a’)作为红外纹理的调制模板,调制生成的红外纹理温度场(b’)(c’),消除了运动物体带来的纹理固化影响,提高了红外纹理温度场分布的真实感,提升了后续仿真的灵活性。
另外,该步骤的分类调制方法是建立在材质的海拔变化不大的基础上的,但是在大规模可见光遥感图像中,地貌特征复杂、环境因素复杂,影响温度场分布的环境因素更多。对于大规模复杂地貌的摇杆图像,需要在此基础上进行进一步调制。
在本发明的一个实施例中,对所述调制图像进行第二调制处理得到红外纹理温度场图像,包括:
对调制图像进行海拔高度调制,得到海拔高度调制模型:
利用地物海拔调制红外纹理温度场时,一般认为场景的经纬度、日期、大气等环境参数是相同的。但是由于地物海拔高度的变化,会导致场景中的大气温度、辐照度及风速发生变化,进而导致相同材质、不同高度的温度场也不同;通过对各种材质温度与热预测模型的环境参数敏感度进行模拟与分析得出,对各种地物材质来说,对环境最为敏感的是大气温度,依次是太阳辐照度,再次风速。所以必须考虑海拔高度对材质温度的影响。
海拔高度低于11km以下的地物材质,其周围环境的大气温度满足:Ta=T0+L0·Z,其中,Ta表示在海拔为Z km的大气气温;L0为对流层大气温度变化的梯度参数,一般为-6.5K/km,T0表示在海拔为0km处的大气气温,Z为海拔高度。
本实施例所采用的遥感图像中的地物海拔高度没有超过1km,根据模拟分析可知,对于海拔低于1km以下的地物材质,虽然遥感图像中的地物材质不同,环境参数不同,但其材质表面温度与海拔高度可近似看作是线性关系,公式如下:T(Z,t)=T(Z0,t)+k4(t)·(Z-Z0),其中,T(Z0,t)为t时刻参考高度Z0处的材质温度预测值,T(Z,t)为t时刻参考高度Z处的材质温度调制值,k4(t)为t时刻材质温度随海拔高度变化的梯度参数。
本步利用分类调制方法生成的温度场T3(t),适用于海拔高度变化不太大的地物材质。如果加上高度的权重因子,调制出来的温度场即能适用于海拔高度范围变化明显的地物材质了。根据此分析,海拔高度调制红外纹理温度场的调制公式为:T(Z,t)=T3(t)+wh·(Z-Z0)。
在实际的温度场调制过程中,首先利用热特征预测模型,计算出某材质在参考高度Z0处与另一特定点的Z1处的温度预测值,代入计算出权重因子wh;其次通过与可见光遥感图像对应的高程数据,获取可见光遥感图像各像素点的高度信息。如果高程数据(DEM)不够精细,不足以得到每个像素点的高度信息,则需要利用内插法得到图像中所有像素点的高度信息;最后将得到的各像素点高度值Z,利用分类调制方法得到的温度场T3(t),权重因子wh及参考高度Z0,代入T(Z,t)=T(Z0,t)+k4(t)·(Z-Z0)即可计算出随海拔高度变化的温度场分布T(Z,t)。
图9给出了红外纹理温度场的调制结果对比图,其中,(a)为源可见光遥感图像,(b)为(a)利用分类调制方法得到的红外纹理温度场调制结果,(c)为(b)进一步加入高度信息后的温度场调制结果;(a’)为(a)增质后的高质量遥感图像,(b’)为(a’)利用分类调制方法得到的温度场调制结果,(c’)为(b’)进一步加入高度信息后的温度场调制结果。图中可见光遥感图像对应的高程数据如图10所示,(a)为数据彩色图,(b)为高程数据灰度图,二者数据一致。
从图9可以看出,首先利用分类调制方法(b)(b’),再进一步利用海拔高度调制出的温度场(c)(c’),能够表现处在不同高度同一材质不同的温度特征,即得到了更高级别真实感的红外纹理温度场。还可以看出,(a)中源遥感图像中存在大片的阴影区域,阴影会导致纹理固化,进而影响调制结果。通过增质预处理、分类调制方法,再利用海拔高度调制后的温度场,不仅解决了纹理固化问题,而且温度场分布更加符合客观物理规律。
对所述海拔高度调制模型进行阳光照射比例调制得到红外纹理温度场:
由于大规模场景中复杂的地物类型结构,高低起伏的地物海拔高度,所处地理位置的不同,不同时刻太阳照射角度的不同,导致材质表面接收的阳光照射比例不同,进而导致材质表面温度不同。此时需要利用太阳照射比对红外纹理温度场进行进一步修正。
在太阳照射角度发生变化时,阴影位置也发生了改变。其中,背阴面的温度低于向阳面温度。
但是,未增质的可见光遥感图像,由于这些固化阴影的存在,若直接进行红外纹理温度场的细节调制,会导致调制结果中,背阴面的温度恒小于向阳面的温度。仿真时间是持续的,固化的阴影显然是我们不需要的。所以,为了提高温度场应用的灵活性,在进行温度场调制生成时,我们必须首先进行阴影去除,以增质预处理后的遥感图像为模板,来进行温度场的调制生成。
但是,利用该调制模板进行温度场的调制生成时,却无法表征出因太阳朝向产生的温度差异。当太阳垂直照射山体,此时,遥感图像中不存在向阳面和背阴面。另,若遥感图像拍摄时恰好在阴天时,由于不同的太阳照射角度对温度影响可以忽略。
在现有的红外纹理温度场调制过程中,未利用增质预处理,只考虑了海拔高度的变化对温度的影响,会导致调制出来的温度场数据固化,不便应用。利用一定的方式去除了阴影,考虑了海拔高度变化,但是未考虑阳光照射比例信息对温度的影响,会导致调制出来的同一材质、同一高度的温度场基本相同,明显与实际不符。
红外场景仿真中,一天内由于太阳照射角度的不同,向阳面和背阴面的位置不断发生变化,向阳面的温度会持续上升,背阴面的温度也会逐渐降低。前步考虑了图像增质和高度信息,未考虑朝向调制出来的温度场,再结合仿真时太阳的具体朝向,即能加工出这次仿真下确定了具体太阳朝向的红外纹理温度场分布。太阳朝向直接影响阳光照射比例,因此大规模的红外纹理温度场生成的关键在于确定阳光的照射比例。
通过上述分析,找到一种能够确定出仿真时刻时,阳光的照射比例的物理量,然后建立该物理量与温度之间的联系即可。山体阴影即使能够反映遥感图像中阳光照射比例信息的物理量。不少工程软件如ENVI,ArcGis等可以进行山体阴影的可视化,山体阴影能够根据一定的太阳朝向,与DEM数据结合,生成可视化的山体阴影效果图。
在本发明的一个实施例中,对所述海拔高度调制模型进行阳光照射比例调制得到红外纹理温度场,包括:
对海拔高度调制模型进行山体阴影计算,得到初级阳光照射比调制模型:
山体阴影与各点的海拔高度Z有关,与太阳的高度角和方位角有关,也与像元所处的坡度和坡向有关。它能够具象的表征由于太阳朝向造成的阳光照射比例,并且能够随着时间的变化而发生变化。图11中各参数的意义为:θ1为太阳天顶角的弧度数;
Figure GDA0004230063140000161
为像元的法向量;θ为太阳高度角的弧度数;θ2为太阳方位角的弧度数;θ3为像元的坡度弧度数;θ4为像元的坡向弧度数;/>
Figure GDA0004230063140000162
为像元法向量/>
Figure GDA0004230063140000163
在地平面上的投影向量;/>
Figure GDA0004230063140000164
为太阳光线在地平面上的投影向量。
山体阴影的计算公式为:Hs=255*(sinθ1*sinθ3*cos(θ24)+cosθ1*cosθ3),其中,上式已将山体阴影Hs量化到0~255范围内,在实际计算时,若Hs<0时,可令Hs=0。为求解上式,必须得到计算公式内的每个角度。
(1)天顶角θ1
太阳天顶角与太阳高度角互为余角,太阳高度角θ的计算公式为:θ=arcsin(sinσ1×sinδ1+cosσ1×cosδ1×cosτ),其中Δ为地理纬度,δ为太阳赤纬角,τ为时角,则天顶角θ1
Figure GDA0004230063140000165
(2)方位角θ2
太阳方位角θ2的计算公式为:
Figure GDA0004230063140000166
(3)像元坡度θ3
在一个以e点像元为中心的3×3邻域内,设每个像素点对应的高度值为a~i,如图12所示。
则e点在x方向的变化率可表示为
Figure GDA0004230063140000171
e点在y方向的变化率为/>
Figure GDA0004230063140000172
其中,z为某像素点的高度值,m为DEM的栅格大小,可通过读取DEM数据获得。则e点像元的坡度θ3为/>
Figure GDA0004230063140000173
其中,z_f为协调z方向的单位与xy平面上的单位的系数,默认为1。
(4)像元坡向θ4
Figure GDA0004230063140000174
时,/>
Figure GDA0004230063140000175
其中atan2(·)是反正切函数,/>
Figure GDA0004230063140000176
或/>
Figure GDA0004230063140000177
为0时,程序不会报错。当根据(3-35)式求解得到的θ4<0,则令θ4=2π+θ4
Figure GDA0004230063140000178
时,/>
Figure GDA0004230063140000179
在实际的温度场调制过程中,首先,通过与可见光遥感图像对应的高程数据,获取可见光遥感图像各像素点的高度信息和经纬度信息。一般处理的可见光遥感图像,经纬度变化不太大,读取第一个像素点的经纬度信息即可。其次,通过确定的仿真时间和获取的高度信息与经纬度信息,计算出太阳天顶角、方位角、坡度与坡向。最后,将四个角度代入到式(4-28)中,即可得出遥感图像每个像素点的山体阴影值Hs。
对初级阳光照射比调制模型进行阳光照射比例调制,得到红外纹理温度场:
图13为可见光遥感图像的山体阴影值,其中,(a)为2011年5月23日上午10时拍摄的可见光遥感图像,(b)为该图像每一个像素点的山体阴影值Hs(灰度值已量化到0~255)。通过与可见光遥感图像相比,可以看出,山体阴影图表征阳光照射比例信息,山体阴影值与阳光照射比例呈正相关。
还可以看出,同一材质(图中的森林材质)中阳光照射比例越低的区域,即山体阴影值越低灰度值越暗的区域,材质表面接收到的太阳辐射越少,甚至只剩下环境光的辐射,对应在遥感图像上的像素点的灰度值越低。从热物理的角度来看,阳光照射比例越低的区域,吸收的辐射能量越少,表面的温度越低。根据山体阴影与温度的正相关关系,建立山体阴影与温度之间的模型关系。
材质表面的温度是由于来自太阳的辐射和环境光的辐射,且主要来源于太阳光的辐射,处于阳光照射比例极低的材质区域由于甚少接收到太阳光照射,材质区域表面的温度基本来源于周围环境光的辐射;另一方面,阴天的时候,由于云彩的遮挡等原因,材质表面接收到的太阳光的辐射很少,材质表面的温度也基本是由于周围环境光的辐射造成的。
通过上述分析,晴天时,阳光照射比例极低的材质区域的温度是基本由于环境光的辐射造成的,阴天时材质表面的温度也基本是由于环境光造成的。阳光照射比例越低的材质区域,其温度值越低,所以,将同一材质温度最低值Ts近似等于阴天时材质表面的平均温度Ty
同一材质的温度与山体阴影Hs呈正相关,将材质表面温度与山体阴影近似看作是线性关系,可得出:
Figure GDA0004230063140000181
Ts(t)为t时刻各材质的最低温度,Ty(t)为t时刻阴天时对应的各材质的平均温度,Ths(t)为t时刻参考山体阴影处材质的温度预测值,hs(t)为t时刻各像素点的Hs的值,hs0(t)为t时刻各材质的Hs最小值,kh(t)为t时刻材质温度随Hs变化的梯度参数。
其中,Ty(t)可根据热预测模型联立的热平衡方程和热传导方程,通过设置和计算获取阴天时的气温、太阳辐射、大气辐射、风速、相对湿度等环境参数,采用有限差分算法,进而计算出阴天时一天的各时刻不同材质的温度近似解。
通过海拔高度调制红外纹理温度场,将海拔高度信息对温度的影响,加入到红外纹理温度场调制方法中,改进了红外纹理温度场调制模型,得到了材质表面温度随高度变化的温度场T(Z,t)。经过对山体阴影Hs与温度之间的关系分析得出,山体阴影Hs与同一材质的温度呈正相关,且山体阴影表征阳光照射比例。将阳光照射比例信息加入到红外纹理温度场调制方法中,再次改进红外纹理温度场调制模型,可使生成的红外纹理温度场分布更加贴合实际。
改进后的红外纹理温度场调制模型如下:
Figure GDA0004230063140000182
其中,whs(t)为增加的山体阴影变化的权重因子,T(Z,t)为t时刻经过增加了高度信息得到的温度场分布,T(hs,t)为t时刻加入了阳光照射比例信息后的温度场分布。其中whs(t)和kh(t)为参数因子。
实际计算时,首先利用实验室软件设置阴天时的大气参数得到随时间变化的温度值Ty(t),然后通过仿真时间、高度角、方位角、经纬度信息,计算出各像素点随时间变化的山体阴影值hs(t),并获得各材质的最低山体阴影值hs0(t)。最后利用海拔高度调制的调制模型公式,计算出仿真时刻的温度场分布T(Z,t),并将所得的T(Z,t)、Ty(t)、hs(t)、hs0(t)和参数因子whs(t)、kh(t)代入
Figure GDA0004230063140000191
计算出既加入高度信息,又加入了阳光照射比例信息的红外纹理温度场分布T(hs,t)。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明,不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干简单推演或替换,都应当视为属于本发明的保护范围。

Claims (3)

1.基于可见光遥感图像的温度场调制方法,其特征在于,步骤如下:
根据所述可见光遥感图像得到阴影区域;
对所述阴影区域进行补偿,得到补偿图像;
去除所述补偿图像中的运动物体,得到红外纹理调制模板;
对所述红外纹理调制模板进行第一调制处理得到调制图像;
对所述调制图像进行第二调制处理得到红外纹理温度场图像;
所述对所述阴影区域进行补偿,得到补偿图像,包括:
(1)RGB通道处理:将RGB可见光遥感图像分为三个通道进行一一处理,将被作为Minkowski范式输入的各通道的灰度值;
(2)分割阴影和非阴影区域:利用阴影检测区域,将各通道的灰度图像分为阴影区域与非阴影区域;
(3)计算光源颜色:将各通道的阴影区域与非阴影区域利用Minkowski范式进行计算,得出各通道的阴影区域与非阴影区域的光源颜色;
(4)全部恢复到标准光照条件下:阴影区域恢复到标准光照条件下的表达式为:
Figure FDA0004230063130000011
其中,yb(i,j)为某通道阴影区域恢复到标准光照条件下的灰度值,
y(i,j)为该通道阴影区域的灰度值,e1表示阴影区域的光照颜色;
非阴影区域恢复到标准光照条件下的表达式可表示为:
Figure FDA0004230063130000012
其中,gb(i,j)为某通道非阴影区域恢复到标准光照条件下的灰度值,g(i,j)为该通道非阴影区域的灰度值,e2表示非阴影区域的光照颜色;
在标准光照条件下,yb(i,j)和gb(i,j)的场景的光照条件一样,所以两者的光源颜色也应该相等;
(5)恢复到非阴影光照条件下:标准光照条件恢复到非阴影区域的表达式可表示为:g(i,j)=gb(i,j)*e2
则阴影区域恢复到非阴影区域的灰度值yg(i,j)应满足:
Figure FDA0004230063130000013
(6)得到补偿图像:将各通道的阴影区域恢复到非阴影区域的灰度值yg(i,j)保持非阴影区域g(i,j)值不变,赋给各通道对应的区域,即可得到各通道的灰度值,然后,将各通道按原来的R,G,B通道输出,即可得到补偿图像;
其中,阴影区域和非阴影区域满足Minkowski范式:
Figure FDA0004230063130000021
其中,e为光源颜色,f为图像三色通道的灰度值,k为比例系数,χ∈[1,∞)且为整数,属于权重因子,M为图像长度,N为图像宽度;
对所述调制图像进行第二调制处理得到红外纹理温度场图像,包括:
对调制图像进行海拔高度调制,得到海拔高度调制模型;
对海拔高度调制模型进行山体阴影计算,得到初级阳光照射比调制模型;
对初级阳光照射比调制模型进行阳光照射比例调制,得到红外纹理温度场。
2.根据权利要求1所述的基于可见光遥感图像的温度场调制方法,其特征在于:对所述红外纹理调制模板进行第一调制处理得到调制图像,包括,
对红外纹理调制模板进行分割处理与材质识别,得到分割图像;
根据分割图像建立热特征预测模型;
将所述热特征预测模型的宏观温度与材质区域分布图进行一一映射,得到基本红外纹理;
对所述宏观温度分布图进行分类调制得到调制图像。
3.根据权利要求2所述的基于可见光遥感图像的温度场调制方法,其特征在于:对所述宏观温度分布图进行分类调制得到调制图像,包括,
对所述基本红外纹理进行平坦表面调制,得到初级分类调制模型;
对所述初级分类调制模型进行起伏表面调制,得到分类调制模型;
对所述分类调制模型进行分析得到调制图像。
CN201811333483.1A 2018-11-09 2018-11-09 基于可见光遥感图像的温度场调制方法 Active CN109671111B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811333483.1A CN109671111B (zh) 2018-11-09 2018-11-09 基于可见光遥感图像的温度场调制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811333483.1A CN109671111B (zh) 2018-11-09 2018-11-09 基于可见光遥感图像的温度场调制方法

Publications (2)

Publication Number Publication Date
CN109671111A CN109671111A (zh) 2019-04-23
CN109671111B true CN109671111B (zh) 2023-07-07

Family

ID=66142660

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811333483.1A Active CN109671111B (zh) 2018-11-09 2018-11-09 基于可见光遥感图像的温度场调制方法

Country Status (1)

Country Link
CN (1) CN109671111B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111402233B (zh) * 2020-03-16 2021-02-12 清华大学 航空航天飞行器热防护部件测试装置
CN113379759A (zh) * 2021-05-11 2021-09-10 中国资源卫星应用中心 一种光学遥感卫星影像水体自动提取方法
CN115077386B (zh) * 2022-08-19 2022-12-16 南京木木西里科技有限公司 一种水溶胶表面全自动测量装置、系统及其测量方法
CN115859841B (zh) * 2023-02-28 2023-05-05 湖南光华防务科技集团有限公司 一种灭火弹挂飞温度仿真方法和系统

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8439503B2 (en) * 2008-08-06 2013-05-14 Disney Enterprises, Inc. Infrared imaging projection
WO2011123483A1 (en) * 2010-03-31 2011-10-06 Omron Scientific Technologies, Inc Method and apparatus for generating texture in a three-dimensional scene
CN104637073B (zh) * 2014-12-30 2017-09-15 华中科技大学 一种基于太阳照射阴影补偿的带状地下结构探测方法
CN106644092B (zh) * 2017-01-05 2019-01-15 西安电子科技大学 基于可见光信息的红外纹理温度场调制方法
CN107516297A (zh) * 2017-07-14 2017-12-26 西安电子科技大学 基于遥感图像的地物红外纹理调制模板生成方法

Also Published As

Publication number Publication date
CN109671111A (zh) 2019-04-23

Similar Documents

Publication Publication Date Title
CN109671111B (zh) 基于可见光遥感图像的温度场调制方法
Li et al. Quantifying the shade provision of street trees in urban landscape: A case study in Boston, USA, using Google Street View
CN108830844B (zh) 一种基于多时相高分辨率遥感影像的设施蔬菜提取方法
CN110569797B (zh) 地球静止轨道卫星影像山火检测方法、系统及其存储介质
CN106644092B (zh) 基于可见光信息的红外纹理温度场调制方法
Li et al. Using Google street view for street-level urban form analysis, a case study in Cambridge, Massachusetts
CN109671038A (zh) 一种基于伪不变特征点分类分层的相对辐射校正方法
CN115272876A (zh) 一种基于深度学习的遥感图像船舶目标检测方法
Wang et al. Atmospheric correction methods for GF-1 WFV1 data in hazy weather
Cui et al. Combined Model Color-Correction Method Utilizing External Low-Frequency Reference Signals for Large-Scale Optical Satellite Image Mosaics.
Wang et al. A PM2. 5 concentration estimation method based on multi-feature combination of image patches
CN116519557A (zh) 一种气溶胶光学厚度反演方法
Song et al. Cloud detection method based on clear sky background under multiple weather conditions
CN110070513A (zh) 遥感影像的辐射校正方法及系统
Aher et al. Rainfall estimation over roof-top using land-cover classification of google earth images
CN116246272A (zh) 针对国产卫星多光谱图像质量标记的云雪区分方法
He et al. Relative radiometric correction of high-resolution remote sensing images based on feature category
Parmar et al. Land use land cover change detection and its impact on land surface temperature of malana watershed kullu, Himachal Pradesh, India
CN109146896A (zh) 一种低对比度情况下太阳能光伏板的有效分割方法
Saputro et al. Study of land surface temperature using remote sensing satellite imagery in Makassar, South Sulawesi
He et al. Spatial and temporal characteristics of surface albedo in Badain Jaran Desert, China
Zhang et al. Basis image decomposition of outdoor time-lapse videos
Volchok et al. Scene to scene radiometric normalization of the reflected bands of the Landsat Thematic Mapper
CN111680659A (zh) 国际空间站rgb夜间灯光图像的相对辐射归一化方法
Li et al. Intelligent recognition of point source target image control points with simulation datasets

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