CN104637073A - 一种基于太阳照射阴影补偿的带状地下结构探测方法 - Google Patents

Abstract

一种基于太阳照射阴影补偿的带状地下结构探测方法，属于遥感技术、自然地理和模式识别的交叉领域，对阴影进行检测后，进行补偿处理，提高探测带状地下结构的识别率，降低虚警率。本发明包括获取指定区域的DEM地形数据步骤、利用DEM和太阳高度角、太阳方位角获取图像阴影位置步骤、对阴影区域进行处理和补偿步骤、阴影区域校正后带状地下结构的探测步骤。本发明利用获取的DEM地形数据对指定区域的阴影进行检测；将检测出来的阴影区域进行处理和补偿操作，减小阴影区域对带下地下结构探测的影响；最后对阴影补偿后的遥感图像进行带状地下结构探测，相对未进行阴影补偿处理的遥感图的带状地下结构探测，提高了带状地下结构探测正确率，降低了虚警率。

Description

一种基于太阳照射阴影补偿的带状地下结构探测方法

技术领域

本发明属于遥感技术、图像处理和模式识别的交叉领域，具体涉及一种基于太阳照射阴影补偿的带状地下结构探测方法。

背景技术

阴影是在自然界中普遍存在的一种物理现象，是由于光源被物体遮挡而产生的。图像中阴影的存在对计算机视觉领域的相关问题有不同的影响，这种影响有有利的也有不利的，比如在虚拟现实、3D游戏中为物体添加阴影，可以提高场景的真实感。但是更多的时候，图像中的阴影会对计算机视觉的相关问题产生不利的影响，如在遥感图像中，阴影的存在会影响后继的图像匹配、模式识别和地物的提取等多种遥感图像处理操作；在视频监控中，阴影和运动目标结合在一起，导致计算机对目标物体的提取和追踪出现错误。因此，有必要对图像中的阴影进行检测和分析，并根据需要进而消除或减弱阴影的影响。

目前，国内外很多专家学者对视频中的阴影去除进行了比较深入的研究，提出了很多有效的算法。在视频监控和车辆追踪等领域得到了广泛的应用，但是无法处理静态的图像。在静态图像中的阴影检测方法，一直是一个比较困难的问题，因为静态图像所包含的信息量更少，不能像在视频中一样利用帧间关系来进行阴影检测。同时，在遥感成像等基于静态图像的研究领域，在检测出阴影之后，还要尽可能的减弱或消除阴影的影响，恢复阴影中物体的本来面目，从而提高红外遥感图像的真实性。

目前国内外并未有使用红外成像技术探测带状地下结构的相关报道，比如暗河，地下隧道等带状地下结构。红外遥感图像中阴影的存在增加了对带状地下结构探测的虚警率，降低了探测率。

发明内容

本发明提出了一种基于太阳照射阴影补偿的带状地下结构探测方法，结合数字高度模型(Digital Elevation Model，DEM)对单幅红外遥感图像进行阴影检测，并对阴影所在的区域的灰度进行了补偿，从而减弱了阴影对探测的影响。

本发明提供的一种基于太阳照射阴影补偿的带状地下结构探测方法，具体步骤如下：

(1)根据区域的地理位置得到区域的经纬度信息，然后通过图像的分辨率计算经纬度步长，获取指定区域的海拔高度信息，最后对所得到的海拔高度信息进行校正处理，具体过程包括以下子步骤：

(1.1)利用所检测的区域确定经纬度步骤

(1.1.1)经纬度信息的确定可以根据需要检测的位置和所需数据的分辨率来确定，待检测区域左上方点经纬度信息为A(x₁，y₁)，右下方点的经纬度信息为B(x₂，y₂)。

(1.1.2)根据仿真的红外遥感图像的设置参数可知红外遥感图像的分辨率为k米，因此需要将DEM地形数据的分辨率也设置为k米。地球上任意两点的实际距离D满足一定的关系，关系式如下：

C＝sin(y₁)*sin(y₂)*cos(x₁-x₂)+cos(y₁)*cos(y₂)

D＝R*arccos(C)*Pi/180

R＝6371.004

上式中(x₁，y₁)和(x₂，y₂)分别代表这两点的经纬度信息。

但是，我们只需求平行于赤道和垂直于赤道方向的两点距离，故平行于赤道时的距离d₁和经纬度的关系如下：

C＝sin(y₁)²*cos(x₁-x₂)+cos(y₁)² (y₁＝y₂)

d₁＝R*arccos(C)*Pi/180

R＝6371.004

lon＝k*(x₂-x₁)/d₁

式中，d₁为平行于赤道的两点距离，lon为经度步长。

垂直于赤道时的距离d₂和经纬度的关系如下：

C＝cos(y₁-y₂) (x₁＝x₂)

d₂＝R*arccos(C)*Pi/180

R＝6371.004

lat＝k*(y₁-y₂)/d₂

式中，d₂为垂直于赤道的两点距离，lat为纬度步长。

根据以上的经度步长和纬度步长得到待检测区域每一点的经纬度信息，其经纬度分辨率为k米。

(1.2)利用Google Earth获取所需区域的地形数据步骤

根据所获取的经纬度信息，利用Google Earth提供的编程接口进行编程，首先根据每一点的经纬度信息在Google Earth中进行定位，然后再定位成功后获取定位点的海拔信息。最后就会得到一张和经纬度图一样大小的海拔信息图，其分辨率同样为k米。

(1.3)对所获取的地形数据进行校正处理步骤

(1.3.1)通过对所获取的海拔高度信息观察发现，其中有许多的奇异点，有些位置的海拔高度明显高于邻近点的海拔高度，通过下式对以上的海拔信息图处理可以求出海拔高度中的奇异点，其中奇异点和非奇异点满足以下条件：

H(i,j)>H(i+k,j+l)+h (k,l＝-1,1)式中(i，j)为奇异点

H(i,j)<H(i+k,j+l)+h (k,l＝-1,1)式中(i，j)为非奇异点

式中H(i，j)为(i，j)点处的海拔高度，H(i+k，j+l)为(i，j)邻近点的海拔高度。

由于中值滤波对奇异点的处理效果很好，故对奇异点进行中值滤波处理，可以去除掉奇异点，得到初步校正后海拔高度图。

(1.3.2)对初步校正后的海拔高度图再次校正，利用m*m大小的模板对海拔高度图进行均值滤波处理。

进行以上的步骤是由于所获取到的区域海拔数据并不是准确无误的，海拔高度信息会存在一定的误差，会给以后的检测带来不利的影响，故需要尽量使海拔高度信息准确。由于误差的出现具有随机性，有时海拔高度比实际的高，有时海拔高度比实际的低，所以利用其随机性消除这种误差，其原理主要是基于以下的公式：

$\begin{array}{c}\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}(H_i+{\mathrm{\&epsiv;}}_i)\\ =\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{H}_i+\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_i\\ =\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{H}_i\end{array}$

式中H_i+ε_i为读取的海拔高度，H_i为实际的海拔高度，ε_i为高度误差。从上式中可以看出，由于误差的随机性，平均误差最后得到比较准确的区域海拔高度信息数据。

(2)利用太阳方位角等图像信息对校正后的海拔高度图进行旋转处理，通过太阳高度角信息计算山体本体和山体投影的比例，确定阴影的大小和位置，再经过一系列的处理得基于DEM的图像阴影位置，具体过程主要包括以下子步骤：

(2.1)利用太阳方位角旋转海拔高度图步骤

利用产生仿真图像时设置的参数确定太阳方位角θ，然后对校正后的海拔高度图旋转太阳方位角θ，旋转可以看做是从一个坐标系变换到了另一个坐标系的仿射变换，为了减小旋转带来的误差，采用双线性插值的方法对图像进行旋转处理。仿射变换的计算公式如下：

最后得到旋转后的海拔高度图，超出的部分海拔高度设置为0。

上述中太阳方位角为竖立在地面上的直线的阴影和正南方的夹角。

(2.2)计算山体本体和山体投影的比例，确定阴影的大小和位置

其中阴影位置和太阳高度角也有关系，因为太阳高度角决定山体本体和山体投影的比例，这个比例直接决定阴影的大小，而太阳高度角的计算公式为：

φ代表当地的地理纬度，代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+。

山体本体和山体阴影的关系为：

由上式可以得到投影的范围大小，从而可以确定具体阴影的大小和位置。

(2.3)阴影的区域的检测标记步骤

(2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：

$H(:,i)\&DoubleLeftRightArrow;H(:,n-i)(i=\mathrm{1,2,3}...\frac{2}{n})$

式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，代表交换该符号前后的数值。

将旋转后得到的图进行左右倒置处理，是为了保证小的山丘在大的山丘的阴影部位的时候不被检测为非阴影区域，从而相当于太阳是从东边照射过来的。

(2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，为了对其中的每一个点进行判断，判断该点是否存在于阴影区，判断的公式如下：

H'(i,j)≠0(i＝1,2,3...m、j＝1,2,3...n)

$\begin{array}{cc}(k-j)*10<(H^{\&prime;}(i,k)-H^{\&prime;}(i,j))*\frac{1}{\tan \mathrm{\&beta;}}& (k=j+1,j+2...n)\end{array}$

式中H’(i，j)代表旋转倒置后的海拔高度图，m和n分别代表海拔高度图的行数和列数。

如果满足以上等式的条件，就标记该点为阴影区域。

(2.3.3)将处理后的标记图左右倒置回来，将倒置回来的标记图进行旋转，旋转的角度就是360-太阳方位角θ。将旋转后得到的标记图大小和山体的红外遥感图像F(x，y)对应起来，就可以根据标记图对山体红外遥感图像中的阴影部分进行标记处理。

(3)通过对得到的阴影部分的标记图进行膨胀处理，对膨胀后阴影区域进行一定的补偿，具体过程主要包括以下子步骤：

(3.1)阴影标记图的膨胀处理步骤

对得到的标记图进行膨胀处理以得到连续的区域。膨胀处理就是用一个结构元素在标记图上移动并进行比较，根据比较的结果做处理，处理的过程为：

a、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点没有一个相同，该点就为黑色，即非阴影区域。

b、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点有一个相同，该点就为白色，即阴影区域。

经过上述的处理过程就可以得到膨胀后的标记图，将不连续的区域变得连续，更加符合实际。

(3.2)对阴影区域进行补偿处理步骤

假定图像是局部平稳的，可以认为阴影区域与其周围一定范围内的非阴影区域的地质结构同质，红外特性相似。在检测出阴影区域之后，对各个独立阴影区域与其邻近的非阴影区域进行匹配，完成阴影区域补偿操作，从而达到去除阴影的目的。算法过程如下：

1)对于阴影区域，通过下式计算其邻近的非阴影区域

Ω_noshadow＝{p|0<d(p,Ω_shadow)<dist}

式中Ω_noshadow表示邻近某个距离阈值dist的非阴影区域集合，

2)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：

${F^{\&prime;}}_{\mathrm{shadow}}(i,j)=A[m_{\mathrm{noshadow}}+\frac{F_{\mathrm{shadow}}(i,j)-m_{\mathrm{shadow}}}{{\mathrm{\&sigma;}}_{\mathrm{shadow}}}*{\mathrm{\&sigma;}}_{\mathrm{noshadow}}]$

D＝F'_shadow(i,j)-F_shadow(i,j)

式中F_shadow(i,j)是补偿之前的阴影区域灰度值，F'_shadow(i,j)是补偿之后阴影区域灰度值，m_shadow和σ_shadow是阴影区域灰度值的均值和方差，m_noshadow和σ_noshadow是邻近非阴影区域灰度值的均值和方差，A为补偿强度系数，D为阴影区域需要的补偿值。

从上述中得到需要补偿的灰度值D后，遍历所有检测到的阴影区域，将阴影区域对应到红外遥感图像上每一个点的灰度值补偿D，得到最终补偿后的红外遥感图像。

(4)对阴影区域补偿后的红外遥感图像进行初步的探测识别，然后利用位置的关联性进行二次的判断，提高识别率，降低虚警率，具体的过程包含以下子步骤：

(4.1)初步判断带状地下结构位置步骤

利用阴影补偿后的红外遥感图像探测带状地下结构，依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别，每段的长度为250米宽度为30米(即25个像素点长3个像素点宽)，对像素点的左右遍历，探测是否存在带状地下结构，最后得到初步的隧道识别结果。同时在垂直于隧道的方向利用相同的方法探测地下结构，计算虚警率。

(4.2)根据位置关联性进行二次判断步骤

对初步识别出的符合脉冲模式的分散段进行统计分析，依据500米长度(50个像素点)上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白。

1)若根据判断可以忽视掉这些分散脉冲模式段中间的小空白，那么该500的长度段就确认为带状目标所在位置。

2)若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白，那么该500的长度段就确认为非带状目标所在位置。

与现有技术相比，本发明具有以下有益效果：

1、利用获取的指定区域的DEM地形数据和太阳方位角、太阳高度角等有用信息，可以有效的将静态红外图像中太阳照射的阴影区域检测出来，并利用阴影区域与其周围一定范围内的非阴影区域的地质结构同质，红外特性相似等知识，对阴影区域进行了一定的补偿处理，这样可以减小阴影对地下结构探测的不利影响。

2、对太阳照射的阴影补偿后的结果图进行地下结构探测，相对于对未进行阴影补偿的红外图进行地下结构探测来说，能增加地下结构的探测正确率，同时降低了虚警率。

附图说明

图1为本发明基于太阳照射阴影补偿的带状地下结构探测方法流程示意图；

图2为本发明实施例中经纬度和距离的示意图；

图3为本发明实施例中未处理奇异点前的海拔高度图；

图4为本发明实施例中处理奇异点后海拔高度图；

图5为本发明实施例中海拔高度旋转过程示意图；

图6为本发明实施例中太阳方位角示意图；

图7为本发明实施例中太阳高度角示意图；

图8(a)为本发明实施例中仿真隧道1区域的阴影标记图；

图8(b)为本发明实施例中仿真隧道1区域的阴影标记膨胀图；

图9(a)为本发明实施例中仿真隧道2区域的阴影标记图；

图9(b)为本发明实施例中仿真隧道2区域的阴影标记膨胀图；

图10(a)为本发明实施例中仿真隧道1区域的红外遥感图像；

图10(b)为本发明实施例中仿真隧道1区域的阴影补偿红外遥感图像；

图11(a)为本发明实施例中仿真隧道2区域的红外遥感图像；

图11(b)为本发明实施例中仿真隧道2区域的阴影补偿红外遥感图像；

图12(a)为本发明实施例中未补偿除阴影的仿真隧道1的检测结果图；

图12(b)为本发明实施例中补偿阴影后的仿真隧道1的检测结果图。

具体实施方式

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

本发明的流程如图1所示，其中具体的实施方法包括以下步骤。包括：获取指定区域的DEM地形数据步骤、利用DEM和太阳高度角、太阳方位角获取图像阴影位置步骤、对阴影区域进行处理和补偿步骤、阴影区域校正后带状地下结构的探测步骤：

(1)根据区域的地理位置得到区域的经纬度信息，然后通过图像的分辨率计算经纬度步长，获取指定区域的海拔高度信息，最后对所得到的海拔高度信息进行校正处理，具体过程包括以下子步骤：

(1.1)确定指定区域的经纬度步骤

(1.1.1)经纬度信息的确定可以根据需要检测的位置和所需数据的分辨率来确定。如图2所示，本实例中，待检测区域左上方点经纬度信息为A(x₁，y₁)＝(116.150049，40.296833)，右下方点的经纬度信息为

B(x₂，y₂)＝(116.194775，40.260787)。

(1.1.2)所仿真的红外遥感图像的分辨率为k米，需要将DEM地形数据的分辨率也设置为k米，地球上任意两点的实际距离d满足一定的关系，关系式如下：

C＝sin(y₁)*sin(y₂)*cos(x₁-x₂)+cos(y₁)*cos(y₂)

D＝R*arccos(C)*Pi/180

R＝6371.004

上式中(x₁，y₁)和(x₂，y₂)分别代表这两点的经纬度信息。

如图2所示，故平行于赤道时的距离d₁和经纬度的关系如下：

C＝sin(y₁)²*cos(x₁-x₂)+cos(y₁)² (y₁＝y₂)

d₁＝R*arccos(C)*Pi/180

R＝6371.004

lon＝k*(x₂-x₁)/d₁

本实例中，d₁＝3800m，lon＝0.0001177。

垂直于赤道时的距离d₂和经纬度的关系如下：

C＝cos(y₁-y₂) (x₁＝x₂)

d₂＝R*arccos(C)*Pi/180

R＝6371.004

lat＝k*(y₁-y₂)/d₂

本实例中，d₂＝4000m，lat＝0.000090115。

(1.2)利用Google Earth获取所需区域的地形数据步骤

根据所获取的经纬度信息，利用Google Earth提供的编程接口进行编程，首先根据每一点的经纬度信息在Google Earth中进行定位，然后再定位成功后获取定位点的海拔信息。最后就会得到一张和经纬度图一样大小的海拔信息图，其分辨率同样为k米。

(1.3)地形数据的矫正处理步骤

(1.3.1)如图3所示，通过对图3的海拔高度信息观察发现，其中有许多的奇异点，有些位置的海拔高度明显高于邻近点的海拔高度，为了减少误差，通过和邻近的比较来确定是否为奇异点，其中奇异点和非奇异点满足以下条件：

H(i,j)>H(i+k,j+l)+h (k,l＝-1,1)式中(i，j)为奇异点

H(i,j)<H(i+k,j+l)+h (k,l＝-1,1)式中(i，j)为非奇异点

式中H(i，j)为(i，j)点处的海拔高度，H(i+k，j+l)为(i，j)邻近点的海拔高度。本实例中，h为20。

由于中值滤波对奇异点的处理效果很好，故对奇异点进行中值滤波处理，可以去除掉奇异点，得到初步校正后海拔高度图的。结果如图4所示。

(1.3.2)对初步校正后的海拔高度图再次校正，利用m*m大小的模板对海拔高度图进行均值滤波处理。

进行以上的步骤是由于所获取到的区域海拔数据并不是准确无误的，海拔高度信息会存在一定的误差，会给以后的检测带来不利的影响，故需要尽量使海拔高度信息准确。由于误差的出现具有随机性，有时海拔高度比实际的高，有时海拔高度比实际的低，所以利用其随机性消除这种误差，其原理主要是基于以下的公式：

$\begin{array}{c}\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}(H_i+{\mathrm{\&epsiv;}}_i)\\ =\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{H}_i+\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}_i\\ =\frac{1}{m*m}\underset{i=1}{\overset{m*m}{\mathrm{\&Sigma;}}}{H}_i\end{array}$

式中H_i+ε_i为读取的海拔高度，H_i为实际的海拔高度，ε_i为高度误差。从上式中可以看出，由于误差的随机性，平均误差最后得到比较准确的区域海拔高度信息数据。本实例中，m为5。

(2)利用太阳方位角等图像信息对校正后的海拔高度图进行旋转处理，通过太阳高度角等图像信息计算山体本体和山体阴影的比例，再经过一系列的处理得基于DEM的图像阴影位置，具体过程主要包括以下子步骤：

(2.1)利用太阳方位角旋转海拔高度图步骤

利用产生仿真图像时设置的参数确定太阳方位角θ，太阳方位角如图6所示，然后对校正后的海拔高度图旋转太阳方位角θ。如图5所示，旋转可以看做是从一个坐标系变换到了另一个坐标系的仿射变换，为了减小旋转带来的误差，采用双线性插值的方法对图像进行旋转处理。仿射变换的计算公式如下：

最后得到旋转后的海拔高度图，超出的部分海拔高度设置为0。

上述中太阳方位角为竖立在地面上的直线的阴影和正南方的夹角。本实例中，太阳方位角θ＝234.65。

(2.2)计算山体本体和山体投影的比例步骤

其中阴影位置和太阳高度角也有关系，因为太阳高度角决定山体本体和山体投影的比例，这个比例直接决定阴影的大小，太阳高度角如图7所示，太阳高度角的计算公式为：

φ代表当地的地理纬度，代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+。本实例中，地理纬度和太阳直射点纬度在同一个半球，得β＝29。

山体本体和山体阴影的关系为：

由上式可以得到投影的范围大小，从而可以确定具体阴影的大小和位置。本实例中，tanβ＝0.554

(2.3)阴影的区域的检测标记步骤

(2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：

$H(:,i)\&DoubleLeftRightArrow;H(:,n-i)(i=\mathrm{1,2,3}...\frac{2}{n})$

式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，代表交换该符号前后的数值。

将旋转后得到的图进行左右倒置处理，是为了保证小的山丘在大的山丘的阴影部位的时候不被检测为非阴影区域，从而相当于太阳是从东边照射过来的。

(2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，为了对其中的每一个点进行判断，判断该点是否存在于阴影区，判断的公式如下：

H'(i,j)≠0 (i＝1,2,3...m、j＝1,2,3...n)

$\begin{array}{cc}(k-j)*10<(H^{\&prime;}(i,k)-H^{\&prime;}(i,j))*\frac{1}{\tan \mathrm{\&beta;}}& (k=j+1,j+2...n)\end{array}$

式中H^’(i，j)代表旋转倒置后的海拔高度图，m和n分别代表海拔高度图的行数和列数。

如果满足以上等式条件，就标记该点为阴影区域的位置：

label(i,j)＝1。

(2.3.3)将处理后的标记图左右倒置回来，将倒置回来的标记图进行旋转，旋转的角度就是360-太阳方位角θ。将旋转后得到的标记图大小和山体的红外遥感图像F(x，y)对应起来，就可以根据标记图对山体红外遥感图像中的阴影部分进行标记处理。本实例中，旋转的角度为125.35。如图8(a)和图9(a)所示，分别是旋转后的仿真图1海拔高度和仿真图2海拔高度阴影标记图。

(3)通过对得到的阴影部分的标记图进行膨胀处理，对膨胀后阴影区域进行一定的补偿，具体过程主要包括以下子步骤：

(3.1)阴影标记图的膨胀处理步骤

对得到的标记图进行膨胀处理以得到连续的区域。膨胀处理就是用一个结构元素在标记图上移动并进行比较，根据比较的结果做处理，处理的过程为：

a、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点没有一个相同，该点就为黑色，即非阴影区域。

b、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点有一个相同，该点就为白色，即阴影区域。

经过上述的处理过程就可以得到膨胀后的标记图，将不连续的区域变得连续，更加符合实际。如图8(b)和图9(b)所示，分别为对仿真图像1的阴影标记图和仿真图像2的阴影标记图膨胀处理后的结果图。

(3.2)对阴影区域进行补偿处理步骤

假定图像是局部平稳的，可以认为阴影区域与其周围一定范围内的非阴影区域的地质结构同质，红外特性相似。如图10(a)和图11(a)所示，仿真隧道1区域和仿真隧道2区域存在明显的阴影分布，故在检测出阴影区域之后，对各个独立阴影区域与其邻近的非阴影区域进行匹配，完成阴影区域补偿操作，从而达到去除阴影的目的。算法过程如下：

1)对于阴影区域，通过下式计算其邻近的非阴影区域

Ω_noshadow＝{p|0<d(p,Ω_shadow)<dist}

式中Ω_noshadow表示邻近某个距离阈值dist的非阴影区域集合，dist＝50。

2)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：

${F^{\&prime;}}_{\mathrm{shadow}}(i,j)=A[m_{\mathrm{noshadow}}+\frac{F_{\mathrm{shadow}}(i,j)-m_{\mathrm{shadow}}}{{\mathrm{\&sigma;}}_{\mathrm{shadow}}}*{\mathrm{\&sigma;}}_{\mathrm{noshadow}}]$

D＝F'_shadow(i,j)-F_shadow(i,j)

式中F_shadow(i,j)是补偿之前的阴影区域灰度值，F'_shadow(i,j)是补偿之后阴影区域灰度值，m_shadow和σ_shadow是阴影区域灰度值的均值和方差，m_noshadow和σ_noshadow是邻近非阴影区域灰度值的均值和方差，A为补偿强度系数，D为阴影区域需要的补偿值。本实例中，m_shadow＝1465，m_noshadow＝1503，σ_shadow＝55.2752，σ_noshadow＝61.9714，D＝44，A＝1.0。

从上述中得到需要补偿的灰度值后，遍历所有检测到的阴影区域，将阴影区域对应到红外遥感图像每一个点的灰度值补偿D，得到最后矫正后的灰度图。如图10(b)和图11(b)所示为仿真隧道1区域和仿真隧道2区域在阴影补偿后的红外遥感图像，从得到的红外遥感图像可以看出，在一定程度上消除了阴影。

(4)对阴影区域补偿后的红外遥感图像进行初步的探测识别，然后利用位置的关联性进行二次的判断，提高识别率，降低虚警率，具体的过程包含以下子步骤：

(4.1)初步判断带状地下结构位置步骤

利用阴影补偿后的红外遥感图像探测带状地下结构，依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别，每段的长度为250米宽度为30米(即25个像素点长3个像素点宽)，对像素点的左右遍历，探测是否存在带状地下结构，最后得到初步的隧道识别结果。同时在垂直于隧道的方向利用相同的方法探测地下结构，计算虚警率。

(4.2)根据位置关联性进行二次判断步骤

对初步识别出的符合脉冲模式的分散段进行统计分析，依据500米长度(50个像素点)上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白。

1)若根据判断可以忽视掉这些分散脉冲模式段中间的小空白，那么该500的长度段就确认为带状目标所在位置。

2)若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白，那么该500的长度段就确认为非带状目标所在位置。

如图12(a)和图12(b)所示，通过对两幅图的对比可以看到对阴影补偿后的图进行带状地下结构检测结果比未补偿的检测结果要好。

通过对为去除阴影和去除阴影的红外遥感图像进行探测识别，得到了以下对比表，可以看出去除阴影后，识别率有所提高，同时虚警率也有所降低。表格如下：

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于太阳照射阴影补偿的带状地下结构探测方法，其特征在于，所述方法包括如下步骤：

(1)根据待检测区域的地理位置确定待检测区域的经纬度信息，然后根据红外遥感图像的分辨率计算经纬度步长，并根据待检测区域的经纬度信息和经纬度步长获取待检测区域的海拔高度信息，最后对所得到的海拔高度信息进行校正处理；

(2)利用太阳方位角信息对校正后的海拔高度图进行旋转处理，通过太阳高度角信息计算山体本体和山体阴影的比例，确定海拔高度图中阴影的大小和位置，对旋转后的海拔高度图进行处理得到图像阴影位置的标记图；

(3)对得到的图像阴影位置的标记图进行膨胀处理，并对膨胀后的阴影区域进行补偿；

(4)对阴影区域补偿后的红外遥感图像进行初步的探测识别，然后利用位置的关联性进行二次的判断，从而探测得到带状地下结构。

2.如权利要求1所述的方法，其特征在于，所述步骤(1)中根据待检测区域的地理位置确定待检测区域的经纬度信息，然后根据红外遥感图像的分辨率计算经纬度步长，具体包括如下子步骤：

(1.1.1)根据待检测区域的地理位置获得待检测区域左上方点经纬度信息为A(x₁，y₁)，以及右下方点的经纬度信息为B(x₂，y₂)；

(1.1.2)获得红外遥感图像的分辨率为k米，分别根据下式求取经度步长lon和纬度步长lat，

C＝sin(y₁)²*cos(x₁-x₂)+cos(y₁)²   (y₁＝y₂)

d₁＝R*arccos(C)*Pi/180

R＝6371.004

lon＝k*(x₂-x₁)/d₁

C＝cos(y₁-y₂)   (x₁＝x₂)

d₂＝R*arccos(C)*Pi/180

R＝6371.004。

lat＝k*(y₁-y₂)/d₂

3.如权利要求1或2所述的方法，其特征在于，所述步骤(1)中根据待检测区域的经纬度信息和经纬度步长获取待检测区域的海拔高度信息，具体为：根据待检测区域的经纬度信息和经纬度步长得到待检测区域每一点的经纬度信息，根据每一点的经纬度信息在Google Earth中进行定位，再获取每一个定位点的海拔信息，得到一张和经纬度图一样大小的海拔信息图，其分辨率同样为k米。

4.如权利要求1或2所述的方法，其特征在于，所述步骤(1)中对所得到的海拔高度信息进行校正处理，具体包括：

(1.3.1)通过下式对海拔信息图处理求出海拔高度中的奇异点，其中奇异点和非奇异点满足以下条件：

H(i,j)>H(i+k,j+l)+h   (k,l＝-1,1)式中(i，j)为奇异点

H(i,j)<H(i+k,j+l)+h   (k,l＝-1,1)式中(i，j)为非奇异点

式中H(i，j)为(i，j)点处的海拔高度，H(i+k，j+l)为(i，j)邻近点的海拔高度；

对奇异点进行中值滤波处理，去除奇异点，得到初步校正后的海拔高度图；

(1.3.2)对初步校正后的海拔高度图进行m*m模板大小的均值滤波，对海拔高度图进行再次校正。

5.如权利要求1或2所述的方法，其特征在于，所述步骤(2)中利用太阳方位角信息对校正后的海拔高度图进行旋转处理，具体包括：

确定太阳方位角θ，对校正后的海拔高度图旋转太阳方位角θ，变换公式如下：

旋转前 $\left\{\begin{array}{c}{x}_0=r\cos \mathrm{\&alpha;}\\ y_0=r\sin \mathrm{\&alpha;}\end{array}\right.$

旋转后 $\left\{\begin{array}{c}{x}_1=r\cos (\mathrm{\&alpha;}-\mathrm{\&theta;})=r\cos \mathrm{\&alpha;}\cos \mathrm{\&theta;}+r\sin \mathrm{\&alpha;}\sin \mathrm{\&theta;}=x_0\cos \mathrm{\&theta;}+y_0\sin \mathrm{\&theta;}\\ y_1=r\sin (\mathrm{\&alpha;}-\mathrm{\&theta;})=r\sin \mathrm{\&alpha;}\cos \mathrm{\&theta;}-r\cos \mathrm{\&alpha;}\sin \mathrm{\&theta;}=-x_0\sin \mathrm{\&theta;}+y_0\cos \mathrm{\&theta;}\end{array}\right.$

对旋转后得到的海拔高度图，超出的部分海拔高度设置为0；

其中，所述中太阳方位角为地面的竖立直线的阴影和正南方的夹角。

6.如权利要求1或2所述的方法，其特征在于，所述步骤(2)中通过太阳高度角信息计算山体本体和山体阴影的比例，确定海拔高度图中阴影的大小和位置，具体包括：

根据公式计算太阳高度角，其中φ代表当地的地理纬度，代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+；

根据山体阴影和山体本体的关系求得阴影的大小和位置。

7.如权利要求1或2所述的方法，其特征在于，所述步骤(2)中对旋转后的海拔高度图进行处理得到图像阴影位置的标记图，具体包括：

(2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：

$H(:,i)\&DoubleLeftRightArrow;H(:,n-i)(i=\mathrm{1,2,3}...\frac{2}{n})$

式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，代表交换该符号前后的数值；

(2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，判断每一个点是否存在于阴影区，判断公式如下：

H'(i,j)≠0(i＝1,2,3...m、j＝1,2,3...n)

$(k-j)*10<(H^{\&prime;}(i,k)-H^{\&prime;}(i,j))*\frac{1}{\tan \mathrm{\&beta;}}(k=j+1,j+2...n)$

式中H^’(i，j)代表旋转倒置后的海拔高度图，m和n分别代表海拔高度图的行数和列数，如果一个点满足以上等式的条件，就将该点标记为阴影区域。

(2.3.3)将处理后的标记图左右倒置回来，将倒置回来的标记图进行旋转，旋转的角度为360-太阳方位角θ，得到图像阴影位置的标记图。

8.如权利要求1或2所述的方法，其特征在于，所述步骤(3)中对得到的阴影部分的标记图进行膨胀处理，具体包括：

对得到的阴影部分的标记图进行膨胀处理以得到连续的区域，其中膨胀处理具体为，用一个结构元素在标记图上移动并进行比较，根据比较的结果做处理，处理过程为：

a、如果结构元素中的所有阴影点与它对应的标记图的阴影区域像素点没有一个相同，该点就为黑色，即非阴影区域；

b、如果结构元素中的所有阴影点与它对应的标记图的阴影区域像素点有一个相同，该点就为白色，即阴影区域。

9.如权利要求1或2所述的方法，其特征在于，所述步骤(3)中对对膨胀后阴影区域进行补偿，具体包括：

1)对于阴影区域，通过下式计算其邻近的非阴影区域

Ω_noshadow＝{p|0<d(p,Ω_shadow)<dist}

式中Ω_noshadow表示邻近某个距离阈值dist的非阴影区域集合；

2)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：

${F^{\&prime;}}_{\mathrm{shadow}}(i,j)=A[m_{\mathrm{noshadow}}+\frac{F_{\mathrm{shadow}}(i,j)-m_{\mathrm{shadow}}}{{\mathrm{\&sigma;}}_{\mathrm{shadow}}}*{\mathrm{\&sigma;}}_{\mathrm{noshadow}}]$

D＝F'_shadow(i,j)-F_shadow(i,j)

式中F_shadow(i,j)是补偿之前的阴影区域灰度值，F'_shadow(i,j)是补偿之后阴影区域灰度值，m_shadow和σ_shadow是阴影区域灰度值的均值和方差，m_noshadow和σ_noshadow是邻近非阴影区域灰度值的均值和方差，A为补偿强度系数，D为阴影区域需要的补偿值；

得到需要补偿的灰度值D后，遍历所有检测到的阴影区域，将阴影区域对应到红外遥感图像上每一个点的灰度值补偿D，得到最终补偿后的红外遥感图像。

10.如权利要求1或2所述的方法，其特征在于，所述步骤(4)具体包括：

(4.1)初步判断带状地下结构位置

利用阴影补偿后的红外遥感图像探测带状地下结构，依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别，对像素点的左右遍历，探测是否存在带状地下结构，最后得到初步的隧道识别结果；同时在垂直于隧道的方向利用相同的方法探测地下结构，计算虚警率；

(4.2)根据位置关联性进行二次判断

对初步识别出的符合脉冲模式的分散段进行统计分析，依据预设长度上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白：若根据判断可以忽视掉这些分散脉冲模式段中间的小空白，那么该预设长度段就确认为带状目标所在位置；若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白，那么该预设长度段就确认为非带状目标所在位置。