CN103884431A - 平面地表环境中地下建筑的红外成像探测定位方法 - Google Patents

Abstract

本发明公开了一种平面地表环境中地下建筑的红外成像探测定位方法，包括：获取地下建筑经地层调制后所形成的原始红外图像g₀，并确定地下建筑在原始红外图像g₀中的大***置的局部红外图像g；设置迭代终止条件，并设定高斯热扩散函数的初始值h₀；以所述局部红外图像g作为初始目标图像f₀，根据所述高斯热扩散函数的初始值h₀，利用最大似然估计算法迭代求解热扩展函数h_n和目标图像f_n；判断是否满足迭代终止条件，如果满足，则本次迭代求解得到的目标图像f_n即为最终的目标图像f；若不满足，则继续迭代计算。本发明方法通过对地下建筑经地层调制后所形成的红外图像进行反调制处理，不仅使原来地下建筑的红外图像显示更清晰，还可以反演地下建筑的真实结构。

Description

平面地表环境中地下建筑的红外成像探测定位方法

技术领域

本发明属于地球物理学学科理论与遥感技术交叉的领域，更具体地，涉及一种平面地表环境中地下建筑的红外成像探测定位方法。

背景技术

目前，随着城市化的进程不断加快，越来越多的人口涌入城市，城市用地显得越来越紧张，这使得更多的建筑不得不选择修建在地下如：地下停车场，地下仓库，地下餐厅等；当然，也有很多设施是出于安全因素的考虑修建在地下，如：大型地下油库，军事设施等。这些都是地下建筑的典型例子，那么一旦这些地下建筑出现故障，就面临着勘察故障出现位置难度大的问题。当然不只是城市中的地下建筑存在勘察难的问题，考古、探矿、工程热物理、水坝探测等方面也面临着这样的问题。由此看来，对地下建筑的探测在民用上有着重要的意义。因此，有必要开展地下建筑探测识别的研究。

目前国内外主要是基于一种被动红外成像探测技术，其物理基础是大量的太阳光能量照射土壤被吸收产生热量，这些被加热的土壤发出红外辐射被热红外传感器探测。自然太阳能经过每日循环的加热和冷却对埋藏的物体和包围其周围的土地的影响是不同的，从而导致可探测的温差。地下建筑的存在可以产生地表水/热分布异常，可用于探测、发现、验证、确认地下建筑。

国内外使用红外成像技术探测地下建筑，但其均直接使用红外成像传感器进行成像，并对其得到的红外图像进行人工判读，这样就存在很大的局限性。其局限性在于：第一，地下建筑的热异常通过其埋入地层的传导调制，到达地表的热分布发生了与地下建筑有很大差异的变化，表现为热扩散严重，温差大大降低，热信号微弱。第二，通过地层的传导调制，地表热信号所反映的地下建筑的位置可能发生变化，难于发现与定位。第三，人工判读困难，不利于地下建筑的准确探测与定位。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提出了一种新的探测地下建筑的方法，通过对存在地下建筑的红外图像进行反调制处理，解决了地下建筑在红外图像中经地层传导调制后信号变弱，人工无法探测和定位的问题，从而对平面环境中地下建筑进行探测和定位。

本发明提供了一种平面地表环境中地下建筑的红外成像探测定位方法，所述方法根据地下建筑的能量扩散的高斯模型，对地下建筑经地层调制后所形成的原始红外图像进行反调制处理，得到地下建筑的目标图像，所述方法包括以下步骤：包括：

（1）获取地下建筑经地层调制后所形成的原始红外图像g₀，并确定地下建筑在原始红外图像g₀中的大***置的局部红外图像g；

（2）设置迭代终止条件，并设定高斯热扩散函数的初始值h₀；

（3）以所述局部红外图像g作为初始目标图像f₀，根据所述高斯热扩散函数的初始值h₀，利用最大似然估计算法迭代求解热扩展函数h_n和目标图像f_n；

（4）判断是否满足迭代终止条件，如果满足，则本次迭代求解得到的目标图像f_n即为最终的目标图像f；若不满足，则返回步骤（3），继续迭代计算。

优选地，所述步骤（1）具体为：

（1.1）将原始红外图像g₀分成N个m×m像素大小的红外图像区域；

（1.2）计算得到N个红外图像区域的平均灰度值，其中第i幅红外图像区域的平均灰度值为

g_ij为第i幅红外图像区域中的第j个像素的灰度值；

（1.3）分别计算出相邻红外图像区域的平均灰度差值ΔG＝G_i+1-G_i；

（1.4）将平均灰度值大且相邻平均灰度值差值ΔG小的区域作为需要反调制处理的局部红外图像g。

优选地，所述步骤（3）中迭代求解热扩展函数h_n和目标图像f_n具体根据下面两式迭代计算：

$h^{n+1}(x,y)=\left\{\right[\frac{g(i,j)}{h^n(x,y)*f^n(i,j)}]*f^n(-i,-j)\}*h^n(x,y)$

$f^{n+1}(i,j)=\left\{\right[\frac{g(i,j)}{h^{n+1}(x,y)*f^n(i,j)\}}]*h^{n+1}(-x,-y)\}*f^n(i,j)$

其中n为当前迭代次数,(x,y)属于热扩散函数支持域，(i,j)属于图像支持域，f(i,j)表示目标图像,g(i,j)表示局部红外图像,h(x,y)是热扩散函数，h(-x,-y)表示h(x,y)的共轭，f(-i,-j)表示f(i,j)的共轭，*为卷积运算符。

优选地，所述步骤（2）中的迭代终止条件为迭代终止次数n＞N₀或误差ε，所述步骤（4）中判断是否满足迭代条件具体为：

判断是否满足|g-hⁿ⁺¹*fⁿ⁺¹|＜ε或者n＞N₀，如果二者中任一个满足，则满足迭代终止条件，否则不满足。

本发明方法通过对地下建筑经地层调制后所形成的红外图像进行反调制处理后，不仅使原来地下建筑的红外图像显示更清晰，还可以反演地下建筑的真实结构。根据所获得的红外图像中已知的红外信息可以推断出地下建筑关键物理特征，达到“看见”并进一步定量测量出地下建筑关键物理特征的目的。

附图说明

图1为本发明流程示意图；

图2为地下建筑（深度100m）的地表热分布及分类结果划分；

图3为二维高斯函数形式示意图；

图4（a）为缩比模型在19:35时刻的红外热成像；

图4（b）为缩比模型在20:00时刻的红外热成像；

图4（c）为缩比模型在20:20时刻的红外热成像；

图4（d）为缩比模型在20:45时刻的红外热成像；

图4（e）为缩比模型在21:00时刻的红外热成像；

图4（f）为缩比模型在21:30时刻的红外热成像；

图5为最大似然估计图象恢复流程图；

图6为典型地下结构在平地，深度100m环境下，地下建筑地表热分布（右上）及其逆变换结果（右下）；

图7（a）为沙埋实验一得到的原始红外图像；

图7（b）为沙埋实验一反调制处理后的结果；

图8（a）为等效16米的地下建筑估计的热扩散函数，热扩散函数标准差为σ＝6.1；

图8（b）为等效32米的地下建筑估计的热扩散函数，热扩散函数标准差为σ＝9.3；

图9（a）为沙埋实验二原始红外图像；

图9（b）为沙埋实验二反调制处理后的结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明针对平面环境中存在地下建筑和不存在地下建筑的红外图像进行初步分类，然后对存在地下建筑的红外图像进行反调制处理，把探测问题归结为数学物理中的求逆问题。将地下建筑的结构信息和位置探测出来，实现对地下建筑目标的探测和定位，在现有的国内外文献中还没有看到与本发明相同的报道。

本发明所提供的一种平面地表环境中地下建筑的探测定位方法，包括下述几个步骤：

（1）平面环境中地下建筑红外特性分析与建模：以地球物理学、热学、遥感技术等学科的基本理论为指导，分析地下建筑与埋藏环境的物质/能量特征关系及相互作用规律，研究并建立目标区和背景区的红外信息映射模型，并对红外仿真的热图像进行分类学习。包括以下几个子步骤：

（1.1）本发明针对地下建筑红外特征及其与环境的关系先期开展了理论分析。利用现有的软件对典型的地下建筑进行建模仿真。为了将地下建筑经地层调制后的地表红外信息显示出来，利用ANSYS对地下建筑在平面地表环境下深度为100米时的地下建筑的热分布进行了仿真计算。地下建筑热辐射仿真是按照张朝晖编著的《ANSYS12.0热分析工程应用实战手册》中的热辐射分析实现。

（1.2）根据上述地下建筑的地表红外图像仿真结果进行聚类分析，所得到的红外图像中包含目标区红外信息和非目标区红外信息，利用近邻几何约束均值聚类算法把输入红外图像分为下有地下建筑红外图像和没有地下建筑的红外图像两类。

（1.2.1）选择红外图像中下有地下建筑的样本区和没有地下建筑的样本区，确定样本区的中心位置，即样本区的类心；

（1.2.2）计算每个样本区的平均灰度值，其中第i个红外图像区域的平均灰度值为：

$G_i=\frac{1}{m\×m}\underset{j=1}{\overset{m\×m}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}$

（1.2.3）计算出不同红外图像目标区的平均灰度差值ΔG＝G_i+1-G_i，由于非目标区之间的类心距离和目标区之间的类心距离较小，而非目标区和目标区的类间距离较大，且目标区的平均灰度值高于非目标去的平均灰度值。

（2）目标区域地层调制处理过程：一般的地下建筑由于各种要求必须保持恒温恒湿的状态且其状态一般不会发生改变，其本身存在的热辐射经过地层的调制，就会导致探测区域的物质、能量迁移状态发生改变，导致地表温度分布产生异常，并产生不同于其它探测区域的独特的红外信息场。其调制形式表现为一种高斯热扩散的调制形式。

我们假设拍摄环境中的地层是相对均匀的，拍摄的是在地下建筑热辐射稳定时的地表红外图像，我们做了试验说明在面积为150m²的单元，恒温25℃的条件下，在不同深度、不同地层介质时，通过能量交换在上覆地表体现出0.5K温差所需时间仿真计算结果如下表所示：

由上面的表格可以知道，不同介质的地层挡不住目标能量交换。同时，高斯热扩散函数标准差参数σ与地下建筑的埋深h和地层材质导热率ε可能存在有如下关系：

(h/ε)＝kσ+l

（3）目标区域红外图像的反调制处理方法：

本方法利用地表热晕图像探测识别地下建筑结构的试验，表明在对能量经地层正向调制变换高斯热扩散的基础上，进行相应的反变换处理，因为拍摄所得到的红外图像是地下建筑经地层调制后稳定的红外图像，我们可以选取一幅效果较好的稳定红外图像进行处理，处理后可以得到较好的目标图像，可揭示被地层遮蔽的地下建筑的结构信息。

（3.1）将两个10cm×10cm发热体分别埋于距正面10cm和20cm厚的沙土瓦砾子和1.5cm厚的大理石板中，拍摄距离为4.5m，烘鞋器发热体温度27℃，埋藏环境温度为22℃，且在表面覆盖一层草。利用中波红外窗口为3.4～5μm的中波红外相机进行拍摄。根据热物理学方程计算得知，该试验得到的目标对应地表与环境温差和10m*10m的地下建筑，埋深16m（混凝土2m，泥土14m）和埋深32m（混凝土2m，泥土30m）时目标上地表温度与环境温度差相当。

（3.2）目标区域的最大似然估计法的反卷积处理步骤：

（3.2.1）读入地下建筑经地层调制后所形成的原始红外图像g₀(i,j)，然后确定地下建筑在红外图像g₀(i,j)中的大***置的局部红外图像。我们利用分类的方法对该局部目标红外图像进行提取和处理，首先将原始红外图像g₀(i,j)中分成N个m×m像素大小的红外图像区域，计算得到N个红外图像区域的平均灰度值，其中第i幅红外图像区域的平均灰度值为：

$G_i=\frac{1}{m\×m}\underset{j=1}{\overset{m\×m}{\mathrm{\Σ}}}{g}_{\mathrm{ij}},i=1~N,$

g_ij为第i幅红外图像区域中的第j个像素的灰度值；然后分别计算出相邻红外图像区域的平均灰度差值ΔG＝G_i+1-G_i，由于相邻非目标区域的平均灰度差值ΔG比目标区域与非目标区域的平均灰度差值ΔG小，且相邻目标区域的平均灰度差值ΔG也比目标区域与非目标区域的平均灰度差值ΔG小，目标区域的平均灰度值也比非目标区域的平均灰度值大，故可以将平均灰度值大且相邻平均灰度值差值ΔG小的区域作为需要反调制处理的红外图像g(i,j)。

（3.2.2）地下建筑经地层调制后所形成的红外图像模型一般可表示为：

g(i,j)＝∫∫h(i,j;α,β)f(i-α,j-β)dαdβ+n(i,j)    (1)

式中g(i,j)为某一时刻的红外图象，f(i,j)为目标原图象，n(i,j)为传感器噪声，h(i,j;α,β)为高斯热扩散函数。

地层调制对地下建筑成像的影响通常可假定为线性移位不变，即模糊算子(热扩散函数)具有空间移不变性，可描述为：

h_t(i,j;α,β)＝h_t(α,β)   (i,j)∈Ω            (2)

将式(2)代入式(1)，则可得出如下卷积形式

g_t(i,j)＝∫∫_Dh_t(α,β)f(i-α,j-β)dαdβ+n(i,j)＝h_t(i,j)*f(i,j)+n(i,j)      (3)

由式(3)可知，在忽略噪声的情况下，只有得到高斯热扩散的热扩散函数，才能从经地层调制后的红外图象中恢复出地下建筑的结构。对近邻几何约束均值聚类后下有地下建筑的红外图象开展反调制处理，通过此图像来估计热扩散函数和目标图像。

（3.2.3）由于高斯热扩散函数是未知的，首先设置迭代终止次数N和误差ε，通过误差ε来检查是否满足反调制处理结果，并初始化最初的高斯热扩散模型参数σ；

（3.2.4）为了简化记号，我们对图像使用一维描述。首先定义目标强度为非负函数：{f(x),x∈X}，X为目标强度的支持域。定义地层对目标图像影响的热扩散函数为：{h(y|x),y∈Y}，Y为热扩散函数的支持域。为了减小误差，提高探测准确性，我们可以选取拍摄得到的一张稳定后的地表红外图像，可以更好的恢复地下建筑的结构。热扩散函数可看成是空间不变的，我们定义i(y)为红外图像对在坐标y处的强度，则

$i\left(y\right)=\underset{x\∈X}{\mathrm{\Σ}}h\left(y\right|x)f\left(x\right)=\underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)---\left(4\right)$

显然，由于h为空间不变的卷积算子，同时式（4）也可表示为

$i\left(y\right)=\underset{x\∈X}{\mathrm{\Σ}}h\left(x\right)f(y-x)y\∈Y---\left(5\right)$

在大多数实际情况下，i不可能被完美地检测到，它总是被一些噪声所污染，红外图象的某像元位置y处实际观测到的图像数据为g(y)。在许多情况下，地层调制后的红外图像具有Poisson分布性质。由此可见，在给定目标强度f和热扩散函数h条件下，可以假定g(y)是一个以i(y;f,h)为均值的服从Poisson分布的独立随机变量，因此，在象元位置y处取整数灰度值g(y)的概率可以表达为

$P\left(g\left(y\right)\right|f,h)=\frac{i^g{e}^{-i}}{g!}---\left(6\right)$

假定观测图象各像元是相互独立的，则其联合概率分布

$P\left(g(y_1{y}_2...y_n)\right|f,h)=\underset{y\∈Y}{\mathrm{\Π}}\frac{i^g{e}^{-i}}{g!}---\left(7\right)$

对(6)式取对数，得到其对数似然函数为

$-\begin{array}{c}\ln P\left(g(y_1{y}_2...y_n)\right|f,h)=\underset{y}{\mathrm{\Σ}}(-i+g\ln i-\ln g!)\\ =(-\underset{y}{\mathrm{\Σ}}i+\underset{y}{\mathrm{\Σ}}g\ln i)-\underset{y}{\mathrm{\Σ}}\ln g!\end{array}---\left(8\right)$

假定红外观测图象g(y₁y₂…y_n)在统计上是互相独立的，则泊松联合概率分布的对数似然函数为：

$\begin{array}{c}\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}(-i+g\ln i-\ln g!)\\ =-\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}i+\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}g\ln i-\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}\ln g!\end{array}---\left(9\right)$

可以看到式(8)中的最后一项为常数，因此它不影响似然函数的变化，为了简化问题可以将其舍去。由式(4)和(8)，我们可以得到对数似然函数为

$\begin{array}{c}L(f,\{h\left\}\right)=\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\\ +\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\left(g\ln \underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\right)\end{array}---\left(10\right)$

我们感兴趣的是从红外图像数据g₁(y₁y₂…y_n)中估计出目标强度f。当热扩散函数h(y|x)已知时，我们可以通过一些普通的图象去模糊处理技术来获得目标强度的估计。然而地下建筑的热扩散函数一般都是未知的，这就增加了去模糊的难度。

为了极大化对数似然函数，可将式(9)分别对各分量f(x)和h(x)求导并令其导数等于零，可以推导出

$\frac{\∂L(f,\{h\left\}\right)}{\∂f\left(x\right)}=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}h(y-x)+\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h(y-z)f\left(z\right)}=0---\left(11\right)$

由于高斯热扩散函数离散值之和为1，因此有

$\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h(y-z)f\left(z\right)}=1---\left(12\right)$

由此，我们可以建立如下迭代关系：

$f^{n+1}\left(x\right)=f^n\left(x\right)\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h^n(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}{h}^n(y-z)f^n\left(z\right)}---\left(13\right)$

对于本发明中的二维图像，式（13）应该为：

$f^{n+1}(i,j)=\left\{\right[\frac{g(i,j)}{h^{n+1}(x,y)*f^n(i,j)\}}]*h^{n+1}(-x,-y)\}*f^n(i,j)$

其中n为当前迭代次数,(x,y)属于热扩散函数支持域，(i,j)属于图像支持域，h(-x,-y)表示h(x,y)的共轭，f(-i,-j)表示f(i,j)的共轭。

为了便于对h(x)求导，将i用式（5）表示，代入式（9），此时，对数似然函数可等效地表示为

$\begin{array}{c}L(f,\{h\left\}\right)=\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\\ +\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\left(g\ln \underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\right)\end{array}---\left(14\right)$

对h(x)求导并令其为0，得

$\frac{\∂L(f,\{h\left\}\right)}{\∂h\left(x\right)}=-\underset{y\∈Y}{\mathrm{\Σ}}f(y-x)+\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h\left(z\right)f(y-z)}=0---\left(15\right)$

目标图象在观察图象的支持域之内，由于对目标原图像进行归一化处理，因此，目标原图像能量值之和为1，于是有

$\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h\left(z\right)f(y-z)}=1---\left(16\right)$

同理，当目标图象fⁿ⁺¹(x)用式（13）估计出来后，我们可以建立如下求解新高斯热扩散函数hⁿ⁺¹(x)的迭代关系：

$h^{n+1}\left(x\right)=h^n\left(x\right)\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f^{n+1}(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}{h}^n\left(z\right)f^{n+1}(y-z)}---\left(17\right)$

对于本发明中的二维图像，式（17）应该为：

$h^{n+1}(x,y)=\left\{\right[\frac{g(i,j)}{h^n(x,y)*f^n(i,j)}]*f^n(-i,-j)\}*h^n(x,y)$

其中n为当前迭代次数,(x,y)属于热扩散函数支持域，(i,j)属于图像支持域，h(-x,-y)表示h(x,y)的共轭，f(-i,-j)表示f(i,j)的共轭。

因此，我们可以通过最大似然估计法来恢复地下建筑的真实红外图像。

（3.2.5）利用最大似然估计的方法反复迭代红外图像g(x,y)，可以得到高斯热扩散函数h(x,y)和处理后的红外图像f(x,y)。

（4）地下建筑的恢复和定位：

存在地下建筑的红外图像经过反调制处理可以将高斯热扩散后的热能量收集起来，反调制处理后的红外图像就是收集能量之后的红外图像，该红外图像可以反映地下建筑的热分布情况，这样就尽可能真实的恢复了地下建筑结构信息，使人准确的对地下建筑进行定位。

本发明的特点在于：第一，本方法将高斯热扩散作为地下建筑经地层调制的形式；第二，本方法通过所做的缩比模型试验发现，对地下建筑和其周围环境的红外图像进行反调制处理后，可以很明显的看到地下建筑所处的位置，揭露地下建筑的本征信息，这样就可以准确的探测地下建筑；第三，本方法通过对仿真结果图反调制处理后，不仅使原来地下建筑的红外图像显示更清晰，而且将不可见的地下管状设施也表现出来了，同时可以反演地下建筑的真实结构。根据所获得的红外图像中已知的红外信息可以推断出地下建筑关键物理特征。达到“看见”并进一步定量测量出地下建筑关键物理特征的目的。

以下以实例具体地对本发明方法进行说明，本发明流程如图1所示，具体实施方法包括以下步骤。包括：平面环境中地下建筑红外特性分析与建模、目标区域地层调制处理过程、目标区域红外图像的反调制处理方法、地下建筑的恢复和定位。以下结合附图和具体实例对本发明进一步说明。

（1）平面环境中地下建筑红外特性分析与建模步骤：

本发明针对地下建筑红外特征及其与环境的关系先期开展了理论分析，为了将地下建筑经地层调制后的地表红外信息显示出来，利用ANSYS对地下建筑在平面地表环境下深度为100米时的热分布进行了仿真计算。地下建筑热辐射仿真是按照张朝晖编著的《ANSYS12.0热分析工程应用实战手册》中的热辐射分析实现的，具体的建模的步骤如下：

（1.1）选择表面效应单元SURF151，其利用实体表面的节点行程单元，并直接覆盖在实体单元的表面；其中，SURF151单元有主要实常数（角系数、Stefan-Boltzmann常数）、材料属性（密度、热辐射）、表面载荷（对流、热流密度）和体载荷（生热率）；

（1.2）设置表面效应单元对应的材料属性；设置密度为7800，比热为465，辐射率为1；

（1.3）设置表面效应单元的实常数；Stefan-Boltzmann常数为角系数的形状参数设置为1；

（1.4）创建几何模型、划分网格；Keypoint number设置为默认值，Global Element Sizes为3，即可生成有限元模型；

（1.5）利用Anasys加载求解；

（1.6）查看求解结果，并利用Anasys画出仿真图像，模型如图2所示。

（1.2）根据上述地下建筑的地表红外图像仿真结果进行聚类分析，所得到的红外图像中包含目标区红外信息和非目标区红外信息，利用近邻几何约束均值聚类算法把输入红外图像粗分为下有地下建筑红外图像和没有地下建筑的红外图像两类。

（1.2.1）选择红外图像中下有地下建筑的样本区和没有地下建筑的样本区，确定样本区的中心位置，即样本区的类心；本实例中，其中数字1代表不存在地下建筑的非目标区域，数字2代表存在地下建筑的非目标区。区域示意图如图3。

（1.2.2）计算每个样本区的平均灰度值，其中第i个红外图像区域的平均灰度值为：

$G_i=\frac{1}{m*m}\underset{j=1}{\overset{m*m}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}$

由于是计算机仿真结果，可以直接读出各个点的温度值，然后计算平均温度：

$T_i=\frac{1}{m*m}\underset{j=1}{\overset{m*m}{\mathrm{\Σ}}}{t}_{\mathrm{ij}}$

（1.2.3）计算出不同红外图像目标区的平均灰度差值ΔG＝G_i+1-G_i，由于非目标区之间的类心距离和目标区之间的类心距离较小，而非目标区和目标区的类间距离较大，且目标区的平均灰度值高于非目标区的平均灰度值。本实例中，区域2类内之间的平均温差为0.1K，区域1类内之间的平均温差为0K，区域1和区域2类间平均温差为0.3K。

（2）目标区域地层调制处理过程：

一般的地下建筑由于各种要求需要保持恒温恒湿的状态且其状态一般不会发生改变，其本身存在的热辐射经过地层的调制，就会导致探测区域的物质、能量迁移状态发生改变，导致地表温度分布产生异常，并产生不同于其它探测区域的独特的红外信息场。其调制形式表现为一种高斯热扩散的调制模型。我们通过试验来验证这一过程，试验具体步骤如下：

（2.1）将10cm×10cm发热体埋于距正面10cm厚沙土和1.5cm的大理石板中，拍摄距离4.5m，发热体温度27℃，埋藏环境的温度为22℃。通过对红外成像实验数据的时间变化情况分析发现，背景沙土区域的灰度变化缓慢，而目标所对应的沙土区域的灰度随着目标温度向表面的扩散而呈现较大的改变，图4（a）（b）（c）（d）（e）（f）是在不同时刻拍摄得到的红外图像。热辐射平衡之后拍摄得到的原始红外图像如图7（a）所示。

（2.2）图7（a）显示的是发热源经过地层的调制处理后形成的热晕图像，从该热晕红外图像中不容易准确定位发热源的准确位置，从热源的形状和相关的知识可以认为其调制形式为高斯热扩散的形式。地下建筑的热传导是空间上的热传导，由于只能探测到地下建筑向地表的热扩散信息，所以可以等效为二维高斯函数模型，模型如图3所示。

（3）目标区域红外图像的反调制处理方法：

地下建筑经地层高斯热扩散调制形成的地表热晕红外图像，经过步骤（1）可以知道地下建筑的大***置，但是地下建筑的具体结构并无法确定，为了得到地下建筑的具体结构信息，开展反调制处理，可揭示被地层遮蔽的地下建筑的结构信息。通过沙埋缩比试验对此进行验证，具体实施过程和结果如下：

（3.1）我们所做的沙埋缩比试验是模拟实场拍摄。拍摄的器材和背景环境如下表：

表1.1拍摄试验器材及相关数据

对试验精心设计后，我们通过多次试验，建立比照与参考以使得数据更加的准确和具有可比性。同时我们对图像进行了定标测温，得到整幅图像场景的温度。

（3.2）将两个10cm×10cm发热体分别埋于距正面10cm和20cm厚的沙土和1.5cm厚的大理石板中，拍摄距离为4.5m，发热体温度27℃，埋藏环境温度为22℃，且在表面覆盖一层草。拍摄得到的原始红外图像如图9（a）所示。根据热物理学方程计算得知，该试验得到的目标对应地表与环境温差和10m*10m的地下建筑，埋深16m（混凝土2m，泥土14m）和埋深32m（混凝土2m，泥土30m）时目标上地表温度与环境温度差相当。

（3.3）目标区域的最大似然估计法的反卷积处理步骤：

（3.3.1）读入地下建筑经地层调制后所形成的原始红外图像g₀(i,j)，然后确定地下建筑在红外图像g₀(i,j)中的大***置的局部红外图像。我们利用分类的方法对该局部目标红外图像进行提取和处理，首先将原始红外图像g₀(i,j)中分成N个m*m像素大小的红外图像区域，计算得到N个红外图像区域的平均灰度值，其中第i幅红外图像区域的平均灰度值为：

$G_i=\frac{1}{m*m}\underset{j=1}{\overset{m*m}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}$

然后分别计算出相邻红外图像区域的平均灰度差值ΔG＝G_i+1-G_i，由于相邻非目标区域的平均灰度差值ΔG比目标区域与非目标区域的平均灰度差值ΔG小，且相邻目标区域的平均灰度差值ΔG也比目标区域与非目标区域的平均灰度差值ΔG小，目标区域的平均灰度值也比非目标区域的平均灰度值大，故可以将平均灰度值大且相邻平均灰度值差值ΔG小的区域作为需要反调制处理的红外图像g(i,j)。

（3.3.2）地下建筑经地层调制后所形成的红外图像模型一般可表示为：

g(i,j)＝∫∫h(i,j;α,β)f(i-α,j-β)dαdβ+n(i,j)     (1)

式中g(i,j)为某一时刻的红外图象，f(i,j)为目标原图象，n(i,j)为传感器噪声，h(i,j;α,β)为高斯热扩散函数。

地层调制对地下建筑成像的影响通常可假定为线性移位不变，即模糊算子(热扩散函数)具有空间移不变性，可描述为：

h_t(i,j;α,β)＝h_t(α,β)  (i,j)∈Ω         (2)

将式(2)代入式(1)，则可得出如下卷积形式

g_t(i,j)＝∫∫_Dh_t(α,β)f(i-α,j-β)dαdβ+n(i,j)＝h_t(i,j)*f(i,j)+n(i,j)    (3)

由式(3)可知，在忽略噪声的情况下，只有得到高斯热扩散的热扩散函数，才能从经地层调制后的红外图象中恢复出地下建筑的结构。对近邻几何约束均值聚类后下有地下建筑的红外图象开展反调制处理，通过此图像来估计热扩散函数和目标图像。原始红外图像如图9（a）。

（3.3.3）设置迭代终止次数N和误差ε，通过误差ε来检查是否满足反调制处理结果，并设置最初的高斯热扩散模型参数σ；本实例中，N=40，ε=0.003，σ=4.0。

（3.3.4）定义目标强度为非负函数：{f(x),x∈X}，其中X为目标强度的支持域，定义地层对目标图像影响的热扩散函数为：{h(y|x),y∈Y}，Y为观察图象的支持域。为了减小误差，提高探测准确性，我们可以选取拍摄得到的一张稳定后的地表红外图像，可以更好的恢复地下建筑的结构。热扩散函数可看成是空间不变的，我们定义i(y)为红外图像对在坐标y处的强度，则

$i\left(y\right)=\underset{x\∈X}{\mathrm{\Σ}}h\left(y\right|x)f\left(x\right)=\underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)---\left(4\right)$

显然，由于h为空间不变的卷积算子，同时式（4）也可表示为

$i\left(y\right)=\underset{x\∈X}{\mathrm{\Σ}}h\left(x\right)f(y-x)y\∈Y---\left(5\right)$

在给定目标强度f和热扩散函数h条件下，可以假定g(y)是一个以i(y;f,h)为均值的服从Poisson分布的独立随机变量，因此，在象元位置y处取整数灰度值g(y)的概率可以表达为

$P\left(g\left(y\right)\right|f,h)=\frac{i^g{e}^{-i}}{g!}---\left(6\right)$

假定观测图象各像元是相互独立的，则其联合概率分布

$P\left(g(y_1{y}_2...y_n)\right|f,h)=\underset{y\∈Y}{\mathrm{\Π}}\frac{i^g{e}^{-i}}{g!}---\left(7\right)$

对(6)式取对数，得到其对数似然函数为

$-\begin{array}{c}\ln P\left(g(y_1{y}_2...y_n)\right|f,h)=\underset{y}{\mathrm{\Σ}}(-i+g\ln i-\ln g!)\\ =(-\underset{y}{\mathrm{\Σ}}i+\underset{y}{\mathrm{\Σ}}g\ln i)-\underset{y}{\mathrm{\Σ}}\ln g!\end{array}---\left(8\right)$

假定红外观测图象g(y₁y₂…y_n)在统计上是互相独立的，则泊松联合概率分布的对数似然函数为：

$\begin{array}{c}\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}(-i+g\ln i-\ln g!)\\ =-\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}i+\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}g\ln i-\mathrm{\Σ}\underset{y}{\mathrm{\Σ}}\ln g!\end{array}---\left(9\right)$

可以看到式(8)中的最后一项为常数，因此它不影响似然函数的变化，为了简化问题可以将其舍去。由式(4)和(8)，我们可以得到对数似然函数

$\begin{array}{c}L(f,\{h\left\}\right)=\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\\ +\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\left(g\ln \underset{x\∈X}{\mathrm{\Σ}}h(y-x)f\left(x\right)\right)\end{array}---\left(10\right)$

我们感兴趣的是从红外图像数据g₁(y₁y₂…y_n)中估计出目标强度f。当热扩散函数h(y|x)已知时，我们可以通过一些普通的图象去模糊处理技术来获得目标强度的估计。然而地下建筑的热扩散函数一般都是未知的，这就增加了去模糊的难度。

为了极大化对数似然函数，可将式(9)分别对各分量f(x)和h(x)求导并令其导数等于零，可以推导出

$\frac{\∂L(f,\{h\left\}\right)}{\∂f\left(x\right)}=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}h(y-x)+\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h(y-z)f\left(z\right)}=0---\left(11\right)$

由于高斯热扩散函数离散值之和为1，因此有

$\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h(y-z)f\left(z\right)}=1---\left(12\right)$

由此，我们可以建立如下迭代关系：

$f^{n+1}\left(x\right)=f^n\left(x\right)\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}g\frac{h^n(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}{h}^n(y-z)f^n\left(z\right)}---\left(13\right)$

为了便于对h(x)求导，将i用式（5）表示，代入式（9），此时，对数似然函数可等效地表示为

$\begin{array}{c}L(f,\{h\left\}\right)=\ln P\left(\right\{g\left\}\right|f,\left\{h\right\})=-\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\underset{x\∈X}{\mathrm{\Σ}}h\left(x\right)f(y-x)\\ +\mathrm{\Σ}\underset{y\∈Y}{\mathrm{\Σ}}\left(g\ln \underset{x\∈X}{\mathrm{\Σ}}h\left(x\right)f(y-x)\right)\end{array}---\left(14\right)$

对h(x)求导并令其为0，得

$\frac{\∂L(f,\{h\left\}\right)}{\∂h\left(x\right)}=-\underset{y\∈Y}{\mathrm{\Σ}}f(y-x)+\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h\left(z\right)f(y-z)}=0---\left(15\right)$

目标图象在观察图象的支持域之内，由于对目标原图像进行归一化处理，因此，目标原图像能量值之和为1，于是有

$\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}h\left(z\right)f(y-z)}=1$

同理，当目标图象fⁿ⁺¹(x)用式（13）估计出来后，我们可以建立如下求解新热扩散函数hⁿ⁺¹(x)的迭代关系：

$h^{n+1}\left(x\right)=h^n\left(x\right)\underset{y\∈Y}{\mathrm{\Σ}}g\frac{f^{n+1}(y-x)}{\underset{z\∈X}{\mathrm{\Σ}}{h}^n\left(z\right)f^{n+1}(y-z)}---\left(16\right)$

上述实施例中，我们利用的是二维高斯函数，由于高斯函数具有可分离性，二维函数的卷积可以用两次一维函数卷积来代替，所以反调制的过程是可以使用一维来描述的，针对具体的二维图像，将公式对应的变换为二维形式即可。

因此，我们可以通过使用最大似然估计法来设法找到红外目标图像和高斯热扩散函数。最大似然估计图像恢复流程图如图5所示。

（3.3.5）利用最大似然估计的方法反复迭代红外图像g(x,y)，可以得到高斯热扩散函数和处理后的红外目标图像，从处理后的红外目标图像中可以清楚地看到地下建筑的具体结构。图6是计算机仿真图像反调制处理前后效果的对比图，同时对实际地下建筑和反调制后的红外图像进行了标记。从图6可以看出反调制处理得到的结果图可以更清晰的探测和定位目标，提高了人工判读的准确性，同时可以将看不到地下建筑的结构反演出来。利用所做的沙埋试验对本发明进行验证，得到图9（a）的等效16米的地下建筑的高斯热扩散函数如图8（a）所示，等效32米的地下建筑的高斯热扩散函数如图8（b）所示，其中σ为获得的高斯热扩散函数标准差。图7（b）是图7（a）经过反调制处理之后所得到的结果图。

(4)地下建筑的恢复和定位：

存在地下建筑的红外图像经过反调制处理可以将高斯热扩散后的热能量收集起来，在反调制处理后的红外图像中可以看到收集能量之后的红外图像，该红外图像可以反映地下建筑的热分布情况，这样就尽可能真实的恢复了地下建筑结构信息，使人准确的对地下建筑进行定位。图9（b）是图9（a）经过反调制处理之后所得到的结果图。对比图9（a）和图9（b），从图9（a）中根本无法看出发热源的位置，而从图9（b）中可以清楚地看到发热源所处的位置。实现了地下建筑的定位。同时从图6可以看出反调制处理得到的结果图可以更清晰的探测和定位目标，提高了人工判读的准确性，同时可以将看不到地下建筑的结构反演出来。

Claims (4)

1.一种平面地表环境中地下建筑的红外成像探测定位方法，其特征在于，所述方法根据地下建筑的能量扩散的高斯模型，对地下建筑经地层调制后所形成的原始红外图像进行反调制处理，得到地下建筑的目标图像，所述方法包括以下步骤：

（1）获取地下建筑经地层调制后所形成的原始红外图像g₀，并确定地下建筑在原始红外图像g₀中的大***置的局部红外图像g；

（2）设置迭代终止条件，并设定高斯热扩散函数的初始值h₀；

（3）以所述局部红外图像g作为初始目标图像f₀，根据所述高斯热扩散函数的初始值h₀，利用最大似然估计算法迭代求解热扩展函数h_n和目标图像f_n；

（4）判断是否满足迭代终止条件，如果满足，则本次迭代求解得到的目标图像f_n即为最终的目标图像f；若不满足，则返回步骤（3），继续迭代计算。

2.如权利要求1所述的方法，其特征在于，所述步骤（1）具体为：

（1.1）将原始红外图像g₀分成N个m×m像素大小的红外图像区域；

（1.2）计算得到N个红外图像区域的平均灰度值，其中第i幅红外图像区域的平均灰度值为

g_ij为第i幅红外图像区域中的第j个像素的灰度值；

（1.3）分别计算出相邻红外图像区域的平均灰度差值ΔG＝G_i+1-G_i；

（1.4）将平均灰度值大且相邻平均灰度值差值ΔG小的区域作为需要反调制处理的局部红外图像g。

3.如权利要求1或2所述的方法，其特征在于，所述步骤（3）中迭代求解热扩展函数h_n和目标图像f_n具体根据下面两式迭代计算：

$h^{n+1}(x,y)=\left\{\right[\frac{g(i,j)}{h^n(x,y)*f^n(i,j)}]*f^n(-i,-j)\}*h^n(x,y)$

$f^{n+1}(i,j)=\left\{\right[\frac{g(i,j)}{h^{n+1}(x,y)*f^n(i,j)\}}]*h^{n+1}(-x,-y)\}*f^n(i,j)$

其中n为当前迭代次数,(x,y)属于热扩散函数支持域，(i,j)属于图像支持域，f(i,j)表示目标图像,g(i,j)表示局部红外图像,h(x,y)是热扩散函数，h(-x,-y)表示h(x,y)的共轭，f(-i,-j)表示f(i,j)的共轭，*为卷积运算符。

4.如权利要求3所述的方法，其特征在于，所述步骤（2）中的迭代终止条件为迭代终止次数n＞N₀或误差ε，所述步骤（4）中判断是否满足迭代条件具体为：

判断是否满足|g-hⁿ⁺¹*fⁿ⁺¹|＜ε或者n＞N₀，如果二者中任一个满足，则满足迭代终止条件，否则不满足。

Images