CN114398812A - 滤除分布式地下建筑背景热流通量的反演探测方法及装置 - Google Patents
滤除分布式地下建筑背景热流通量的反演探测方法及装置 Download PDFInfo
- Publication number
- CN114398812A CN114398812A CN202111680054.3A CN202111680054A CN114398812A CN 114398812 A CN114398812 A CN 114398812A CN 202111680054 A CN202111680054 A CN 202111680054A CN 114398812 A CN114398812 A CN 114398812A
- Authority
- CN
- China
- Prior art keywords
- mountain
- radiation
- solar
- underground building
- model
- 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
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 51
- 238000001914 filtration Methods 0.000 title claims abstract description 45
- 230000004907 flux Effects 0.000 title claims abstract description 37
- 230000005855 radiation Effects 0.000 claims abstract description 331
- 238000000034 method Methods 0.000 claims abstract description 94
- 238000012546 transfer Methods 0.000 claims abstract description 24
- 230000007797 corrosion Effects 0.000 claims abstract description 15
- 238000005260 corrosion Methods 0.000 claims abstract description 15
- 238000004364 calculation method Methods 0.000 claims description 24
- 239000011435 rock Substances 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 16
- 230000008878 coupling Effects 0.000 claims description 8
- 238000010168 coupling process Methods 0.000 claims description 8
- 238000005859 coupling reaction Methods 0.000 claims description 8
- 230000005484 gravity Effects 0.000 claims description 8
- 238000002834 transmittance Methods 0.000 claims description 8
- 101000641239 Homo sapiens Synaptic vesicular amine transporter Proteins 0.000 claims description 7
- 102100034333 Synaptic vesicular amine transporter Human genes 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 230000003595 spectral effect Effects 0.000 claims description 6
- 230000008030 elimination Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 claims description 3
- 230000003628 erosive effect Effects 0.000 claims 1
- 238000005516 engineering process Methods 0.000 abstract description 5
- 239000010410 layer Substances 0.000 description 27
- 239000002689 soil Substances 0.000 description 18
- 238000010586 diagram Methods 0.000 description 12
- 230000008859 change Effects 0.000 description 8
- 230000011218 segmentation Effects 0.000 description 7
- 238000004088 simulation Methods 0.000 description 6
- 230000005540 biological transmission Effects 0.000 description 5
- 239000002344 surface layer Substances 0.000 description 5
- 238000013178 mathematical model Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000009933 burial Methods 0.000 description 3
- 101100372758 Danio rerio vegfaa gene Proteins 0.000 description 2
- 241000196324 Embryophyta Species 0.000 description 2
- 101150030763 Vegfa gene Proteins 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000008033 biological extinction Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 241000127225 Enceliopsis nudicaulis Species 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 239000003673 groundwater Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 229910052755 nonmetal Inorganic materials 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
- G06T5/30—Erosion or dilatation, e.g. thinning
-
- 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/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Computer Hardware Design (AREA)
- Computer Graphics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Image Processing (AREA)
Abstract
本发明提供了滤除分布式地下建筑地质背景热流通量的反演探测方法及装置，属于地球物理学学科理论与遥感技术交叉的领域，方法包括：从地下建筑对应山体的红外遥感图像中滤除山体表面太阳辐射总能量和山体表面植被热辐射能量；结合DEM数据，在红外图像中分层滤除山体地质背景热流通量，获取分布式地下建筑扰动信号图像；利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；本发明实现了对分布式地下建筑的反演探测定位。
Description
技术领域
本发明属于地球物理学学科理论与遥感技术交叉的领域，更具体地，涉及一种滤除分布式地下建筑地质背景热流通量的反演探测方法及装置。
背景技术
随着城市化的进程不断加快，越来越多的人口涌入城市，城市用地显得越来越紧张，这使得更多的建筑不得不选择修建在地下如：地下停车场，地下仓库，地下餐厅等；有很多建筑是出于安全因素的考虑修建在地下，如：大型地下油库，军事设施等。这些都是分布式地下建筑的典型例子，一旦这些分布式地下建筑出现故障，就面临勘察故障出现位置难度大的问题。地下建筑的探测在民用上有着重要的意义，对分布式地下建筑物的反演探测和方法研究有着重要的意义。
目前，国内外利用遥感技术对地下建筑探测方法的研究采用的手段多为基于多时相红外图像的地下建筑探测方法，该方法数据样本量大且获取难度较高，研究对象目前也仅限于埋深较浅且具有一定特性的小型地下建筑，存在很大的局限性；基于微波红外增强技术的地下建筑探测方法的研究对象为浅埋金属地下建筑，根据微波功率传输数学模型增强红外图像，突出地下建筑/背景信号强度对比，难以适用于对其他非金属类型地下建筑的探测，并且环境状况和背景特性对该方法的探测结果具有一定影响。
针对较深层地下建筑红外遥感探测方法的研究，此前提出一种利用双时相探测定位山体中地下建筑的背景滤波方法，该方法避免了多时相红外图获取困难的问题，但该方法未考虑山体表面植被对山体表面接受太阳辐射的影响，基于图像处理的手段计算阴影区域，同时计算山顶起每一层厚度山体所含的热量来滤除山体背景能量理论依据不够完整。
本发明通过结合高光谱图像仿真计算植被传热模型来滤除山体表面植被影响，通过真实DEM数据获得实际海拔信息，结合经纬度信息计算出每个数据栅格点的太阳光线遮挡角度，从而得到地表阴影区域，然后通过计算山体地表太阳辐射模型去除太阳辐射影响，采用试错法对不同的温度条件和模型参数进行仿真计算，建立山体分布式地下建筑温度场模型，然后通过分层迭代滤除地质背景热流通量，降低探测结果误差，更加符合真实情况，有利于分布式地下建筑的准确反演探测与定位。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种滤除分布式地下建筑地质背景热流通量的反演探测方法及装置，旨在解决现有的分布式地下建筑在红外图像中经地层传导调制后信号变弱，无法有效反演探测定位的问题。
为实现上述目的，一方面，本发明提供了一种滤除分布式地下建筑地质背景热流通量的反演探测方法，包括以下步骤：
(1)根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；且基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
(2)从地下建筑对应山体的红外遥感图像中滤除山体表面太阳辐射总能量和山体表面植被热辐射能量；
(3)结合DEM数据，在经过步骤(2)处理后的红外图像中分层滤除山体地质背景热流通量，获取扰动信号图像；
(4)利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
(5)根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；
(6)对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
进一步优选地，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建地下建筑块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述地下建筑块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
进一步优选地，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影区域；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布；
获取山体表面太阳散射辐射分布的方法，包括以下步骤：
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，当扇区足够小时，该扇区内的太阳散射辐射可认为是相同，将若干区域中的太阳散射辐射量分别计算，叠加获取整个扇区的太阳散射辐射量；
将所有要计算散射辐射的位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布对可见视域范围内的散射辐射进行累加，得到该观测处的总散射辐射；
对山体地区所有位置进行上述操作，获取所述山体表面太阳散射辐射分布。
进一步优选地，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
进一步优选地，地下建筑的位置信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置。
另一方面，本发明提供了一种滤除分布式地下建筑地质背景热流量的反演探测装置，包括：
太阳总辐射能量分布计算模块，用于根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；
地表植被辐射分布计算模块，用于基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；
山体三维温度场计算模块，用于基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
扰动信号图像获取模块，用于从山体的红外图像中滤除山体表面太阳辐射总能量、山体地表植被辐射能量；再结合DEM数据，分层滤除山体地质背景热流通量，获取地下建筑的扰动信号图像；
噪声消除模块，用于利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
地下建筑位置获取模块，用于根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
进一步优选地，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建地下建筑块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述地下建筑块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
进一步优选地，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布。
进一步优选地，获取山体表面太阳散射辐射分布的方法，包括以下步骤：
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，当扇区足够小时，该扇区内的太阳散射辐射可认为是相同，将若干区域中的太阳散射辐射量分别计算，叠加获取整个扇区的太阳散射辐射量；
将所有要计算散射辐射的位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布对可见视域范围内的散射辐射进行累加，得到该观测处的总散射辐射。对山体地区所有位置进行上述操作，获取所述山体表面太阳散射辐射分布。
进一步优选地，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
进一步优选地，地下建筑的位置信息和埋深信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置信息。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
本发明采用分布式地下建筑所在山体与空气层之间热辐射模型，结合DEM数据，计算山体表面太阳辐射总能量和结合高光谱图像计算山体地表植被热辐射模型，山体分层迭代计算地质背景热流通量；滤除山体背景热流量得到地下建筑热传播能量扰动信号；充分考虑外部环境所带给山体的影响，如山体岩层热传播能量、地下建筑热传播扩散能量、太阳辐射能量以及地形因素带给红外数据的扰动差异；以本发明解决了分布式地下建筑在红外图像中经地层传导调制后信号变弱，无法有效反演探测定位等问题，从而实现对分布式地下建筑的反演探测定位。
本发明考虑了山体表面太阳辐射能量，在DEM各个数据栅格上充分考虑太阳直接辐射和散射辐射的影响，相比于多时相红外图像分析地下建筑信号随时间变化规律，建立数学模型计算求解地下建筑位置的方法相比，更加精确的考虑了山体外部环境对地下建筑能量扰动的影响，使得探测反演结果更加精确；相比于对山体红外图像进行统计分析，寻找有地下建筑的山体与旁边下无地下建筑的山体的形成的特征模式，利用对可能含有地下建筑的山体环境的红外图像中进行一定方向的遍历搜寻该特征模式，从而实现分布式地下建筑的探测、定位的方法，对分布式地下建筑反演探测定位的精准度很高。
附图说明
图1是本发明实施例提供的地下建筑反演探测方法流程图；
图2是本发明实施例提供的喻家山地区红外图像示意图；
图3是本发明实施例提供的喻家山地区DEM数据；
图4是本发明实施例提供的喻家山区高光谱图像示意图；
图5是本发明实施例提供的喻家山区域山体三维温度场分布；
图6是本发明实施例提供的地表某点接收的太阳辐射示意图；
图7是本发明实施例提供的山体阴影示意图；
图8是本发明实施例提供的喻家山地区太阳直接辐射图；
图9是本发明实施例提供的视域计算中观测点设置；
图10是本发明实施例提供的观测点对应半球视域(灰色表示可见，黑色表示不可见)；
图11是本发明实施例提供的观测点对应半球视域散射辐射分布图；
图12是本发明实施例提供的喻家山散射辐射分布图；
图13是本发明实施例提供的喻家山地区地表接收到的太阳总辐射能量分布图；
图14是本发明实施例提供的喻家山区域山体地表植被热辐射能量图；
图15(a)是本发明实施例提供的滤除20m山体厚度所包含热流场能量的结果图像；
图15(b)是本发明实施例提供的滤除40m山体厚度所包含热流场能量的结果图像；
图15(c)是本发明实施例提供的滤除60m山体厚度所包含热流场能量的结果图像；
图15(d)是本发明实施例提供的滤除80m山体厚度所包含热流场能量的结果图像；
图16是本发明实施例提供的探测结果图像示意图；
图17是本发明实施例提供的探测结果二值化图像；
图18是本发明实施例提供的闭运算和开运算结果示意图；
图19是本发明实施例提供的连通域分割结果示意图；
图20是本发明实施例提供的探测结果上显示连通域分割结果示意图；
图21是本发明实施例提供的地下建筑所在位置区域示意图；
图22是本发明实施例提供的地下建筑位置结果示意图；
图23是本发明实施例提供的可见光图像上的地下建筑位置示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明总体思路为：针对山体环境中，利用红外遥感图像得到的山体表面真实红外辐射对不同厚度的山体地质背景热流通量依次进行分层滤除，有效的保护地下建筑信息，得到较为准确的地下建筑扰动信号分布，实现对地下建筑的反演定位。
本发明提供了一种分布式地下建筑地质背景热流量的滤除方法，包括以下步骤：
一方面，本发明提供了一种滤除分布式地下建筑地质背景热流通量的反演探测方法，包括以下步骤：
(1)根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；且基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
(2)从地下建筑对应山体的红外遥感图像中滤除山体表面太阳辐射总能量和山体表面植被热辐射能量；
(3)结合DEM数据，在经过步骤(2)处理后的红外图像中分层滤除山体地质背景热流通量，获取扰动信号图像；
(4)利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
(5)根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；
(6)对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
进一步优选地，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
进一步优选地，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影区域；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布；
获取山体表面太阳散射辐射分布的方法，包括以下步骤：
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，当扇区足够小时，该扇区内的太阳散射辐射可认为是相同，将若干区域中的太阳散射辐射量分别计算，叠加获取整个扇区的太阳散射辐射量；
将所有要计算散射辐射的位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布对可见视域范围内的散射辐射进行累加，得到该观测处的总散射辐射。对山体地区所有位置进行上述操作，获取所述山体表面太阳散射辐射分布。
进一步优选地，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
进一步优选地，地下建筑的位置信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置。
另一方面，本发明提供了一种滤除分布式地下建筑地质背景热流量的反演探测装置，包括：
太阳总辐射能量分布计算模块，用于根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；
地表植被辐射分布计算模块，用于基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；
山体三维温度场计算模块，用于基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
扰动信号图像获取模块，用于从山体的红外图像中滤除山体表面太阳辐射总能量、山体地表植被辐射能量；再结合DEM数据，分层滤除山体地质背景热流通量，获取地下建筑的扰动信号图像；
噪声消除模块，用于利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
地下建筑位置获取模块，用于根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
进一步优选地，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述地下建筑块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
进一步优选地，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布。
进一步优选地，获取山体表面太阳散射辐射分布的方法，包括以下步骤：
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，当扇区足够小时，该扇区内的太阳散射辐射可认为是相同，将若干区域中的太阳散射辐射量分别计算，叠加获取整个扇区的太阳散射辐射量；
将所有要计算散射辐射的位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布对可见视域范围内的散射辐射进行累加，得到该观测处的总散射辐射。对山体地区所有位置进行上述操作，获取所述山体表面太阳散射辐射分布。
进一步优选地，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
进一步优选地，地下建筑的位置信息和埋深信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置信息；对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
实施例
如图1所示，以喻家山为例进行说明，本发明所提供的一种滤除分布式地下建筑地质背景热流通量的反演探测方法，包括下述几个步骤：
(1)获取包含地下建筑的可见光/红外图像以及对应区域的DEM数据和高光谱图像；
喻家山地区红外图像如图2所示，喻家山DEM数据图像如图3所示，图4是喻家山高光谱图像；
(2)地下建筑块状热能辐射模型的建立
以地球物理学、传热学、遥感技术等学科的基本理论为指导，分析地下建筑与埋藏环境的物质/能量特征关系及相互作用规律，研究并建立地下建筑区的地下建筑热流场模型，包括以下几个子步骤：
理论研究地下建筑传热特性与地下建筑环境的关系；结合对应喻家山地区DEM数据，利用现有的软件对山体背景进行建模仿真；为了将山体背景传热信息显示出来，利用COMSOL有限元仿真软件对山体背景三维温度场分布进行了仿真计算；
山体地下建筑耦合模型传热过程中的热量守恒可描述如下：
Qd+Qv＝Q
其中，Qv和Qd分别表示对流热通量和传导热通量，Q表示外界输入的热通量，ρ和 Cρ分别表示介质的密度和恒压比热容，w表示流体速度，k表示介质的热导率，h表示海拔梯度，
表示温度在垂直方向上的梯度变化，计算时可采用差分近似，
地下建筑内部温度恒定作为地下建筑表面边界条件，可表示为：
Tin＝f(y)＝T0
通过联立求解所述公式可得到喻家山山体耦合模型的表面温度场分布，本发明结合相关参数利用COMSOL软件进行建模仿真计算，得到最终的喻家山区域山体三维温度场分布情况如图5所示；
(3)计算并滤除地下建筑表面太阳辐射能量
太阳辐射对山体分布式地下建筑的探测有着较大的干扰，同时喻家山地区地形起伏明显，在直射太阳光的影响下会产生山体阴影，从而造成山体表面各部分接收到的太阳辐射不同，从而改变了该地区的热量分布；在地下建筑热流场分析与建模过程中，通过建模仿真计算得到的山体表面温度数据，是没有考虑外界环境影响的理想情况下的理论数据，与利用红外遥感图像获得的山体表面真实红外辐射数据具有一定的差异；在考虑喻家山地区的干扰因素时候，需要对太阳辐射、山体地表植被热辐射能量和山体背景热辐射能量这三大外界干扰因素进行研究计算并去除，在滤除太阳辐射时考虑阴影的影响，以保证探测结果的真实性和有效性；
地表接收到的太阳辐射主要包括直接辐射、散射辐射以及反射辐射，由于反射辐射占比小且仅在特殊条件下才会对地表接收到的太阳辐射产生较大影响，通常太阳辐射计算模型中不考虑反射辐射部分；地表某点接收到的太阳辐射如图6所示。
(3.1)山体阴影的计算方法
本发明所使用的红外图像拍摄于10月08日的下午14点47分，拍摄时天气晴朗无云，此地表存在着山体阴影；阴影区域产生的原因在于太阳直接辐射无法照射到地表，因此阴影区域只有太阳散射辐射；
山体阴影的主要影响因素是太阳方位角、太阳高度角和山体的地形；在计算山体阴影的时候，较高的山峰会在其后的区域形成阴影，阴影区域内部分无法接收来自太阳的直接辐射；
山体阴影可由下式计算得出：
H＝cos(Azenith)*cos(Aslope)+sin(Azenith)*sin(Aslope)*cos(Aazimuth-Aaspect)
其中，Aazimuth为太阳入射方向；Azenith为太阳入射角度，数值上等于π/2-Azenith；Aslope为坡度，即山体表面每个像元开始的最陡坡降；Aaspect为坡向，最陡坡降的方向就是坡向的朝向；
其中，坡度和坡向的计算可以用3*3的窗口依次扫过高程图的每一个像素，通过中心差分的方式计算在x和y方向上的高程变化率，进而推出在该像元位置的最陡坡降、坡度和坡向；计算得到的山体阴影如图7所示；
(3.2)太阳直接辐射分布计算
晴天无云的条件下，地表某点的太阳直接辐射Ebh的计算公式如下：
Ebh＝E0×τb×sinh
其中，E0为大气层上界所接收到的太阳辐射；
其中，Esc为太阳常数，世界气象组织(WHO)公布的参考值为1367W/m2，n为计算日与1月1日之间相隔的天数；
τb为太阳直接辐射的大气透过率，具体为：
τb＝0.56×(e-0.56m(h)+e-0.096m(h))×k1
其中，h为太阳高度角，在获取的红外图像信息中已经给出，参数m为大气质量，参数k1为经验值，一般取0.8～0.9；
将山体阴影叠加在太阳直接辐射分布图上，阴影区域的太阳直接辐射置零，非阴影区域的太阳直接辐射不变，得到的喻家山地区太阳直接辐射如图8所示；
(3.3)地表接收的太阳散射辐射计算
参考半球视域的计算方法，通过计算太阳散射辐射的半球视域分布从而获得山体接收的太阳散射辐射；半球视域可以理解成在某个位置向天空仰视范围内未被遮挡的区域；对视域的定义如下：识别DEM数据，通过计算获得DEM上某观测位置所能看到的像元位置；在计算视域的算法中，将输出定义为所在位置的可见性，是否可在观测位置直接看到地下建筑位置，以二值图表示可见性，即将可直接看到的位置设置为1，其他位置设置为0；
使用地球曲率以对视域进行校正，用以进行校正的公式如下：
其中，Zsurface为不考虑地球曲面的视域；Zactual为校正后的视域；Dist为观测点与地下建筑点之间的平面距离，即观测点与地下建筑点的连线在XY平面(地平面)上投影的距离；Diam为地球直径，为12740km；Rrefr为光的折射系数，这是反应大气质量的参数，鉴于晴朗无云能见度佳的拍摄情况，取折射系数为0.13；
以喻家山地区DEM图像中心点为例(海拔高度145m，位于山脊上)，计算DEM图像中心点的半球视域，如图9所示；
在计算太阳散射辐射的时候，图10所示，首先针对观测位置建立起该位置的天空仰视半球视域，并构建这一半球视域内的散射辐射分布；同时将半球域划分为多个扇区，当扇区划分地足够小时，该扇区内的太阳散射辐射可以认为相同；本发明选择5°为扇区分割的角度，并对每个扇区依半径分为100个小区域；将获得的散射辐射所要计算的特定位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布，如图11所示；
对视域可见范围内的散射辐射进行累加，便可以得到该观测位置处的总散射辐射；对喻家山地区所有位置依次进行上述操作，得到喻家山地区的散射辐射分布图如图12所示；
将地表接收到的太阳直接辐射能量和散射辐射能量求和，可得到喻家山地区地表所接收到的太阳总辐射能量分布如图13所示；
(4)结合高光谱图像仿真计算地表植被辐射模型
对于土壤-植被混合的真实地表，各组分之间交错分布，并与外界大气存在能量传输，其能量平衡过程是相互耦合的；
根据能量能量守恒定律，自然地表的能量热平衡方程式可表述如下：
αEsun+εEsky-G-M-LE-H＝0
其中，α表示地表的短波红外吸收率，Esun表示地表接收到的太阳辐照度，ε表示地表的长波红外发射率，Esky表示地表接收到的大气辐照度，G、M、LE和H分别表示地表的下行热传导通量、地表的辐射出射度、潜热和显热交换通量；
在气象学中，土壤-植被-大气传输模型的广义定义是陆面过程中考虑植被在土壤-植被 -大气系统各界面之间能量、物质传输和交换过程中重要作用的物理-化学-生物联合模型，利用该模型研究土壤-植被-大气系统中大气、植物、地表、土壤和地下水层中的水及其相互作用和相互关系十分有效；
具体内容如下：
根据SVAT模型原理提出了土壤-植被混合地表红外辐射温度场的耦合求解模型，利用该模型计算地表仿真温度场和红外辐射场；参考该模型原理，可作如下假设：
①在较短时间间隔内，地表与环境能量交换和传输引起的温度变化十分微小，可认为地表温度场处于稳态；
②忽略植物短时间内因生理作用引起的内部能量变化；
③为简化计算，仅考虑土壤与植被冠层在垂直方向上的能量交换；
④不考虑土壤与植被冠层之间间接产生的潜热和显热交换；
以下公式中，LE表示潜热交换，H表示显热交换，M表示热辐射，G表示下行导热，Qsky和Qsun分别表示达到植被冠层的大气辐射和太阳辐射，各变量下标中cd表示植被冠层吸收的直接辐射量，cr表示植被冠层吸收的土壤反射辐射量，c表示植被冠层吸收的辐射量或自身热辐射量，g表示土壤吸收的辐射量或自身热辐射量；
针对土壤和植被分别建立热平衡方程如下式：
VegF·Rnc＝Hca+Hcg+LEca+LEcg+Gc
(1-VegF)·Rng＝Hga+Hgc+LEga+LEgc+Gg
其中，VegF为植被覆盖率，Rnc为植被冠层净辐射，Rng为土壤净辐射，H、LE、G 分别表示显热交换、潜热交换以及下行热传导通量，下标a、c和g分别代表空气、植被和地表；
Rnc＝Qcdsun+Qcdsky+Qcrsun+Qcrsky-2ηlMc+ηlεcMg
Rng＝Qgsun+Qgsky+ηlεgMc-Mg
Gg≈kl(Tg-T0)/ξz0
其中，kl表示土壤表层的热导率，Tg表示土壤表层辐射温度，T0表示土壤表层向下一定深度处的温度，ξz0为土壤表层厚度；
在太阳辐射和大气辐射经过植被冠层到达土壤表层的过程中，由Lambert-Beer定律可知，植被叶片的阻挡会使辐射量出现衰减，通常采用消光系数描述衰减程度；对于均匀分布的植被群体，植被冠层对太阳辐射和大气辐射的消光系数可分别取0.4和0.8，太阳辐射和大气辐射的冠层透过率的计算方法为：
ηs＝1-exp(-0.4LAI)
ηl＝1-exp(-0.8LAI)
其中，LAI为叶面积指数；叶面积指数采用的是随机数生成的方法进行计算，导致最后得到的地表植被红外辐射存在极大的误差；
考虑到叶面积指数实际采样的难度和估算精度的问题，采用归一化植被指数(normalized difference vegetation index，NDVI)估算叶面积指数；首先对获取的多光谱图像预处理，然后利用多光谱图像可见光红色和近红外波段反射率按照下式计算归一化植被指数，最后采用非线性模型估算叶面积指数；
其中，ρNIR和ρR分别表示多光谱图像的近红外和红色波段的反射率；
冠层吸收的直接太阳辐射和大气辐射以及冠层吸收的土壤反射的太阳辐射和大气辐射计算如下式：
Qcdsun＝(ηs-αc)Qsun
Qcdsky＝ηlεcQsky
Qcrsun＝(1-ηs)αgQsun
Qcrsky＝(1-εg)(1-ηl)εcQsky
其中，αc和εc分别为植被冠层的吸收率和发射率；αg和εg分别为土壤表层的吸收率和发射率；
冠层自身热辐射和土壤自身热辐射计算如下：
Mc＝εcσTc 4
其中，σ为斯蒂芬-玻尔兹曼常数；Tc为待求解的植被冠层辐射温度，Tg为土壤表层辐射温度；
植被冠层与土壤以及空气之间的潜热交换量和显热交换量计算如下：
其中，ρa为空气密度；cp为空气定压比热容；γ为干湿表常数；kd和Da分别为空气温度为Ta时的饱和比湿和空气比湿差；Ra和Rs分别为植被群体的边界层阻力和气孔阻力；根据相关气象数据和公式可得到上述参数的取值；es为太阳辐照度；ea为空气辐照度；
通过联立求解土壤和植被冠层热平衡方程式组得到混合地表温度场分布，进而利用辐射公式得到混合地表红外辐射场分布如图14所示，由于山脊具有能量聚焦效应，并且山脊周围植被稀疏，因而计算出山脊附近红外辐射较强；
(5)分层滤除地下建筑的地质背景热流通量，实现反演探测定位
由于山体中分布式地下建筑的传热规律符合热传导的数学模型，可根据传热学定理和规律，获得喻家山山体和地下建筑的温度场分布，进而计算得到二者叠加的热流场分布；鉴于山体辐射面与空气层之间热辐射模型为三维结构，因此能量场的分布公式也应表示成三维形式，假设在某时刻t0山体表面(x0,y0,z0)处的辐射能量为 BT(x0,y0,z0,t)，根据能量守恒定律，此辐射能量的值可通过地下建筑的热流场能量、山体背景热流场能量、它们与环境之间交换的辐射能量以及接收到的太阳辐射能量计算得到；
BT(x0,y0,z0,t)＝B(x,y,z,t)+T(x,y,z,t)+δ(x,y,z,t)
+S(x0,y0,z0,t)-Bm(x,y,z,t)
其中，BT(x0,y0,z0,t)为山体表面辐射能量，B(x,y,z,t)为山体背景热流场能量，T(x,y,z,t)为地下建筑热流场能量，δ(x,y,z,t)为山体与空气接触面的辐射能量， S(x0,y0,z0,t)为山体表面接收到的太阳辐射能量，Bm(x,y,z,t)为地下建筑对应位置岩体的热流场能量；
在真实环境中，山体与空气接触面的辐射能量很小，相比山体和地下建筑的热流场能量，基本可以忽略；在本发明中，只考虑山体背景热流场能量、地下建筑热流场能量以及接收到的太阳辐射能量之间的关系；
由上述公式可知，地下建筑产生的热流场能量、山体背景热流场能量以及太阳辐射能量叠加形成了可观测到的山体表面红外辐射能量分布；将上式加以变形，可以得到地下建筑热流场能量的计算公式：
T(x,y,z,t)＝BT(x0,y0,z0,t)-B(x,y,z,t)-S(x0,y0,z0,t)
计算从山顶起每一厚度层山体所含的热流量，结合红外遥感图像对每一层背景热流通量进行依次剥离，当部分凸显的扰动随着剥离厚度的增加而不发生变化时，可认为是地下建筑的扰动信号，得到地下建筑所在的埋深信息；
(5.1)计算地质背景各层热流量
理想情形下，单位面积的热流量计算如下公式所示：
其中，cm为岩层介质的容积热容量，k为岩层介质的热导率，
为垂直方向的温度梯度，可用差分近似计算；温度梯度需要利用温度场计算得到，可以看出二者具有近似的线性关系；
在分层滤除地质背景热流通量过程中，需要结合山体背景热流场能量分布，对不同厚度的山体背景所包含的能量进行计算：
其中，A(z)和ε(z)分别为山体内部某一海拔高度处的横截面积和热流的衰减系数； z1(x,y)和z2(x,y)分别为山体表面高度坐标和从表面向下某一深度的高度坐标，通过调整 z2(x,y)实现计算不同厚度的山体背景包含的热流场能量；
(5.2)分层滤除地下建筑的地质背景热流量
伴随着山体内部的热量在岩层中的传导，会形成有方向的热流，类似于X射线具有的穿透性，热流可以从山体底部“穿透”岩层达到地下建筑表面，从而形成矢量热流场；热流在不同岩层中传递时由于地质情况的差异会有不同程度的衰减，并且热流场与山体内部温度场具有一定的关联，可近似为线性关系；从红外遥感图像获得的山体表面红外辐射可以看作是山体各层热流场的叠加；
位于同一海拔高度的地层热场分布大致相同，同样地，根据岩层分布情况，可认为热流场也是相同的，计算每一层的热流场能量，利用通过红外遥感图像得到的山体表面真实红外辐射对不同厚度的山体热流场能量依次进行滤除，可以有效的保护地下建筑信息，得到较为准确的地下建筑扰动信号分布；
喻家山海拔最高为145m，以10m的高度间隔进行分割，向下逐层剥离背景能量；通过不断向下分层，获得滤除10m、20m、30m、40m…80m深度的山体背景热流程能量后的地下建筑扰动信号分布如图15(a)、15(b)、15(c)和15(d)所示；
在探测反演定位过程中，当剥离背景能量至埋深80m时，地下建筑的扰动基本获取完整；剥离背景能量至90m的时候，可以看到此时扰动图像的变化已经很小了，其与上一个高度层的探测反演定位结果相比，扰动区域范围几乎没有发生变化，此时增加的部分主要是80m～90m高度层内的区域干扰，经过上述实验证明，这一类干扰将在逐层的探测反演定位过程中被逐渐削弱；根据可见光图像对山体表面存在的岩石和房屋等干扰进一步滤除；最终得到喻家山区域的探测反演定位结果，如图16所示；
(6)利用膨胀腐蚀加连通域分割算法得到地下建筑位置
根据上述步骤得到的地下建筑探测反演的结果图像，利用膨胀腐蚀操作消除噪声点干扰，根据连通域分割方法结合热传导在时间和空间上的连续性，获得地下建筑的探测位置；
通过形态学运算定位地下建筑的位置；将探测结果图像二值化如图17，然后应用开运算和闭运算消除噪声点的干扰，使细窄处连接断开并消除毛刺，闭运算和开运算结果如图 18；
其中，用半径为6的圆形结构元素进行腐蚀；应用连通域分割，可以根据上图结果分割得到三个区域位置如图19；并在探测结果上显示连通域分割结果，如图20；
由于热传递在时间和空间上具有连续性，判断区域图像重心周围8邻域像素是否都是扰动点，可排除掉非地下建筑区域；地下建筑所在位置区域如图21所示。
在地下建筑区域内，从重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的80％以上，此时的矩形可视为地下建筑的位置如图22。并在可见光图像上标记对应地下建筑位置如图23。
综上所述，本发明与现有技术相比，存在以下优势：
本发明采用分布式地下建筑所在山体与空气层之间热辐射模型，结合DEM数据，计算山体表面太阳辐射总能量和结合高光谱图像计算山体地表植被辐射模型，山体分层迭代计算地质背景热流通量；滤除背景热流量得到地下建筑热传播能量扰动信号；充分考虑外部环境所带给山体的影响，如山体岩层热传播能量、地下建筑热传播扩散能量、太阳辐射能量以及地形因素带给红外数据的扰动差异；以本发明解决了分布式地下建筑在红外图像中经地层传导调制后信号变弱，无法有效反演探测定位等问题，从而实现对分布式地下建筑的反演探测定位。
本发明考虑了山体表面太阳辐射能量，在DEM各个数据栅格上充分考虑太阳直接辐射和散射辐射的影响，相比于多时相红外图像分析地下建筑信号随时间变化规律，建立数学模型计算求解地下建筑位置的方法相比，更加精确的考虑了山体外部环境对地下建筑能量扰动的影响，使得探测反演结果更加精确；相比于对山体红外图像进行统计分析，寻找有地下建筑的山体与旁边下无地下建筑的山体的形成的特征模式，利用对可能含有地下建筑的山体环境的红外图像中进行一定方向的遍历搜寻该特征模式，从而实现分布式地下建筑的探测、定位的方法，对分布式地下建筑反演探测定位的精准度很高。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种滤除分布式地下建筑背景热流通量的反演探测方法，其特征在于，包括以下步骤：
(1)根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；且基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
(2)从地下建筑对应山体的红外遥感图像中滤除山体表面太阳辐射总能量和山体地表植被热辐射能量；
(3)结合DEM数据，在经过步骤(2)处理后的红外图像中分层滤除山体地质背景热流通量，获取扰动信号图像；
(4)利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
(5)根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；
(6)对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
2.根据权利要求1所述的反演探测方法，其特征在于，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述地下建筑块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
3.根据权利要求1或2所述的滤除方法，其特征在于，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影区域；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布；
获取山体表面太阳散射辐射分布的方法，包括以下步骤：
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，分别对若干区域中的太阳散射辐射量进行计算，叠加获取整个扇区的太阳散射辐射量；
将所有栅格对应的所有扇区的太阳散射辐射量叠加，获取所述山体表面太阳散射辐射分布。
4.根据权利要求1所述的滤除方法，其特征在于，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
5.根据权利要求1或4所述的滤除方法，其特征在于，地下建筑的位置信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置。
6.一种滤除分布式地下建筑地质背景热流量的反演探测装置，其特征在于，包括：
太阳总辐射能量分布计算模块，用于根据DEM数据，采用太阳总辐射能量分布模型计算山体表面太阳辐射总能量；
地表植被辐射分布计算模块，用于基于高光谱数据，采用地表植被辐射模型获取山体表面植被热辐射能量；
山体三维温度场计算模块，用于基于DEM数据，采用山体三维温度场模型获取山体地质背景热流通量；
扰动信号图像获取模块，用于从山体的红外图像中滤除山体表面太阳辐射总能量、山体地表植被辐射能量；再结合DEM数据，分层滤除山体地质背景热流通量，获取地下建筑的扰动信号图像；
噪声消除模块，用于利用膨胀腐蚀方法消除扰动信号图像的噪声点干扰后，采用连通域分割方法，将扰动信号图像分割成若干区域图像；
地下建筑位置获取模块，用于根据热传递在时间和空间上的连续性，排除没有地下建筑的区域图像，获取含有地下建筑的位置信息；对含有地下建筑的扰动信号进行反演定位，确定地下建筑的埋深信息；
其中，所述太阳总辐射能量分布模型为通过太阳直接辐射分布和太阳散射辐射分布构建的模型；所述地表植被辐射模型为根据SVAT模型结合高光谱数据构建；所述山体三维温度场模型为基于分布式地下建筑的几何形貌、山体背景热辐射特性和地下建筑的热辐射特性构建的模型。
7.根据权利要求6所述的反演探测装置，其特征在于，所述山体三维温度场模型的建立方法，包括以下步骤：
利用COMSOL有限元仿真软件，结合DEM数据，根据地下建筑耦合模型传热过程中的热量守恒以及地下建筑表面边界温度条件构建地下建筑块状热能辐射模型；
建立山体中岩层向有温度梯度方向发出的热辐射方程，构建山体背景热辐射模型；其中山体中的岩层为不包含地下建筑的岩层；
将所述地下建筑块状热能辐射模型和山体背景热辐射模型叠加，获取含有地下建筑的山体与空气层之间的热辐射模型；
将山体与空气层之间的热辐射模型与DEM数据相结合，构建山体三维温度场模型。
8.根据权利要求6或7所述的滤除装置，其特征在于，获取太阳直接辐射分布的方法，包括以下步骤：
基于太阳方位角、太阳高度角和山体的地形，获取山体阴影；
根据太阳高度角、大气层上界接收的太阳辐射和太阳直接辐射的大气透过率，获取太阳直接辐射的表达式；
将山体阴影区域的太阳直接辐射量设置为0，非阴影区域的太阳直接辐射量通过太阳直接辐射的表达式计算，获取太阳直接辐射分布；
根据DEM数据的各栅格，计算半球视域；
采用地球曲率对所述半球视域进行校正；
将所述半球视域与太阳位置和天空信息相结合，将所述半球视域分割成若干扇区；
将若干所述扇区分割成若干区域，当扇区足够小时，该扇区内的太阳散射辐射可认为是相同，将若干区域中的太阳散射辐射量分别计算，叠加获取整个扇区的太阳散射辐射量；
将所有要计算散射辐射的位置点的视域叠加，从而得到该位置所能接收到的散射辐射分布对可见视域范围内的散射辐射进行累加，得到该观测处的总散射辐射。对山体地区所有位置进行上述操作，获取所述山体表面太阳散射辐射分布。
9.根据权利要求5所述的滤除装置，其特征在于，所述膨胀腐蚀的方法为：将扰动信号图像二值化处理后，应用开运算和闭运算消除噪声点的干扰。
10.根据权利要求5或9所述的滤除装置，其特征在于，地下建筑的位置信息获取方法为：
在区域图像中排除没有地下建筑的区域图像，获取地下建筑区域；
从地下建筑区域的重心点开始做矩形，依次扩大矩形面积，直至包含地下建筑区域内像素总数的预设比例以上，将扩大后的矩形面积作为地下建筑的位置信息。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111680054.3A CN114398812B (zh) | 2021-12-31 | 2021-12-31 | 滤除分布式地下建筑背景热流通量的反演探测方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111680054.3A CN114398812B (zh) | 2021-12-31 | 2021-12-31 | 滤除分布式地下建筑背景热流通量的反演探测方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114398812A true CN114398812A (zh) | 2022-04-26 |
| CN114398812B CN114398812B (zh) | 2024-06-21 |
Family
ID=81229164
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111680054.3A Active CN114398812B (zh) | 2021-12-31 | 2021-12-31 | 滤除分布式地下建筑背景热流通量的反演探测方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114398812B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105426881A (zh) * | 2015-12-24 | 2016-03-23 | 华中科技大学 | 山体背景热场模型约束的地下热源昼间遥感探测定位方法 |
| WO2018028191A1 (zh) * | 2016-08-10 | 2018-02-15 | 福州大学 | 一种基于波段比模型和太阳高度角的tavi计算方法 |
-
2021
- 2021-12-31 CN CN202111680054.3A patent/CN114398812B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105426881A (zh) * | 2015-12-24 | 2016-03-23 | 华中科技大学 | 山体背景热场模型约束的地下热源昼间遥感探测定位方法 |
| WO2018028191A1 (zh) * | 2016-08-10 | 2018-02-15 | 福州大学 | 一种基于波段比模型和太阳高度角的tavi计算方法 |
Non-Patent Citations (1)
| Title |
|---|
| 吕震华;潘晓丽;龚光红;周荣坤;: "地面背景红外辐射特性建模与仿真研究", 中国电子科学研究院学报, no. 04, 20 August 2017 (2017-08-20) * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114398812B (zh) | 2024-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN103942420B (zh) | 一种建筑物尺度的太阳直射辐射能快速估算方法 | |
| KR101943759B1 (ko) | 도시 상세 복사량모의를 통한 고온지역 열저감 위한 활용정보 추출 방법 및 장치 | |
| Jiao et al. | Evaluation of four sky view factor algorithms using digital surface and elevation model data | |
| CN105426881B (zh) | 山体背景热场模型约束的地下热源昼间遥感探测定位方法 | |
| Araya-Muñoz et al. | Assessing the solar potential of roofs in Valparaíso (Chile) | |
| CN110990511B (zh) | 一种顾及城市二维和三维形态的局部气候区分类方法 | |
| CN109871637A (zh) | 一种云天条件下近地面气温估算方法 | |
| Miller et al. | A dynamic enhancement with background reduction algorithm: Overview and application to satellite‐based dust storm detection | |
| Effat | Mapping solar energy potential zones, using SRTM and spatial analysis, application in Lake Nasser Region, Egypt | |
| Duulatov et al. | Current and future trends of rainfall erosivity and soil erosion in Central Asia | |
| CN114330074B (zh) | 一种山体条带状地下隧道反演探测定位方法及装置 | |
| CN114112069B (zh) | 地质约束的城市深埋条带通道红外成像探测方法及系统 | |
| CN114398812B (zh) | 滤除分布式地下建筑背景热流通量的反演探测方法及装置 | |
| CN108132096B (zh) | 一种基于激光雷达的林窗太阳辐射监测方法 | |
| Rojo et al. | Development of a dynamic model to estimate canopy par interception | |
| He et al. | Retrieval of rugged mountainous areas Land Surface Temperature from high-spatial-resolution thermal infrared remote sensing data | |
| Reda et al. | Assessing the effect of volcano-tectonic activity on the spatial distribution of surface temperature in the Main Ethiopian Rift, using integrated approach of remote sensing data | |
| Siombone et al. | Land Cover, Land Surface Temperature and Geomorphology Structure at Tulehu Geothermal Area, Ambon, Indonesia | |
| Chowdhury et al. | Application of remote sensing and GIS in groundwater studies: an overview | |
| Zhang et al. | Analysis of influence mechanism of spatial distribution of incoming solar radiation based on DEM | |
| Samanta et al. | Geospatial modeling of solar radiation to explore solar energy potential in Papua New Guinea | |
| Tootchifatidehi | Development of a global wetland map and application to describe hillslope hydrology in the ORCHIDEE land surface model | |
| Palmer et al. | Detection of roof shading for PV based on LiDAR data using a multi-modal approach | |
| Díaz et al. | Solar radiation maps | |
| Apke et al. | On the origin of rotation derived from super rapid scan satellite imagery at the cloud tops of severe deep convection |
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 |