一种山体条带状地下隧道反演探测定位方法及装置

1.本发明属于地球物理学学科理论与遥感技术交叉领域，更具体地，涉及一种山体条带状地下隧道反演探测定位方法及装置。


背景技术：

2.随着城市现代化的不断推进，为保证土地资源的可持续利用，地下空间的开发利用受到了广泛重视。地下空间是珍贵的自然资源和战略资源，其类型有地下建筑物和地下工程、人工/自然洞穴等，相比于地面设施，其具有更好的隐蔽性、热稳定性和经济性。地下空间的开发主要指人工地下建筑和隧道工程等典型地下目标的建设，在日常生活中随处可见，例如交通隧道、地下通道、地下超市和地下停车场等。针对地下目标探测的研究不仅可以探明地下空间内部情况，还可以为地下空间规划建设提供帮助，促进完善地下设施设计方法。对于地下管线等小型地下目标探测的研究同样具有重要意义，通过探测可以掌握管线走向和分布，及时定位管线故障位置。通过探测私接的地下管线，可以及时发现并制止违法行为，促进城市管理的完善。因此十分有必要开展山体中条带状地下隧道探测反演的研究。
3.目前，国内外利用红外成像技术对山体中条带状地下隧道探测方法的研究采用的手段多为根据多时相红外图像分析目标信号随时间变化规律，建立数学模型计算求解隧道位置。或是利用微波红外增强技术，根据微波功率传输数学模型增强红外图像，突出目标/背景信号强度对比，以此实现对山体中条带状地下隧道位置的探测。
4.现有技术公开了一种山地环境中隧道目标红外成像探测定位方法，通过对山地红外图像进行统计分析，寻找地下有隧道山体与旁边下无隧道山体的形成的特征模式，利用对可能含有隧道的山体环境的红外图像中进行一定方向的遍历搜寻该特征模式，从而实现隧道目标的探测、定位。
5.也有现有技术公开了一种山体中条带状地下隧道的探测定位方法，建立条带状目标所在的仿真山体背景热辐射场，将仿真山体背景热辐射场与目标所在山体红外遥感图进行线性映射，得到山体背景热辐射模型。从目标所在山体红外遥感图中滤除仿真山体背景热辐射场，得到条带状目标信息场。
6.但是上述方法局限性在于：1、所需红外遥感数据样本量大，而且获取难度较高2、方法难以适用于对其他类型目标的探测，环境状况和背景特性对方法的探测结果具有一定影响。3、仿真计算出的是山体背景热辐射场，而真实环境中隧道所产生是扰动是以热流的形式传导。4、没有充分考虑外部环境所带给隧道的影响，如地表植被、太阳辐射以及地形因素带给红外数据的扰动差异。


技术实现要素：

7.针对现有技术的缺陷，本发明的目的在于提供一种山体条带状地下隧道反演探测定位方法及装置，旨在解决现有的隧道在红外图像中经地层传导调制后信号变弱，无法有
效反演探测定位的问题。
8.为实现上述目的，一方面，本发明提供了一种山体条带状地下隧道反演探测定位方法，包括以下步骤：
9.(1)采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量，迭代计算山体的各层背景热流场能量；结合高光谱数据，计算山体背景热传播能量；
10.(2)在红外遥感图像中滤除山体表面太阳辐射能量和山体背景热传播能量；
11.(3)在步骤(2)的基础上，采用地下目标反演模型，在红外遥感图像中滤除山体的各层背景热流场能量，获取山体中条带状地下隧道的最优化海拔高度和各层山体中条带状地下隧道热流场能量构建的扰动信号分布图像；
12.(4)利用霍夫变换检测方法在所述扰动信号分布图像中检测直线，根据隧道工程设计的关联性原则，对条带状不连续的扰动信号图进行处理，拟合得到隧道目标的走向及形貌；
13.其中，所述山体与空气层之间热辐射模型包括隧道沿线热能辐射模型、太阳总辐射能量分布模型和山体背景热传播模型；
14.所述隧道沿线热能辐射模型为利用comsol有限元仿真软件，结合dem数据和气象参数构建，用于计算隧道沿线热能辐射分布；
15.所述太阳总辐射能量分布模型为利用半球视域方法，结合所述dem数据和山体隧道地理位置信息，计算山体表面太阳辐射能量的模型；
16.所述山体背景热传播模型为根据svat模型结合高光谱数据，计算地表植被红外辐射分布的模型。
17.进一步优选地，构建所述隧道沿线热能辐射模型的方法为：
18.基于dem数据，利用comsol有限元仿真软件，由山体隧道耦合模型传热过程中的能量守恒方程式，结合隧道表面边界条件、隧道外表面与山体岩层的接触面存在的热平衡以及引入形状因子的导热流量构建。
19.进一步优选地，所述山体表面太阳辐射能量的计算方法，包括以下步骤：
20.根据dem数据的各栅格计算半球视域；
21.将所述半球视域与太阳位置和天空信息相结合，根据栅格处的太阳直接辐射量和太阳散射辐射量，获取栅格处太阳直接辐射分布和太阳散射辐射分布；
22.将所述栅格处太阳直接辐射分布和太阳散射辐射分布进行叠加，获取栅格处地表接收的太阳辐射总量；
23.将所有所述栅格处地表接收的太阳辐射总量相加，得到山体表面太阳辐射能量。
24.进一步优选地，所述步骤(3)中各层背景热流场能量计算方法为：
[0025][0026]
其中，hf(x,y,z,t)为山体中单位面积的热流量，a(z)和ε(z)分别为山体内部某一海拔高度处的横截面积和热流的衰减系数；z1(x,y)和z2(x,y)分别为山体表面高度坐标和从表面向下某一深度的高度坐标，通过调整z2(x,y)实现计算不同厚度的山体背景包含的热流场能量；
[0027]
单位面积的热流量如下：
[0028][0029]
其中，cm为岩层介质的容积热容量，k为岩层介质的热导率，为垂直方向的温度梯度。
[0030]
进一步优选地，所述步骤(3)具体包括以下步骤：
[0031]
a.根据山体地表的实际总辐射量和红外遥感图像上的总辐射量，计算山体总辐射量误差；
[0032][0033]
其中,ε表示误差范围，是一个极小的正数eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；在红外遥感成像过程中，传感器和外界环境的影响会使测量的地表真实辐射量存在误差，因而利用红外遥感图像获得的山体表面辐射量是地表真实辐射量与误差项的叠加，如式：
[0034]
eti＝emi+σ
[0035]
其中，emi表示第i个像素点的山体地表的实际总辐射量，σ表示测量误差项。
[0036]
将测量误差纳入考虑的范围内后，收敛条件可表示为：
[0037][0038]
b.以隧道所在海拔高度ak为待反演变量，令k的初始值为1，初始化a1；
[0039]
c.采用当前红外遥感图像上各像素点总辐射量减去各像素点对应山体在海拔高度ak处的背景热流场能量求平方和，获取山体在海拔高度ak处的条带状地下隧道热流场能量；
[0040][0041]
其中，ak和n分别表示隧道海拔高度和计算区域的像素点数，eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；
[0042]
d.采用当前迭代的隧道所在海拔高度ak、当前迭代的搜索步长和搜索方向，确定下一次迭代的海拔高度a
k+1
；其中，当前迭代的搜索方向为利用上一次迭代的搜索方向、山体在海拔高度ak处的条带状地下隧道热流场能量和共轭参数计算获取；
[0043]ak+1
＝ak+λ
kdk
[0044]
其中，ak表示隧道海拔高度a的第k次运算值；λk表示搜索步长，dk表示搜索方向；
[0045][0046]
其中，表示山体在海拔高度ak处的条带状地下隧道热流场能量的梯度，可采用中心差分近似计算，βk表示共轭参数；
[0047]
e.将下一次迭代的海拔高度a
k+1
作为当前迭代的海拔高度，返回至步骤c，直至山
体在海拔高度ak处的条带状地下隧道热流场能量与山体总辐射量误差相等，则结束迭代；
[0048]
f.根据每次迭代中的山体中条带状地下隧道热流场能量，构建扰动信号分布图像；且将海拔高度ak作为山体中条带地下隧道的最优化海拔高度。
[0049]
进一步优选地，所述太阳直接辐射量e
dir
的计算方法为：
[0050]edir
＝∑e
dir
(θ,α)
[0051]edir
(θ,α)＝s
const
*τ
bm(θ)
*sd
θ,α
*sung
θ,α
*cos(angi(θ,α))
[0052][0053]
angi(θ,α)＝arccos(cos(θ)*cos(gz)+sin(θ)*sin(gz)*cos(α-g
α
))
[0054]
其中，s
const
为太阳常数；τb为最短路径上的大气透过率；m(θ)为光学路径长度；sd
θ,α
为持续时间，sung
θ,α
为太阳图扇区的孔隙度；angi(θ,α)为天空扇区质心与山体表面垂线入射角度；gz和g
α
分别为山体表面的天顶角和方位角，θ1和θ2分别为天空扇区两边界处天顶角；h为海拔高程值；
[0055]
所述太阳散射辐射量e
dif
的计算方法为：
[0056]edif
＝∑e
dif
(θ,α)
[0057]edif
(θ,α)＝e
glb
*p*tc*skyg
θ,α
*w(θ,α)*cos(angi(θ,α))
[0058]eglb
＝(s
const
∑(τ
bm(θ)
))/(1-p)
[0059]
w(θ,α)＝(cosθ
2-cosθ1)/div
azi
[0060]
其中，s
const
为太阳常数；angi(θ,α)为天空扇区质心与表面垂线入射角度，e
glb
为正常总辐射量，p为散射辐射通量所占的比例，skyg
θ,α
为天空扇区的孔隙度，h为海拔高程值；θ1和θ2分别为天空扇区两边界处天顶角；div
azi
为天空图方位分割数量；w(θ,α)天空扇区散射辐射量和总扇区散射辐射量的比值；tc为太阳散射辐射的持续时间；τb为最短路径上的大气透过率；e
dir
(θ,α)和e
dif
(θ,α)分别为观测点所在的质心位于天顶角为θ和方位角为α的直接辐射能量和散射辐射能量。
[0061]
另一方面，本发明提供了一种山体条带状地下隧道反演探测定位装置，包括：
[0062]
山体背景热流场计算模块，用于采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量，迭代计算山体的各层背景热流场能量；结合高光谱数据，计算山体背景热传播能量；
[0063]
扰动信号分布图像构建模块，用于在红外遥感图像中滤除山体表面太阳辐射能量、山体背景热传播能量；再采用地下目标反演模型，在红外遥感图像中滤除山体的各层背景热流场能量，获取山体中条带状地下隧道的最优化海拔高度和各层山体中条带状地下隧道热流场能量构建的扰动信号分布图像；
[0064]
目标探测位置拟合模块，用于利用霍夫变换检测方法在所述扰动信号分布图像中检测直线，根据隧道工程设计的关联性原则，对山体中条带状地下隧道不连续的扰动信号进行处理，拟合得到隧道的走向及形貌；
[0065]
其中，山体背景热流场计算模块包括隧道沿线热能辐射分布构建单元、太阳总辐射能量分布构建单元和山体背景热传播分布构建单元；
[0066]
所述隧道沿线热能辐射分布构建单元，用于利用comsol有限元仿真软件，结合dem数据和气象参数构建隧道沿线热能辐射模型；
[0067]
所述太阳总辐射能量分布构建单元，用于利用半球视域方法，结合所述dem数据和
山体隧道地理位置信息，构建太阳总辐射能量分布模型，所述太阳总辐射能量分布模型用于计算山体表面太阳辐射能量；
[0068]
所述山体背景热传播分布构建单元，用于根据svat模型结合高光谱数据，构建山体背景热传播模型。
[0069]
进一步优选地，构建所述隧道沿线热能辐射模型的方法为：
[0070]
利用comsol有限元仿真软件，基于dem数据，由山体隧道耦合模型传热过程中的能量守恒方程式，结合隧道表面边界条件、隧道外表面与山体岩层的接触面存在的热平衡以及引入形状因子的导热流量构建。
[0071]
进一步优选地，所述山体表面太阳辐射能量的计算方法，包括以下步骤：
[0072]
根据dem数据的各栅格计算半球视域；
[0073]
将所述半球视域与太阳位置和天空信息相结合，根据栅格处的太阳直接辐射能量和太阳散射辐射能量，获取栅格处太阳直接辐射分布和太阳散射辐射分布；
[0074]
将所述栅格处太阳直接辐射分布和太阳散射辐射分布进行叠加，获取栅格处地表接收的太阳辐射总能量；
[0075]
将所有所述栅格处地表接收的太阳辐射总能量相加，得到山体表太阳辐射能量。
[0076]
进一步优选地，获取山体中各层背景热流场能量计算方法为：
[0077][0078]
其中，hf(x,y,z,t)为山体中单位面积的热流量，a(z)和ε(z)分别为山体内部某一海拔高度处的横截面积和热流的衰减系数；z1(x,y)和z2(x,y)分别为山体表面高度坐标和从表面向下某一深度的高度坐标，通过调整z2(x,y)实现计算不同厚度的山体背景包含的热流场能量；
[0079]
单位面积的热流量如下：
[0080][0081]
其中，cm为岩层介质的容积热容量，k为岩层介质的热导率，为垂直方向的温度梯度。
[0082]
进一步优选地，获取山体中条带状地下隧道的最优化海拔高度和干扰信号分布图像的方法，具体包括以下步骤：
[0083]
a.根据山体地表的实际总辐射量和红外遥感图像上的总辐射量，计算山体总辐射量误差；
[0084][0085]
其中,ε表示误差范围，是一个极小的正数eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；在红外遥感成像过程中，传感器和外界环境的影响会使测量的地表真实辐射量存在误差，因而利用红外遥感图像获得的山体表面辐射量是地表真实辐射量与误差项的叠加，如式：
[0086]
eti＝emi+σ
[0087]
其中，emi表示第i个像素点的山体地表的实际总辐射量，σ表示测量误差项。
[0088]
将测量误差纳入考虑的范围内后，收敛条件可表示为：
[0089][0090]
b.以隧道所在海拔高度ak为待反演变量，令k的初始值为1，初始化a1；
[0091]
c.采用当前红外遥感图像上各像素点总辐射量减去各像素点对应山体在海拔高度ak处的背景热流场能量求平方和，获取山体在海拔高度ak处的条带状地下隧道热流场能量；
[0092][0093]
其中，ak和n分别表示隧道海拔高度和计算区域的像素点数，eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；
[0094]
d.采用当前迭代的隧道所在海拔高度ak、当前迭代的搜索步长和搜索方向，确定下一次迭代的海拔高度a
k+1
；其中，当前迭代的搜索方向为利用上一次迭代的搜索方向、山体在海拔高度ak处的条带状地下隧道热流场能量和共轭参数计算获取；
[0095]ak+1
＝ak+λ
kdk
[0096]
其中，ak表示隧道海拔高度a的第k次运算值；λk表示搜索步长，dk表示搜索方向；
[0097][0098]
其中，表示山体在海拔高度ak处的条带状地下隧道热流场能量的梯度，可采用中心差分近似计算，βk表示共轭参数；
[0099]
e.将下一次迭代的海拔高度a
k+1
作为当前迭代的海拔高度，返回至步骤c，直至山体在海拔高度ak处的条带状地下隧道热流场能量与山体总辐射量误差相等，则结束迭代；
[0100]
f.根据每次迭代中的山体中条带状地下隧道热流场能量，构建扰动信号分布图像；且将海拔高度ak作为山体中条带地下隧道的最优化海拔高度。
[0101]
进一步优选地，所述太阳直接辐射量e
dir
的计算方法为：
[0102]edir
＝∑e
dir
(θ,α)
[0103]edir
(θ,α)＝s
const
*τ
bm(θ)
*sd
θ,α
*sung
θ,α
*cos(angi(θ,α))
[0104][0105]
angi(θ,α)＝arccos(cos(θ)*cos(gz)+sin(θ)*sin(gz)*cos(α-g
α
))
[0106]
其中，s
const
为太阳常数；τb为最短路径上的大气透过率；m(θ)为光学路径长度；sd
θ,α
为持续时间，sung
θ,α
为太阳图扇区的孔隙度；angi(θ,α)为天空扇区质心与山体表面垂线入射角度；gz和g
α
分别为山体表面的天顶角和方位角，θ1和θ2分别为天空扇区两边界处天顶角；h为海拔高程值；
[0107]
所述太阳散射辐射能量e
dif
的计算方法为：
[0108]edif
＝∑e
dif
(θ,α)
[0109]edif
(θ,α)＝e
glb
*p*tc*skyg
θ,α
*w(θ,α)*cos(angi(θ,α))
[0110]eglb
＝(s
const
∑(τ
bm(θ)
))/(1-p)
[0111]
w(θ,α)＝(cosθ
2-cosθ1)/div
azi
[0112]
其中，s
const
为太阳常数；angi(θ,α)为天空扇区质心与表面垂线入射角度，e
glb
为正常总辐射能量，p为散射辐射通量所占的比例，skyg
θ,α
为天空扇区的孔隙度，h为海拔高程值；θ1和θ2分别为天空扇区两边界处天顶角；div
azi
为天空图方位分割数量；w(θ,α)为天空扇区散射辐射量和总扇区散射辐射量的比值；tc为太阳散射辐射的持续时间；τb为最短路径上的大气透过率；e
dir
(θ,α)和e
dif
(θ,α)分别为观测点所在的质心位于天顶角为θ和方位角为α的直接辐射量和散射辐射量。
[0113]
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
[0114]
本发明采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量、迭代计算山体的各层山体背景热流场、隧道和山体背景热传播能量；充分考虑外部环境所带给山体的影响，如地表植被、太阳辐射以及地形因素带给红外数据的扰动差异；以本发明解决了山体中条带状地下隧道在红外图像中经地层传导调制后信号变弱，无法有效反演探测定位等问题，从而实现对山地环境中隧道的反演探测定位。
[0115]
本发明考虑了山体表面太阳辐射能量，在dem各个数据栅格上充分考虑太阳直接辐射和散射辐射的影响，相比于多时相红外图像分析目标信号随时间变化规律，建立数学模型计算求解隧道位置的方法相比，更加精确的考虑了外部环境(太阳直接辐射、散射辐射、地表植被等)对山体中条带状地下隧道能量扰动的影响，使得探测反演结果更加精确，相对于目标区域多时相红外数据的需求，此方法对数据获取更加容易；相比于对山地红外图像进行统计分析，寻找地下有隧道山体与旁边下无隧道山体的形成的特征模式，利用对可能含有隧道的山体环境的红外图像中进行一定方向的遍历搜寻该特征模式，从而实现隧道目标的探测、定位的方法，本发明无需大量选取山体中有无隧道的样本段，此方法效率更高；多源信息融合滤除环境干扰，使得山体中条带状地下隧道反演探测定位的精准度很高。
附图说明
[0116]
图1是本发明实施例提供的滤除山体中条带状地下隧道环境干扰的反演探测定位方法流程图；
[0117]
图2是本发明实施例提供的探测反演方法流程图；
[0118]
图3是本发明实施例提供的潭峪沟地区的dem数据；
[0119]
图4是本发明实施例提供的山体背景温度场分布情况示意图；
[0120]
图5是本发明实施例提供的观测点所在位置示意图；
[0121]
图6是本发明实施例提供的观测点位置的半球视域；
[0122]
图7是本发明实施例提供的观测点位置接收到的太阳直接辐射分布示意图；
[0123]
图8是本发明实施例提供的潭峪沟地区地表接收到的太阳直接辐射能量分布图；
[0124]
图9是本发明实施例提供的观测点位置接收到的太阳散射辐照分布示意图；
[0125]
图10是本发明实施例提供的潭峪沟地区地表接收到的太阳散射辐照能量分布图；
[0126]
图11是本发明实施例提供的潭峪沟地区地表接收到的太阳总辐照能量分布图；
[0127]
图12是本发明实施例提供的潭峪沟隧道区域山体植被仿真红外辐照图；
[0128]
图13(a)是本发明实施例提供的滤除40m山体厚度所包含热流畅能量的结果图像；
[0129]
图13(b)是本发明实施例提供的滤除80m山体厚度所包含热流畅能量的结果图像；
[0130]
图13(c)是本发明实施例提供的滤除140m山体厚度所包含热流畅能量的结果图像；
[0131]
图13(d)是本发明实施例提供的滤除180m山体厚度所包含热流畅能量的结果图像；
[0132]
图13(e)是本发明实施例提供的滤除200m山体厚度所包含热流畅能量的结果图像；
[0133]
图14是本发明实施例提供的潭峪沟隧道区域的反演探测定位结果图；
[0134]
图15是本发明实施例提供的潭峪沟隧道区域直线检测结果示意图；
[0135]
图16是本发明实施例提供的潭峪沟隧道最终检测结果和真实位置比较。
具体实施方式
[0136]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0137]
本发明提供的滤除环境干扰的山体中条带状地下隧道反演探测定位方法，包括以下步骤：
[0138]
利用红外图像对应区域的dem数据，计算山体表面接收到的太阳辐射总能量；
[0139]
结合dem数据、山体隧道中介质材料特性参数和气象参数，计算隧道沿线热能辐射分布；
[0140]
利用svat模型结合高光谱图像数据，计算山体背景热传播能量分布；
[0141]
实现滤除环境干扰的条带状目标反演探测定位。
[0142]
具体如下：
[0143]
一方面，本发明提供了一种山体条带状地下隧道反演探测定位方法，包括以下步骤：
[0144]
(1)采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量，迭代计算山体的各层背景热流场能量；结合高光谱数据，计算山体背景热传播能量；
[0145]
(2)在红外遥感图像中滤除山体表面太阳辐射能量和山体背景热传播能量；
[0146]
(3)在步骤(2)的基础上，采用地下目标反演模型，在红外遥感图像中滤除山体的各层背景热流场能量，获取山体中条带状地下隧道的最优化海拔高度和各层山体中条带状地下隧道热流场能量构建的扰动信号分布图像；
[0147]
(4)利用霍夫变换检测方法在所述扰动信号分布图像中检测直线，根据隧道工程设计的关联性原则，对条带状不连续的扰动信号图进行处理，拟合得到隧道目标的走向及形貌；
[0148]
其中，所述山体与空气层之间热辐射模型包括隧道沿线热能辐射模型、太阳总辐射能量分布模型和山体背景热传播模型；
[0149]
所述隧道沿线热能辐射模型为利用comsol有限元仿真软件，结合dem数据和气象
参数构建，用于计算隧道沿线热能辐射分布；
[0150]
所述太阳总辐射能量分布模型为利用半球视域方法，结合所述dem数据和山体隧道地理位置信息，计算山体表面太阳辐射能量的模型；
[0151]
所述山体背景热传播模型为根据svat模型结合高光谱数据，计算地表植被红外辐射分布的模型。
[0152]
进一步优选地，构建所述隧道沿线热能辐射模型的方法为：
[0153]
利用comsol有限元仿真软件，结合dem数据，由山体隧道耦合模型传热过程中的能量守恒方程式，结合隧道表面边界条件、隧道外表面与山体岩层的接触面存在的热平衡以及引入形状因子的导热流量构建。
[0154]
进一步优选地，所述山体表面太阳辐射能量的计算方法，包括以下步骤：
[0155]
根据dem数据的各栅格计算半球视域；
[0156]
将所述半球视域与太阳位置和天空信息相结合，根据栅格处的太阳直接辐射量和太阳散射辐射量，获取栅格处太阳直接辐射分布和太阳散射辐射分布；
[0157]
将所述栅格处太阳直接辐射分布和太阳散射辐射分布进行叠加，获取栅格处地表接收的太阳辐射总量；
[0158]
将所有所述栅格处地表接收的太阳辐射总量相加，得到山体表面太阳辐射能量。
[0159]
进一步优选地，所述步骤(3)中各层背景热流场能量计算方法为：
[0160][0161]
其中，hf(x,y,z,t)为山体中单位面积的热流量，a(z)和ε(z)分别为山体内部某一海拔高度处的横截面积和热流的衰减系数；z1(x,y)和z2(x,y)分别为山体表面高度坐标和从表面向下某一深度的高度坐标，通过调整z2(x,y)实现计算不同厚度的山体背景包含的热流场能量；
[0162]
单位面积的热流量如下：
[0163][0164]
其中，cm为岩层介质的容积热容量，k为岩层介质的热导率，为垂直方向的温度梯度。
[0165]
进一步优选地，所述步骤(3)具体包括以下步骤：
[0166]
a.根据山体地表的实际总辐射量和红外遥感图像上的总辐射量，计算山体总辐射量误差；
[0167][0168]
其中,ε表示误差范围，是一个极小的正数eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；在红外遥感成像过程中，传感器和外界环境的影响会使测量的地表真实辐射量存在误差，因而利用红外遥感图像获得的山体表面辐射量是地表真实辐射量与误差项的叠加，如式：
[0169]
eti＝emi+σ
[0170]
其中，emi表示第i个像素点的山体地表的实际总辐射量，σ表示测量误差项。
[0171]
将测量误差纳入考虑的范围内后，收敛条件可表示为：
[0172][0173]
b.以隧道所在海拔高度ak为待反演变量，令k的初始值为1，初始化a1；
[0174]
c.采用当前红外遥感图像上各像素点总辐射量减去各像素点对应山体在海拔高度ak处的背景热流场能量求平方和，获取山体在海拔高度ak处的条带状地下隧道热流场能量；
[0175][0176]
其中，ak和n分别表示隧道海拔高度和计算区域的像素点数，eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；
[0177]
d.采用当前迭代的隧道所在海拔高度ak、当前迭代的搜索步长和搜索方向，确定下一次迭代的海拔高度a
k+1
；其中，当前迭代的搜索方向为利用上一次迭代的搜索方向、山体在海拔高度ak处的条带状地下隧道热流场能量和共轭参数计算获取；
[0178]ak+1
＝ak+λ
kdk
[0179]
其中，ak表示隧道海拔高度a的第k次运算值；λk表示搜索步长，dk表示搜索方向；
[0180][0181]
其中，表示山体在海拔高度ak处的条带状地下隧道热流场能量的梯度，可采用中心差分近似计算，βk表示共轭参数；
[0182]
e.将下一次迭代的海拔高度a
k+1
作为当前迭代的海拔高度，返回至步骤c，直至山体在海拔高度ak处的条带状地下隧道热流场能量与山体总辐射量误差相等，则结束迭代；
[0183]
f.根据每次迭代中的山体中条带状地下隧道热流场能量，构建扰动信号分布图像；且将海拔高度ak作为山体中条带地下隧道的最优化海拔高度。
[0184]
进一步优选地，所述太阳直接辐射量e
dir
的计算方法为：
[0185]edir
＝∑e
dir
(θ,α)
[0186]edir
(θ,α)＝s
const
*τ
bm(θ)
*sd
θ,α
*sung
θ,α
*cos(angi(θ,α))
[0187][0188]
angi(θ,α)＝arccos(cos(θ)*cos(gz)+sin(θ)*sin(gz)*cos(α-g
α
))
[0189]
其中，s
const
为太阳常数；τb为最短路径上的大气透过率；m(θ)为光学路径长度；sd
θ,α
为持续时间，sung
θ,α
为太阳图扇区的孔隙度；angi(θ,α)为天空扇区质心与山体表面垂线入射角度；gz和g
α
分别为山体表面的天顶角和方位角，θ1和θ2分别为天空扇区两边界处天顶角；h为海拔高程值；
[0190]
所述太阳散射辐射量e
dif
的计算方法为：
[0191]edif
＝∑e
dif
(θ,α)
[0192]edif
(θ,α)＝e
glb
*p*tc*skyg
θ,α
*w(θ,α)*cos(angi(θ,α))
[0193]eglb
＝(s
const
∑(τ
bm(θ)
))/(1-p)
[0194]
w(θ,α)＝(cosθ
2-cosθ1)/div
azi
[0195]
其中，s
const
为太阳常数；angi(θ,α)为天空扇区质心与表面垂线入射角度，e
glb
为正常总辐射量，p为散射辐射通量所占的比例，skyg
θ,α
为天空扇区的孔隙度，h为海拔高程值；θ1和θ2分别为天空扇区两边界处天顶角；div
azi
为天空图方位分割数量；w(θ,α)天空扇区散射辐射量和总扇区散射辐射量的比值；tc为太阳散射辐射的持续时间；τb为最短路径上的大气透过率；e
dir
(θ,α)和e
dif
(θ,α)分别为观测点所在的质心位于天顶角为θ和方位角为α的直接辐射能量和散射辐射能量。
[0196]
另一方面，本发明提供了一种山体条带状地下隧道反演探测定位装置，包括：
[0197]
山体背景热流场计算模块，用于采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量，迭代计算山体的各层背景热流场能量；结合高光谱数据，计算山体背景热传播能量；
[0198]
扰动信号分布图像构建模块，用于在红外遥感图像中滤除山体表面太阳辐射能量和山体背景热传播能量；再采用地下目标反演模型，在红外遥感图像中滤除山体的各层背景热流场能量，获取山体中条带状地下隧道的最优化海拔高度和各层山体中条带状地下隧道热流场能量构建的扰动信号分布图像；
[0199]
目标探测位置拟合模块，用于利用霍夫变换检测方法在所述扰动信号分布图像中检测直线，根据隧道工程设计的关联性原则，对山体中条带状地下隧道不连续的扰动信号进行处理，拟合得到隧道的走向及形貌；
[0200]
其中，山体背景热流场计算模块包括隧道沿线热能辐射分布构建单元、太阳总辐射能量分布构建单元和山体背景热传播分布构建单元；
[0201]
所述隧道沿线热能辐射分布构建单元，用于利用comsol有限元仿真软件，结合dem数据和气象参数构建隧道沿线热能辐射模型；
[0202]
所述太阳总辐射能量分布构建单元，用于利用半球视域方法，结合所述dem数据和山体隧道地理位置信息，构建太阳总辐射能量分布模型，所述太阳总辐射能量分布模型用于计算山体表面太阳辐射能量；
[0203]
所述山体背景热传播分布构建单元，用于根据svat模型结合高光谱数据，构建山体背景热传播模型。
[0204]
进一步优选地，构建所述隧道沿线热能辐射模型的方法为：
[0205]
利用comsol有限元仿真软件，基于dem数据，由山体隧道耦合模型传热过程中的能量守恒方程式，结合隧道表面边界条件、隧道外表面与山体岩层的接触面存在的热平衡以及引入形状因子的导热流量构建。
[0206]
进一步优选地，所述山体表面太阳辐射能量的计算方法，包括以下步骤：
[0207]
根据dem数据的各栅格计算半球视域；
[0208]
将所述半球视域与太阳位置和天空信息相结合，根据栅格处的太阳直接辐射能量和太阳散射辐射能量，获取栅格处太阳直接辐射分布和太阳散射辐射分布；
[0209]
将所述栅格处太阳直接辐射分布和太阳散射辐射分布进行叠加，获取栅格处地表接收的太阳辐射总能量；
[0210]
将所有所述栅格处地表接收的太阳辐射总能量相加，得到山体表太阳辐射能量。
[0211]
进一步优选地，获取山体中条带状地下隧道的最优化海拔高度和干扰信号分布图像的方法，具体包括以下步骤：
[0212]
a.根据山体地表的实际总辐射量和红外遥感图像上的总辐射量，计算山体总辐射量误差；
[0213][0214]
其中,ε表示误差范围，是一个极小的正数eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；在红外遥感成像过程中，传感器和外界环境的影响会使测量的地表真实辐射量存在误差，因而利用红外遥感图像获得的山体表面辐射量是地表真实辐射量与误差项的叠加，如式：
[0215]
eti＝emi+σ
[0216]
其中，emi表示第i个像素点的山体地表的实际总辐射量，σ表示测量误差项。
[0217]
将测量误差纳入考虑的范围内后，收敛条件可表示为：
[0218][0219]
b.以隧道所在海拔高度ak为待反演变量，令k的初始值为1，初始化a1；
[0220]
c.采用当前红外遥感图像上各像素点总辐射量减去各像素点对应山体在海拔高度ak处的背景热流场能量求平方和，获取山体在海拔高度ak处的条带状地下隧道热流场能量；
[0221][0222]
其中，ak和n分别表示隧道海拔高度和计算区域的像素点数，eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；
[0223]
d.采用当前迭代的隧道所在海拔高度ak、当前迭代的搜索步长和搜索方向，确定下一次迭代的海拔高度a
k+1
；其中，当前迭代的搜索方向为利用上一次迭代的搜索方向、山体在海拔高度ak处的条带状地下隧道热流场能量和共轭参数计算获取；
[0224]ak+1
＝ak+λ
kdk
[0225]
其中，ak表示隧道海拔高度a的第k次运算值；λk表示搜索步长，dk表示搜索方向；
[0226][0227]
其中，表示山体在海拔高度ak处的条带状地下隧道热流场能量的梯度，可采用中心差分近似计算，βk表示共轭参数；
[0228]
e.将下一次迭代的海拔高度a
k+1
作为当前迭代的海拔高度，返回至步骤c，直至山体在海拔高度ak处的条带状地下隧道热流场能量与山体总辐射量误差相等，则结束迭代；
[0229]
f.根据每次迭代中的山体中条带状地下隧道热流场能量，构建扰动信号分布图像；且将海拔高度ak作为山体中条带地下隧道的最优化海拔高度。
[0230]
进一步优选地，所述太阳直接辐射量e
dir
的计算方法为：
[0231]edir
＝∑e
dir
(θ,α)
[0232]edir
(θ,α)＝s
const
*τ
bm(θ)
*sd
θ,α
*sung
θ,α
*cos(angi(θ,α))
[0233][0234]
angi(θ,α)＝arccos(cos(θ)*cos(gz)+sin(θ)*sin(gz)*cos(α-g
α
))
[0235]
其中，s
const
为太阳常数；τb为最短路径上的大气透过率；m(θ)为光学路径长度；sd
θ,α
为持续时间，sung
θ,α
为太阳图扇区的孔隙度；angi(θ,α)为天空扇区质心与山体表面垂线入射角度；gz和g
α
分别为山体表面的天顶角和方位角，θ1和θ2分别为天空扇区两边界处天顶角；h为海拔高程值；
[0236]
所述太阳散射辐射能量e
dif
的计算方法为：
[0237]edif
＝∑e
dif
(θ,α)
[0238]edif
(θ,α)＝e
glb
*p*tc*skyg
θ,α
*w(θ,α)*cos(angi(θ,α))
[0239]eglb
＝(s
const
∑(τ
bm(θ)
))/(1-p)
[0240]
w(θ,α)＝(cosθ
2-cosθ1)/div
azi
[0241]
其中，s
const
为太阳常数；angi(θ,α)为天空扇区质心与表面垂线入射角度，e
glb
为正常总辐射能量，p为散射辐射通量所占的比例，skyg
θ,α
为天空扇区的孔隙度，h为海拔高程值；θ1和θ2分别为天空扇区两边界处天顶角；div
azi
为天空图方位分割数量；w(θ,α)为天空扇区散射辐射量和总扇区散射辐射量的比值；tc为太阳散射辐射的持续时间；τb为最短路径上的大气透过率；e
dir
(θ,α)和e
dif
(θ,α)分别为观测点所在的质心位于天顶角为θ和方位角为α的直接辐射量和散射辐射量。
[0242]
实施例
[0243]
如图1所示，本发明提供了一种滤除环境干扰的山体条带状地下隧道反演探测定位方法，包括以下步骤：
[0244]
(1)隧道沿线热能辐射分布模型的建立
[0245]
以地球物理学、传热学、感技术等学科的基本理论为知道，分析地下隧道与埋藏环境的物质/能量特征关系及相互作用规律，建立目标区的隧道沿线热能辐射模型；更为具体地，包括以下步骤：
[0246]
结合dem数据(如图3所示)，利用现有的软件对隧道线条几何建模仿真，为了将隧道沿线热能辐射信息显示，利用comsol有限元仿真软件对隧道沿线热能辐射场进行仿真计算；
[0247]
(1.1)山体隧道耦合模型传热过程中的热量守恒可描述为：
[0248]
qd+qv＝q
[0249][0250][0251]
其中，qv和qd分别表示对流热通量和传导热通量；q为外界输入的热通量；ρ和c
ρ
分别表示介质的密度和恒压比热容；w为流体速度；k为介质的热导率；h为海拔梯度；为温度在山体隧道的垂直方向上的梯度变化，计算时可采用差分近似，这里的
介质为山体隧道中包含的介质；
[0252]
(1.2)隧道内部温度恒定作为隧道表面边界条件，可表示为：
[0253]
t
in
＝f(y)＝t0[0254]
(1.3)对于隧道外表面与山体岩层的接触面，存在如下热平衡关系：
[0255][0256]
其中，t
out
和t
in
分别表示隧道外表面和内部的温度；ts表示山体地表温度；λc和δc分别表示隧道内外表面之间混凝土的厚度和导热系数，λg和δg分别表示隧道外表面到地表之间花岗岩的厚度和导热系数；
[0257]
(1.4)隧道的纵向长度远大于自身的横截面半径，可引入形状因子简化计算过程，引入形状因子后的导热量计算如下：
[0258]
φ＝λs(t
1-t2)
[0259][0260]
其中，φ为隧道的导热量；λ为等温面之间介质材料的导热系数；t1为上表面导热热量；t2为下表面导热热量；d和l分别表示隧道横截面半径和整体长度，h表示隧道的埋深；
[0261]
通过联立求解公式可得到山体隧道耦合模型的表面温度场分布，本发明结合相关参数利用comsol软件进行建模仿真计算，得到最终的潭峪沟隧道区域山体三维温度场分布情况如图4所示；
[0262]
(2)太阳总辐射能量分布模型的构建
[0263]
在以上过程中通过建模仿真计算得到的山体表面温度数据，是没有考虑外界环境影响的理想情况下的理论数据，与利用红外遥感图像获得的山体表面真实红外辐射数据相比具有一定的差异；因此，需要对太阳辐射和地表植被分布这两大外界干扰因素进行研究计算并去除，以保证探测结果的真实性和有效性；在此利用半球视域方法计算地表接收到的太阳辐射；
[0264]
地表接收到的太阳辐射主要包括直接辐射、散射辐射以及反射辐射，由于反射辐射占比小且仅在特殊条件下才会对地表接收到的太阳辐射产生较大影响，通常太阳辐射计算模型中不考虑反射辐射部分；因此，利用红外手段探测地下目标时要充分考虑太阳辐射的散射辐射和直接辐射影响；
[0265]
本实施例采用半球视域算法计算地表接收到的太阳辐射，该算法的基本思想为：首先根据dem数据(digital elevation model，以绝对高程或海拔表示的地形模型)栅格计算半球视域；然后，将半球视域与太阳位置和天空信息结合起来，判断栅格所受的直接辐射和散射辐射情况；最后，逐个对dem所有数据栅格进行计算得到该区域地表所接收的太阳辐射总量；
[0266]
太阳总辐射量由太阳直接辐射量和散射辐射量相加得到：
[0267]etotal
＝e
dir
+e
dif
[0268]
其中，e
total
为太阳辐射总量，e
dir
和e
dif
分别为直接辐射总量和散射辐射总量；
[0269]
太阳直接辐射量是所有太阳图扇区中直接辐射的总和；散射辐射量是所有天空图扇区中散射辐射的总和；
[0270]edir
＝∑e
dir
(θ,α)
[0271]edif
＝∑e
dif
(θ,α)
[0272]
其中，e
dir
(θ,α)和e
dif
(θ,α)分别为质心位于天顶角为θ和方位角为α的直接辐射量和散射辐射量；
[0273]
选取潭峪沟隧道地区dem图像上一点为观测点，所在位置海拔高度为980m，观测点位置及其半球视域计算结果分别如图5和图6所示；
[0274]
(2.1)太阳直接辐射计算：
[0275]edir
＝∑e
dir
(θ,α)
[0276]edir
(θ,α)＝s
const
*τ
bm(θ)
*sd
θ,α
*sung
θ,α
*cos(angi(θ,α))
[0277][0278]
angi(θ,α)＝arccos(cos(θ)*cos(gz)+sin(θ)*sin(gz)*cos(α-g
α
))
[0279]
其中，s
const
为太阳常数；τb为最短路径上的大气透过率；m(θ)为光学路径长度；sd
θ,α
为持续时间，sung
θ,α
为太阳图扇区的孔隙度；angi(θ,α)为天空扇区质心与表面垂线入射角度；gz和g
α
分别为表面的天顶角和方位角，θ1和θ2分别为天空扇区两边界处天顶角；h为海拔高程值；
[0280]
各个天空方向的太阳直接辐射可以通过太阳图计算得到，太阳图可表示太阳的运动轨迹；从半球视域角度来看，太阳图可分为若干个离散的扇区，每个扇区都有各自对应的天顶角和方位角，可认为扇区内直接辐射相同，计算每个扇区内的直接辐射并与半球视域叠加，可获得该观测点处的直接辐射分布，进而通过累加运算得到观测点处所接收到的总太阳直接辐射能量；
[0281]
计算从零点至下午三点该观测点处的太阳图，结合计算得到的半球视域可得到该观测点接收到的太阳直接辐射分布如图7所示；
[0282]
依次对每个位置计算其接收到的太阳直接辐射，可得到潭峪沟地区地表接收到的太阳直接辐射能量的分布，如图8所示；
[0283]
(2.2)太阳散射辐射计算
[0284]edif
＝∑e
dif
(θ,α)
[0285]edif
(θ,α)＝e
glb
*p*tc*skyg
θ,α
*w(θ,α)*cos(angi(θ,α))
[0286]eglb
＝(s
const
∑(τ
bm(θ)
))/(1-p)
[0287]
w(θ,α)＝(cosθ
2-cosθ1)/div
azi
[0288]
其中，s
const
为太阳常数；angi(θ,α)为天空扇区质心与表面垂线入射角度，e
glb
为正常总辐射量，p为散射辐射通量所占的比例，skyg
θ,α
为天空扇区的孔隙度，h为海拔高程值；θ1和θ2分别为天空扇区两边界处天顶角；div
azi
为天空图方位分割数量；w(θ,α)为天空扇区散射辐射量和总扇区散射辐射量的比值；tc为持续时间；τb为最短路径上的大气透过率；
[0289]
计算观测点位置处散射辐射时，需创建一个表示整个天空的半球视图，按照天顶角和方位角划分成一系列扇区，每个扇区内散射辐射可认为相同，然后根据各自方向计算每个扇区的散射辐射并与半球视域叠加，获得该观测点处的散射辐射分布，进而可得到观测点所接收到的总太阳散射辐射量；按照36个天顶分割(10
°
分割角)和50个方位分割，定义
[0311]
其中，r
nc
为植被冠层净辐射，r
ng
为土壤净辐射，q
cdsun
为植被冠层吸收的太阳直接辐射能量，q
cdsky
为被植冠层冠吸收的大气直接辐射能量，q
crsun
为植被冠层吸收的土壤反射太阳辐射能量，q
crsky
为植被冠层吸收的土壤反射大气辐射能量，η
l
为植被到土壤的能量传输效率，mc植被冠层吸收的热辐射能量，η
l
εc为土壤到植被的能量传输效率，mg土壤吸收的热辐射能量，q
gsun
为土壤吸收的太阳辐射能量，q
gsky
为土壤吸收的大气辐射能量，k
l
表示土壤表层的热导率，tg表示土壤表层辐射温度，t0表示土壤表层向下一定深度处的温度，ξz0为土壤表层厚度；
[0312]
在太阳辐射和大气辐射经过植被冠层到达土壤表层的过程中，由lambert-beer定律可知，植被叶片的阻挡会使辐射量出现衰减，通常采用消光系数描述衰减程度；对于均匀分布的植被群体，植被冠层对太阳辐射和大气辐射的消光系数可分别取0.4和0.8，太阳辐射和大气辐射的冠层透过率的计算方法为：
[0313]
ηs＝1-exp(-0.4lai)
[0314]
η
l
＝1-exp(-0.8lai)
[0315]
其中，lai为叶面积指数；叶面积指数采用的是随机数生成的方法进行计算，导致最后得到的地表植被红外辐射存在极大的误差；
[0316]
考虑到叶面积指数实际采样的难度和估算精度的问题，采用归一化植被指数(normalized difference vegetation index，ndvi)估算叶面积指数；首先对获取的多光谱图像预处理，然后利用多光谱图像可见光红色和近红外波段反射率按照下式计算归一化植被指数，最后采用非线性模型估算叶面积指数；
[0317][0318]
其中，ρ
nir
和ρr分别表示多光谱图像的近红外和红色波段的反射率；
[0319]
冠层吸收的直接太阳辐射和大气辐射以及冠层吸收的土壤反射的太阳辐射和大气辐射计算如下式：
[0320]qcdsun
＝(η
s-αc)q
sun
[0321]qcdsky
＝η
l
εcq
sky
[0322]qcrsun
＝(1-ηs)αgq
sun
[0323]qcrsky
＝(1-εg)(1-η
l
)εcq
sky
[0324]
其中，αc和εc分别为植被冠层的吸收率和发射率；αg和εg分别为土壤表层的吸收率和发射率；q
sun
为达到植被冠层的太阳辐射；q
sky
为达到植被冠层的大气辐射；
[0325][0326][0327]
其中，σ为斯蒂芬-玻尔兹曼常数；tc为待求解的植被冠层辐射温度，tg为土壤表层辐射温度；
[0328]
植被冠层与土壤以及空气之间的潜热交换量和显热交换量计算如下：
[0329][0330][0331][0332][0333][0334][0335]
其中，ρa为空气密度；c
p
为空气定压比热容；γ为干湿表常数；kd和da分别为空气温度为ta时的饱和比湿和空气比湿差；ra和rs分别为植被群体的边界层阻力和气孔阻力；根据相关气象数据和公式可得到上述参数的取值；es为太阳辐照度；ea为空气辐照度；
[0336]
通过联立求解土壤和植被冠层热平衡方程式组得到混合地表温度场分布，进而利用辐射公式得到混合地表红外辐射能量分布，如图12所示；
[0337]
(4)地下目标反演探测定位
[0338]
对山体背景热传播模型和太阳总辐射能量分布模型的分析计算，可获得山体表面植被发射的仿真红外辐射图像和山体表面接收的太阳总辐射能量分布图像；实际上，目标产生的扰动在传输过程中受到了山体介质的畸变作用，其热扰动为受介质畸变作用后的热扰动；这一畸变的扰动场叠加在山体背景热流场(由隧道沿线热能辐射模型，结合山体背景热传播模型和太阳总辐射能量分布模型，构建成的山体与空气层之间热辐射模型，山体与空气层之间热辐射模型输出即为山体背景热流场)上，形成了所观测到的红外辐射扰动图像；在探究山体中条带状地下隧道的扰动时，应该滤除山体背景热流场的干扰，并对畸变过程进行反变换，从而得到目标扰动的真实分布情况；
[0339]
伴随着山体内部的热量在岩层中的传导，会形成有方向的热流，类似于x射线具有的穿透性，热流可以从山体底部“穿透”岩层达到山体表面，从而形成矢量热流场；热流在不同岩层中传递时由于地质情况的差异会有不同程度的衰减，并且热流场与山体内部温度场具有一定的关联，可近似为线性关系；从红外遥感图像获得的山体表面红外辐射可以看作是山体各层热流场的叠加；
[0340]
理想情形下，单位面积的热流量如下：
[0341][0342]
其中，cm为岩层介质的容积热容量，k为岩层介质的热导率，为垂直方向的温度
梯度，可用差分近似计算；温度梯度需要利用温度场计算得到，可以看出二者具有近似的线性关系；
[0343]
位于同一海拔高度的地层温度场分布大致相同，同样地，根据岩层分布情况，可认为热流场也是相同的，计算每一层的热流场能量，利用通过红外遥感图像得到的山体表面真实红外辐射对不同厚度的山体热流场能量依次进行滤除，可以有效的保护目标信息，得到较为准确的目标扰动信号分布；
[0344]
由于山体中条带状地下隧道的传热规律符合热传导的数学模型，可根据传热学定理和规律，获得山体和目标的温度场分布，进而计算得到二者叠加的热流场分布；鉴于山体与空气层之间热辐射模型为三维结构，因此能量场的分布公式也应表示成三维形式，假设在某时刻t0山体表面(x0,y0,z0)处的辐射能量为bt(x0,y0,z0)，根据能量守恒定律，此辐射能量的值可通过地下条带状目标的热流场能量、山体背景热流场能量、它们与环境之间交换的辐射能量以及接收到的太阳辐射能量计算得到；
[0345]
bt(x0,y0,z0,t0)＝b(x,y,z,t)+t(x,y,z,t)+δ(x,y,z,t)+s(x,y,z,t)-bm(x,y,z,t)
[0346]
其中，bt(x0,y0,z0,t0)为山体表面辐射能量，b(x,y,z,t)为山体背景热流场能量，t(x,y,z,t)为条带状目标热流场能量，δ(x,y,z,t)为山体与空气接触面的辐射能量，s(x,y,z,t)为山体表面接收到的太阳辐射能量，bm(x,y,z,t)为条带状目标对应位置岩体的热流场能量；
[0347]
在真实环境中，山体与空气接触面的辐射能量很小，相比山体和目标的热流场能量，基本可以忽略；在本发明中，只考虑山体背景热流场能量、条带状目标热流场能量以及接收到的太阳辐射能量之间的关系；
[0348]
由上述公式可知，山体中条带状地下隧道产生的热流场能量、山体背景热流场能量以及太阳辐射能量叠加形成了可观测到的山体表面红外辐射能量分布；将上式加以变形，可以得到条带状目标热流场能量的计算公式如下：
[0349]
t(x,y,z,t)＝bt(x,y,z,t)+bm(x,y,z,t)-b(x,y,z,t)-s(x,y,z,t)
[0350]
在探测反演定位过程中，需要结合山体背景热流场能量分布，对不同厚度的山体背景所包含的能量进行计算，具体如下：
[0351][0352]
其中，a(z)和ε(z)分别为山体内部某一海拔高度处的横截面积和热流的衰减系数；z1(x,y)和z2(x,y)分别为山体表面高度坐标和从表面向下某一深度的高度坐标，通过调整z2(x,y)实现计算不同厚度的山体背景包含的热流场能量；
[0353]
潭峪沟隧道所在的山体海拔最高为1127m，为了获得更好的观测效果，从海拔高度为730m的地方开始反演探测定位，以20m的高度间隔进行分割，向下逐层剥离背景能量。通过不断向下分层，获得滤除20m、40m、60m、
…
深度的山体背景热流场能量后山体中条带状地下隧道扰动信号分布如图13(a)、图13(b)、图13(c)、图13(d)和图13(e)所示。
[0354]
反演地下隧道最优海拔高度的方法为：
[0355]
a.根据山体地表的实际总辐射量和红外遥感图像上的总辐射量，计算山体总辐射量误差；
[0356][0357]
其中,ε表示误差范围，是一个极小的正数eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；在红外遥感成像过程中，传感器和外界环境的影响会使测量的地表真实辐射量存在误差，因而利用红外遥感图像获得的山体表面辐射量是地表真实辐射量与误差项的叠加，如式：
[0358]
eti＝emi+σ
[0359]
其中，emi表示第i个像素点的山体地表的实际总辐射量，σ表示测量误差项。
[0360]
将测量误差纳入考虑的范围内后，收敛条件可表示为：
[0361][0362]
b.以隧道所在海拔高度ak为待反演变量，令k的初始值为1，初始化a1；
[0363]
c.采用当前红外遥感图像上各像素点总辐射量减去各像素点对应山体在海拔高度ak处的背景热流场能量求平方和，获取山体在海拔高度ak处的条带状地下隧道热流场能量；
[0364][0365]
其中，ak和n分别表示隧道海拔高度和计算区域的像素点数，eci表示第i个像素点对应的山体在海拔高度ak处的背景热流场能量，eti表示利用红外遥感图像上第i个像素点总辐射量；
[0366]
d.采用当前迭代的隧道所在海拔高度ak、当前迭代的搜索步长和搜索方向，确定下一次迭代的海拔高度a
k+1
；其中，当前迭代的搜索方向为利用上一次迭代的搜索方向、山体在海拔高度ak处的条带状地下隧道热流场能量和共轭参数计算获取；
[0367]ak+1
＝ak+λ
kdk
[0368]
其中，ak表示隧道海拔高度a的第k次运算值；λk表示搜索步长，dk表示搜索方向；
[0369][0370]
其中，表示山体在海拔高度ak处的条带状地下隧道热流场能量的梯度，可采用中心差分近似计算，βk表示共轭参数；
[0371]
e.将下一次迭代的海拔高度a
k+1
作为当前迭代的海拔高度，返回至步骤c，直至山体在海拔高度ak处的条带状地下隧道热流场能量与山体总辐射量误差相等，则结束迭代；
[0372]
f.根据每次迭代中的山体中条带状地下隧道热流场能量，构建扰动信号分布图像；且将海拔高度ak作为山体中条带地下隧道的最优化海拔高度。
[0373]
在探测反演定位过程中，当剥离背景能量至埋深180m时，隧道的扰动已经基本获取完整；剥离背景能量至200m的时候，可以看到此时扰动图像的变化已经很小了，与其上一个高度层的探测反演定位结果相比，扰动区域范围几乎没有发生变化；此时增加的部分主要是180～200m高度层内的区域干扰，经过上述实验证明，这一类干扰将在逐层的探测反演定位过程中被逐渐削弱。根据可见光图像对山体表面存在的岩石和房屋等干扰进一步滤
除，同时，由于计算地表接收的太阳辐射能量过程中所用的dem模型与真实地理位置相比存在一定的偏差，导致目标区域存在一部分负值点，对此进行保留，并将负值用该点在探测反演定位过程中最后一次出现的正值代替；最终，得到潭峪沟隧道区域的探测反演定位结果如图14所示。
[0374]
选取探测反演定位结果图像中扰动信号较为连续集中分布的区域，利用霍夫变换直线检测方法检测可能存在的隧道路段，处理山体中条带状地下隧道不连续的扰动信号，结果如图15所示；
[0375]
从结果图像中可以看出，受扰动点分布范围的影响检测出了多条直线，从中选取连贯性最好的直线进行连接，中间由于地势影响造成的扰动信号不明显的区域，将所选取的线段首尾连接近似表示，最终检测结果和隧道真实位置比较如图16所示。
[0376]
本发明通过反演探测定位的方法可以从红外遥感图像中获得较为准确的地下条带状目标埋深和位置信息。而不经过本发明，由于环境的强干扰因素以及目标信号的微弱性，很难实现对条带状目标的反演探测定位，通过此发明可以实现对山体环境中山体中条带状地下隧道的反演探测。
[0377]
本发明与现有技术相比，存在以下优势：
[0378]
本发明采用山体与空气层之间热辐射模型，结合dem数据，计算山体表面太阳辐射能量、迭代计算山体的各层山体背景热流场、隧道和山体背景热传播能量；充分考虑外部环境所带给山体的影响，如地表植被、太阳辐射以及地形因素带给红外数据的扰动差异；以本发明解决了山体中条带状地下隧道在红外图像中经地层传导调制后信号变弱，无法有效反演探测定位等问题，从而实现对山地环境中隧道的反演探测定位。
[0379]
本发明考虑了山体表面太阳辐射能量，在dem各个数据栅格上充分考虑太阳直接辐射和散射辐射的影响，相比于多时相红外图像分析目标信号随时间变化规律，建立数学模型计算求解隧道位置的方法相比，更加精确的考虑了外部环境(太阳直接辐射、散射辐射、地表植被等)对山体中条带状地下隧道能量扰动的影响，使得探测反演结果更加精确，相对于目标区域多时相红外数据的需求，此方法对数据获取更加容易；相比于对山地红外图像进行统计分析，寻找地下有隧道山体与旁边下无隧道山体的形成的特征模式，利用对可能含有隧道的山体环境的红外图像中进行一定方向的遍历搜寻该特征模式，从而实现隧道目标的探测、定位的方法，本发明无需大量选取山体中有无隧道的样本段，此方法效率更高；多源信息融合滤除环境干扰，使得山体中条带状地下隧道反演探测定位的精准度很高。
[0380]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。