[0075] 式中Ωnoshadow表示邻近某个距离阈值dist的非阴影区域集合,
[0076] 2)在得到阴影区域和其邻近的非阴影区域之后,采用如下的映射策略对阴影区域的灰度值进行补偿:
[0077]
[0078] D=F'shadow(i,j)-Fshadow(i,j)
[0079] 式中Fshadow(i,j)是补偿之前的阴影区域灰度值,F'shadow(i,j)是补偿之后阴影区域灰度值,mshadow和σshadow是阴影区域灰度值的均值和方差,mnoshadow和σnoshadow是邻近非阴影区域灰度值的均值和方差,A为补偿强度系数,D为阴影区域需要的补偿值。
[0080] 从上述中得到需要补偿的灰度值D后,遍历所有检测到的阴影区域,将阴影区域对应到红外遥感图像上每一个点的灰度值补偿D,得到最终补偿后的红外遥感图像。
[0081] (4)对阴影区域补偿后的红外遥感图像进行初步的探测识别,然后利用位置的关联性进行二次的判断,提高识别率,降低虚警率,具体的过程包含以下子步骤:
[0082] (4.1)初步判断带状地下结构位置步骤
[0083] 利用阴影补偿后的红外遥感图像探测带状地下结构,依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别,每段的长度为250米宽度为30米(即25个像素点长3个像素点宽),对像素点的左右遍历,探测是否存在带状地下结构,最后得到初步的隧道识别结果。同时在垂直于隧道的方向利用相同的方法探测地下结构,计算虚警率。
[0084] (4.2)根据位置关联性进行二次判断步骤
[0085] 对初步识别出的符合脉冲模式的分散段进行统计分析,依据500米长度(50个像素点)上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白。
[0086] 1)若根据判断可以忽视掉这些分散脉冲模式段中间的小空白,那么该500的长度段就确认为带状目标所在位置。
[0087] 2)若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白,那么该500的长度段就确认为非带状目标所在位置。
[0088] 与现有技术相比,本发明具有以下有益效果:
[0089] 1、利用获取的指定区域的DEM地形数据和太阳方位角、太阳高度角等有用信息,可以有效的将静态红外图像中太阳照射的阴影区域检测出来,并利用阴影区域与其周围一定范围内的非阴影区域的地质结构同质,红外特性相似等知识,对阴影区域进行了一定的补偿处理,这样可以减小阴影对地下结构探测的不利影响。
[0090] 2、对太阳照射的阴影补偿后的结果图进行地下结构探测,相对于对未进行阴影补偿的红外图进行地下结构探测来说,能增加地下结构的探测正确率,同时降低了虚警率。
附图说明
[0091] 图1为本发明基于太阳照射阴影补偿的带状地下结构探测方法流程示意图;
[0092] 图2为本发明实施例中经纬度和距离的示意图;
[0093] 图3为本发明实施例中未处理奇异点前的海拔高度图;
[0094] 图4为本发明实施例中处理奇异点后海拔高度图;
[0095] 图5为本发明实施例中海拔高度旋转过程示意图;
[0096] 图6为本发明实施例中太阳方位角示意图;
[0097] 图7为本发明实施例中太阳高度角示意图;
[0098] 图8(a)为本发明实施例中仿真隧道1区域的阴影标记图;
[0099] 图8(b)为本发明实施例中仿真隧道1区域的阴影标记膨胀图;
[0100] 图9(a)为本发明实施例中仿真隧道2区域的阴影标记图;
[0101] 图9(b)为本发明实施例中仿真隧道2区域的阴影标记膨胀图;
[0102] 图10(a)为本发明实施例中仿真隧道1区域的红外遥感图像;
[0103] 图10(b)为本发明实施例中仿真隧道1区域的阴影补偿红外遥感图像;
[0104] 图11(a)为本发明实施例中仿真隧道2区域的红外遥感图像;
[0105] 图11(b)为本发明实施例中仿真隧道2区域的阴影补偿红外遥感图像;
[0106] 图12(a)为本发明实施例中未补偿除阴影的仿真隧道1的检测结果图;
[0107] 图12(b)为本发明实施例中补偿阴影后的仿真隧道1的检测结果图。
具体实施方式
[0108] 为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。此外,下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0109] 本发明的流程如图1所示,其中具体的实施方法包括以下步骤。包括:获取指定区域的DEM地形数据步骤、利用DEM和太阳高度角、太阳方位角获取图像阴影位置步骤、对阴影区域进行处理和补偿步骤、阴影区域校正后带状地下结构的探测步骤:
[0110] (1)根据区域的地理位置得到区域的经纬度信息,然后通过图像的分辨率计算经纬度步长,获取指定区域的海拔高度信息,最后对所得到的海拔高度信息进行校正处理,具体过程包括以下子步骤:
[0111] (1.1)确定指定区域的经纬度步骤
[0112] (1.1.1)经纬度信息的确定可以根据需要检测的位置和所需数据的分辨率来确定。如图2所示,本实例中,待检测区域左上方点经纬度信息为A(x1,y1)=(116.150049,40.296833),右下方点的经纬度信息为
[0113] B(x2,y2)=(116.194775,40.260787)。
[0114] (1.1.2)所仿真的红外遥感图像的分辨率为k米,需要将DEM地形数据的分辨率也设置为k米,地球上任意两点的实际距离d满足一定的关系,关系式如下:
[0115] C=sin(y1)*sin(y2)*cos(x1-x2)+cos(y1)*cos(y2)
[0116] D=R*arccos(C)*Pi/180
[0117] R=6371.004
[0118] 上式中(x1,y1)和(x2,y2)分别代表这两点的经纬度信息。
[0119] 如图2所示,故平行于赤道时的距离d1和经纬度的关系如下:
[0120] C=sin(y1)2*cos(x1-x2)+cos(y1)2 (y1=y2)
[0121] d1=R*arccos(C)*Pi/180
[0122] R=6371.004
[0123] lon=k*(x2-x1)/d1
[0124] 本实例中,d1=3800m,lon=0.0001177。
[0125] 垂直于赤道时的距离d2和经纬度的关系如下:
[0126] C=cos(y1-y2) (x1=x2)
[0127] d2=R*arccos(C)*Pi/180
[0128] R=6371.004
[0129] lat=k*(y1-y2)/d2
[0130] 本实例中,d2=4000m,lat=0.000090115。
[0131] (1.2)利用Google Earth获取所需区域的地形数据步骤
[0132] 根据所获取的经纬度信息,利用Google Earth提供的编程接口进行编程,首先根据每一点的经纬度信息在Google Earth中进行定位,然后再定位成功后获取定位点的海拔信息。最后就会得到一张和经纬度图一样大小的海拔信息图,其分辨率同样为k米。
[0133] (1.3)地形数据的矫正处理步骤
[0134] (1.3.1)如图3所示,通过对图3的海拔高度信息观察发现,其中有许多的奇异点,有些位置的海拔高度明显高于邻近点的海拔高度,为了减少误差,通过和邻近的比较来确定是否为奇异点,其中奇异点和非奇异点满足以下条件:
[0135] H(i,j)>H(i+k,j+l)+h (k,l=-1,1)式中(i,j)为奇异点
[0136] H(i,j)
[0137] 式中H(i,j)为(i,j)点处的海拔高度,H(i+k,j+l)为(i,j)邻近点的海拔高度。本实例中,h为20。
[0138] 由于中值滤波对奇异点的处理效果很好,故对奇异点进行中值滤波处理,可以去除掉奇异点,得到初步校正后海拔高度图的。结果如图4所示。
[0139] (1.3.2)对初步校正后的海拔高度图再次校正,利用m*m大小的模板对海拔高度图进行均值滤波处理。
[0140] 进行以上的步骤是由于所获取到的区域海拔数据并不是准确无误的,海拔高度信息会存在一定的误差,会给以后的检测带来不利的影响,故需要尽量使海拔高度信息准确。由于误差的出现具有随机性,有时海拔高度比实际的高,有时海拔高度比实际的低,所以利用其随机性消除这种误差,其原理主要是基于以下的公式:
[0141]
[0142] 式中Hi+εi为读取的海拔高度,Hi为实际的海拔高度,εi为高度误差。从上式中可以看出,由于误差的随机性,平均误差 最后得到比较准确的区域海拔高度信息数据。本实例中,m为5。
[0143] (2)利用太阳方位角等图像信息对校正后的海拔高度图进行旋转处理,通过太阳高度角等图像信息计算山体本体和山体阴影的比例,再经过一系列的处理得基于DEM的图像阴影位置,具体过程主要包括以下子步骤:
[0144] (2.1)利用太阳方位角旋转海拔高度图步骤
[0145] 利用产生仿真图像时设置的参数确定太阳方位角θ,太阳方位角如图6所示,然后对校正后的海拔高度图旋转太阳方位角θ。如图5所示,旋转可以看做是从一个坐标系变换到了另一个坐标系的仿射变换,为了减小旋转带来的误差,采用双线性插值的方法对图像进行旋转处理。仿射变换的计算公式如下:
[0146]
[0147]
[0148] 最后得到旋转后的海拔高度图,超出的部分海拔高度设置为0。
[0149] 上述中太阳方位角为竖立在地面上的直线的阴影和正南方的夹角。本实例中,太阳方位角θ=234.65。
[0150] (2.2)计算山体本体和山体投影的比例步骤
[0151] 其中阴影位置和太阳高度角也有关系,因为太阳高度角决定山体本体和山体投影的比例,这个比例直接决定阴影的大小,太阳高度角如图7所示,太阳高度角的计算公式为:
[0152]
[0153] φ代表当地的地理纬度,代表太阳直射点地理纬度,(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球:如果在同一半球符号为-,在南北半球符号为+。本实例中,地理纬度和太阳直射点纬度在同一个半球,得β=29。
[0154] 山体本体和山体阴影的关系为:
[0155]
[0156] 由上式可以得到投影的范围大小,从而可以确定具体阴影的大小和位置。本实例中,tanβ=0.554
[0157] (2.3)阴影的区域的检测标记步骤
[0158] (2.3.1)将旋转后得到的海拔高度图左右倒置,具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置,遍历的过程如下:
[0159]
[0160] 式中H(:,i)代表第i列的所有行,n为海拔高度图的列数, 代表交换该符号前后的数值。
[0161] 将旋转后得到的图进行左右倒置处理,是为了保证小的山丘在大的山丘的阴影部位的时候不被检测为非阴影区域,从而相当于太阳是从东边照射过来的。
[0162] (2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历,为了对其中的每一个点进行判断,判断该点是否存在于阴影区,判断的公式如下:
[0163] H'(i,j)≠0 (i=1,2,3...m、j=1,2,3...n)
[0164]
[0165] 式中H’(i,j)代表旋转倒置后的海拔高度图,m和n分别代表海拔高度图的行数和列数。
[0166] 如果满足以上等式条件,就标记该点为阴影区域的位置:
[0167] label(i,j)=1。
[0168] (2.3.3)将处理后的标记图左右倒置回来,将倒置回来的标记图进行旋转,旋转的角度就是360-太阳方位角θ。将旋转后得到的标记图大小和山体的红外遥感图像F(x,y)对应起来,就可以根据标记图对山体红外遥感图像中的阴影部分进行标记处理。本实例中,旋转的角度为125.35。如图8(a)和图9(a)所示,分别是旋转后的仿真图1海拔高度和仿真图2海拔高度阴影标记图。
[0169] (3)通过对得到的阴影部分的标记图进行膨胀处理,对膨胀后阴影区域进行一定的补偿,具体过程主要包括以下子步骤:
[0170] (3.1)阴影标记图的膨胀处理步骤
[0171] 对得到的标记图进行膨胀处理以得到连续的区域。膨胀处理就是用一个结构元素在标记图上移动并进行比较,根据比较的结果做处理,处理的过程为:
[0172] a、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点没有一个相同,该点就为黑色,即非阴影区域。
[0173] b、如果结构元素中的所有白色点(阴影点)与它对应的标记图的阴影区域像素点有一个相同,该点就为白色,即阴影区域。
[0174] 经过上述的处理过程就可以得到膨胀后的标记图,将不连续的区域变得连续,更加符合实际。如图8(b)和图9(b)所示,分别为对仿真图像1的阴影标记图和仿真图像2的阴影标记图膨胀处理后的结果图。
[0175] (3.2)对阴影区域进行补偿处理步骤
[0176] 假定图像是局部平稳的,可以认为阴影区域与其周围一定范围内的非阴影区域的地质结构同质,红外特性相似。如图10(a)和图11(a)所示,仿真隧道1区域和仿真隧道2区域存在明显的阴影分布,故在检测出阴影区域之后,对各个独立阴影区域与其邻近的非阴影区域进行匹配,完成阴影区域补偿操作,从而达到去除阴影的目的。算法过程如下:
[0177] 1)对于阴影区域,通过下式计算其邻近的非阴影区域
[0178] Ωnoshadow={p|0
[0179] 式中Ωnoshadow表示邻近某个距离阈值dist的非阴影区域集合,dist=50。
[0180] 2)在得到阴影区域和其邻近的非阴影区域之后,采用如下的映射策略对阴影区域的灰度值进行补偿:
[0181]
[0182] D=F'shadow(i,j)-Fshadow(i,j)
[0183] 式中Fshadow(i,j)是补偿之前的阴影区域灰度值,F'shadow(i,j)是补偿之后阴影区域灰度值,mshadow和σshadow是阴影区域灰度值的均值和方差,mnoshadow和σnoshadow是邻近非阴影区域灰度值的均值和方差,A为补偿强度系数,D为阴影区域需要的补偿值。本实例中,mshadow=1465,mnoshadow=1503,σshadow=55.2752,σnoshadow=61.9714,D=44,A=1.0。
[0184] 从上述中得到需要补偿的灰度值后,遍历所有检测到的阴影区域,将阴影区域对应到红外遥感图像每一个点的灰度值补偿D,得到最后矫正后的灰度图。如图10(b)和图11(b)所示为仿真隧道1区域和仿真隧道2区域在阴影补偿后的红外遥感图像,从得到的红外遥感图像可以看出,在一定程度上消除了阴影。
[0185] (4)对阴影区域补偿后的红外遥感图像进行初步的探测识别,然后利用位置的关联性进行二次的判断,提高识别率,降低虚警率,具体的过程包含以下子步骤:
[0186] (4.1)初步判断带状地下结构位置步骤
[0187] 利用阴影补偿后的红外遥感图像探测带状地下结构,依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别,每段的长度为250米宽度为30米(即25个像素点长3个像素点宽),对像素点的左右遍历,探测是否存在带状地下结构,最后得到初步的隧道识别结果。同时在垂直于隧道的方向利用相同的方法探测地下结构,计算虚警率。
[0188] (4.2)根据位置关联性进行二次判断步骤
[0189] 对初步识别出的符合脉冲模式的分散段进行统计分析,依据500米长度(50个像素点)上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白。
[0190] 1)若根据判断可以忽视掉这些分散脉冲模式段中间的小空白,那么该500的长度段就确认为带状目标所在位置。
[0191] 2)若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白,那么该500的长度段就确认为非带状目标所在位置。
[0192] 如图12(a)和图12(b)所示,通过对两幅图的对比可以看到对阴影补偿后的图进行带状地下结构检测结果比未补偿的检测结果要好。
[0193] 通过对为去除阴影和去除阴影的红外遥感图像进行探测识别,得到了以下对比表,可以看出去除阴影后,识别率有所提高,同时虚警率也有所降低。表格如下:
[0194]
[0195] 本领域的技术人员容易理解,以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。