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

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

Info

Publication number
WO2016106950A1
WO2016106950A1 PCT/CN2015/072663 CN2015072663W WO2016106950A1 WO 2016106950 A1 WO2016106950 A1 WO 2016106950A1 CN 2015072663 W CN2015072663 W CN 2015072663W WO 2016106950 A1 WO2016106950 A1 WO 2016106950A1
Authority
WO
WIPO (PCT)
Prior art keywords
shadow
area
altitude
latitude
map
Prior art date
Application number
PCT/CN2015/072663
Other languages
English (en)
French (fr)
Inventor
张天序
郝龙伟
马文绚
鲁岑
药珩
王岳环
桑农
杨卫东
Original Assignee
华中科技大学
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 华中科技大学 filed Critical 华中科技大学
Priority to US15/106,686 priority Critical patent/US9582885B2/en
Publication of WO2016106950A1 publication Critical patent/WO2016106950A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V9/00—Prospecting or detecting by methods not provided for in groups G01V1/00 - G01V8/00
    • G01V9/005—Prospecting or detecting by methods not provided for in groups G01V1/00 - G01V8/00 by thermal methods, e.g. after generation of heat by chemical reactions
    • 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
    • G06T5/00—Image enhancement or restoration
    • G06T5/90—Dynamic range modification of images or parts thereof
    • G06T5/94—Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/60—Extraction of image or video features relating to illumination properties, e.g. using a reflectance or lighting model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10048—Infrared image

Definitions

  • the invention belongs to the intersection of remote sensing technology, image processing and pattern recognition, and particularly relates to a method for detecting a strip-shaped underground structure based on sun-shadow compensation.
  • Shadow is a physics phenomenon that is ubiquitous in nature and is caused by the light source being blocked by an object.
  • the presence of shadows in images has different effects on related issues in the field of computer vision. This effect has advantages and disadvantages. For example, adding shadows to objects in virtual reality and 3D games can improve the realism of the scene. But more often, shadows in images can adversely affect computer vision related issues. For example, in remote sensing images, the presence of shadows can affect subsequent image matching, pattern recognition, and feature extraction. Processing operations; in video surveillance, shadows and moving targets are combined, causing errors in the extraction and tracking of the target object by the computer. Therefore, it is necessary to detect and analyze the shadows in the image and to eliminate or reduce the effects of the shadows as needed.
  • the invention provides a method for detecting a strip-shaped underground structure based on sun-shadow compensation, and combines a digital height model (DEM) to perform shadow detection on a single infrared remote sensing image, and performs gray scale on the region where the shadow is located.
  • the compensation which reduces the effect of shadows on the detection.
  • the invention provides a method for detecting a strip-shaped underground structure based on sun-shadow compensation, and the specific steps are as follows:
  • the specific process includes the following sub step:
  • the determination of the latitude and longitude information can be determined according to the position to be detected and the resolution of the required data.
  • the latitude and longitude information of the left upper point of the to-be-detected area is A(x 1 , y 1 ), and the latitude and longitude information of the lower right point is B(x 2 , y 2 ).
  • the resolution of the infrared remote sensing image is k meters, so the resolution of the DEM terrain data needs to be set to k meters.
  • the actual distance D of any two points on the earth satisfies a certain relationship, and the relationship is as follows:
  • (x 1 , y 1 ) and (x 2 , y 2 ) represent the latitude and longitude information of these two points, respectively.
  • d 1 is the distance between two points parallel to the equator and lon is the longitude step.
  • d 2 is the distance between two points perpendicular to the equator and lat is the latitude step.
  • the latitude and longitude information of each point of the area to be detected is obtained, and the latitude and longitude resolution is k meters.
  • the programming interface provided by Google Earth is used for programming. Firstly, the latitude and longitude information of each point is used for positioning in Google Earth, and then the altitude information of the anchor point is obtained after the positioning is successful. Finally, you will get an elevation map of the same size as the latitude and longitude map, with a resolution of k meters.
  • H(i,j) is the altitude at the point of (i,j)
  • H(i+k,j+l) is the altitude of the adjacent point of (i,j).
  • the median filtering has a good effect on the processing of singular points, the median filtering process is performed on the singular points, and the singular points can be removed to obtain the preliminary corrected altitude map.
  • the above steps are performed because the acquired altitude data of the area is not accurate, and there is a certain error in the altitude information, which will adversely affect the subsequent detection. Therefore, it is necessary to make the altitude information as accurate as possible. Due to the randomness of the error, sometimes the altitude is higher than the actual height, and sometimes the altitude is lower than the actual height, so the randomness is used to eliminate this error.
  • the principle is mainly based on the following formula:
  • the solar azimuth angle ⁇ is determined by the parameters set when the simulated image is generated, and then the solar azimuth angle ⁇ is rotated for the corrected altitude map, and the rotation can be regarded as an affine transformation from one coordinate system to another coordinate system, in order to reduce The error caused by small rotation is rotated by the method of bilinear interpolation.
  • the calculation formula of the affine transformation is as follows:
  • the above-mentioned medium azimuth angle is the shadow of a straight line standing on the ground and the angle between the south and the south.
  • the shadow position is also related to the sun height angle, because the sun height angle determines the proportion of the mountain body and the mountain projection. This ratio directly determines the size of the shadow, and the solar height angle is calculated as:
  • ⁇ represents the local geographic latitude
  • (+/-) is whether the geographic latitude of the requested geographic latitude and the direct sun point is in the same hemisphere: if the symbol in the same hemisphere is -, the symbol in the northern and southern hemispheres is +.
  • the range of the projection can be obtained, so that the size and position of the specific shadow can be determined.
  • H(:,i) represents all rows of the i-th column
  • n is the number of columns in the altitude map. Represents the value before and after the exchange of the symbol.
  • the image obtained by the rotation is inverted to the left and right to ensure that the small hill is not detected as a non-shaded area when it is in the shadow of the large hill, so that the sun is irradiated from the east.
  • H'(i,j) represents the altitude map after the rotation is inverted
  • m and n represent the number of rows and columns of the altitude map, respectively.
  • the shaded area after the expansion is compensated to some extent, and the specific process mainly includes the following substeps:
  • the resulting marker image is expanded to obtain a continuous region.
  • the expansion process uses a structural element to move and compare on the marker map, and processes it according to the result of the comparison.
  • the process of processing is:
  • the point is black, that is, the non-shaded area.
  • the expanded mark pattern can be obtained, and the discontinuous area becomes continuous, which is more practical.
  • each independent shadow area is matched with its adjacent non-shadow area to complete the shadow area compensation operation, thereby achieving the purpose of removing the shadow.
  • the algorithm process is as follows:
  • ⁇ noshadow ⁇ p
  • ⁇ noshadow represents a set of non-shaded regions adjacent to a certain distance threshold dist
  • mapping strategy is used to compensate the gray value of the shaded area:
  • F shadow (i,j) is the gray value of the shadow area before compensation
  • F' shadow (i,j) is the gray value of the shadow area after compensation
  • m shadow and ⁇ shadow are the mean values of the gray value of the shadow area.
  • the variance, m noshadow and ⁇ noshadow are the mean and variance of the gray values of the adjacent non-shadow areas
  • A is the compensation intensity coefficient
  • D is the compensation value required for the shadow area.
  • the specific process includes the following sub-steps:
  • the infrared remote sensing image after shadow compensation is used to detect the strip-shaped underground structure, and the thermal infrared remote sensing image is detected by pixel-by-pixel rolling segmentation according to the approximate position and direction information of the tunnel.
  • the length of each segment is 250 meters and the width is 30 meters (ie 25
  • the pixel points are 3 pixels wide, and the left and right traversals of the pixels are detected to detect whether there is a band-like underground structure, and finally the preliminary tunnel recognition result is obtained.
  • the same method is used to detect the underground structure in the direction perpendicular to the tunnel, and the false alarm rate is calculated.
  • the present invention has the following beneficial effects:
  • the shadow area of the sun illumination in the static infrared image can be effectively detected, and the shadow area and a certain range around it can be utilized.
  • the geological structure of the shaded area is homogenous and the infrared characteristics are similar.
  • the shadow area is compensated to a certain extent, which can reduce the adverse effects of shadow on the detection of underground structures.
  • the underground structure detection is performed on the result of the shadow compensation after the sun illumination, which can increase the detection accuracy of the underground structure and reduce the false alarm rate as compared with the underground structure detection of the infrared image without shadow compensation.
  • FIG. 1 is a schematic flow chart of a method for detecting a strip-shaped underground structure based on sun-shadow compensation according to the present invention
  • FIG. 2 is a schematic diagram of latitude and longitude and distance in an embodiment of the present invention.
  • FIG 3 is an elevational elevation diagram of an unprocessed singular point in an embodiment of the present invention.
  • FIG. 5 is a schematic diagram of an altitude rotation process according to an embodiment of the present invention.
  • FIG. 6 is a schematic view showing a solar azimuth angle according to an embodiment of the present invention.
  • Figure 7 is a schematic view showing a solar height angle in an embodiment of the present invention.
  • Figure 8 (a) is a hatched diagram of the area of the simulated tunnel 1 in the embodiment of the present invention.
  • FIG. 8(b) is a view showing a shadow mark expansion of a region of the simulation tunnel 1 in the embodiment of the present invention.
  • Figure 9 (a) is a shadow mark diagram of a region of the simulated tunnel 2 in the embodiment of the present invention.
  • FIG. 9(b) is a diagram showing a shadow mark expansion of a region of the simulated tunnel 2 in the embodiment of the present invention.
  • 10(a) is an infrared remote sensing image of a simulated tunnel 1 region in an embodiment of the present invention
  • 11(a) is an infrared remote sensing image of a simulated tunnel 2 area in an embodiment of the present invention
  • 11(b) is a shadow compensation infrared remote sensing image of a simulated tunnel 2 region in an embodiment of the present invention
  • FIG. 12(a) is a diagram showing detection results of the simulated tunnel 1 without compensation for shading in the embodiment of the present invention.
  • Fig. 12 (b) is a view showing the result of detection of the simulated tunnel 1 after the shading is compensated in the embodiment of the present invention.
  • FIG. 1 The flow of the present invention is shown in Figure 1, wherein the specific implementation method comprises the following steps. Including: obtaining the DEM terrain data step of the specified area, using the DEM and the solar elevation angle, the solar azimuth angle to obtain the image shadow position step, the processing and compensation step of the shadow area, and the detection step of the strip-shaped underground structure after the shadow area correction:
  • the specific process includes the following sub step:
  • the resolution of the simulated infrared remote sensing image is k meters.
  • the resolution of the DEM terrain data needs to be set to k meters.
  • the actual distance d of any two points on the earth satisfies a certain relationship. The relationship is as follows:
  • (x 1 , y 1 ) and (x 2 , y 2 ) represent the latitude and longitude information of these two points, respectively.
  • the programming interface provided by Google Earth Programming firstly locate the latitude and longitude information of each point in Google Earth, and then obtain the altitude information of the anchor point after the positioning is successful. Finally, you will get an elevation map of the same size as the latitude and longitude map, with a resolution of k meters.
  • H(i,j) is the altitude at the point of (i,j)
  • H(i+k,j+l) is the altitude of the adjacent point of (i,j).
  • h is 20.
  • the above steps are performed because the acquired altitude data of the area is not accurate, and there is a certain error in the altitude information, which will adversely affect the subsequent detection. Therefore, it is necessary to make the altitude information as accurate as possible. Due to the randomness of the error, sometimes the altitude is higher than the actual height, and sometimes the altitude is lower than the actual height, so the randomness is used to eliminate this error.
  • the principle is mainly based on the following formula:
  • H i + ⁇ i is the altitude of the reading
  • H i is the actual altitude
  • ⁇ i is the height error
  • the solar azimuth angle ⁇ is determined using the parameters set when the simulated image is generated, and the solar azimuth angle is as shown in Fig. 6, and then the solar azimuth angle ⁇ is rotated for the corrected altitude map.
  • the rotation can be regarded as an affine transformation from one coordinate system to another.
  • the image is rotated by bilinear interpolation.
  • the calculation formula of the affine transformation is as follows:
  • the above-mentioned medium azimuth angle is the shadow of a straight line standing on the ground and the angle between the south and the south.
  • the solar azimuth angle ⁇ 234.65.
  • the shade position is also related to the sun height angle, because the sun height angle determines the proportion of the mountain body and the mountain projection. This ratio directly determines the size of the shadow.
  • the solar height angle is shown in Figure 7. The solar height angle is calculated as:
  • ⁇ represents the local geographic latitude
  • (+/-) is whether the geographic latitude of the requested geographic latitude and the direct sun point is in the same hemisphere: if the symbol in the same hemisphere is -, the symbol in the northern and southern hemispheres is +.
  • the range of the projection can be obtained, so that the size and position of the specific shadow can be determined.
  • tan ⁇ 0.554
  • H(:,i) represents all rows of the i-th column
  • n is the number of columns in the altitude map. Represents the value before and after the exchange of the symbol.
  • the image obtained after the rotation is inverted to the left and right to ensure that the small hills are large.
  • the shadows of the hills are not detected as non-shaded areas, which is equivalent to the sun shining from the east.
  • H'(i,j) represents the altitude map after the rotation is inverted
  • m and n represent the number of rows and columns of the altitude map, respectively.
  • the angle of rotation is the 360-sun azimuth angle ⁇ .
  • the shadow portion of the mountain infrared remote sensing image can be marked according to the marker image.
  • the angle of rotation is 125.35.
  • Fig. 8(a) and Fig. 9(a) respectively, the altitude of the simulated Fig. 1 after the rotation and the altitude hatching map of the simulation Fig. 2 are shown.
  • the shaded area after the expansion is compensated to some extent, and the specific process mainly includes the following substeps:
  • the resulting marker image is expanded to obtain a continuous region.
  • the expansion process uses a structural element to move and compare on the marker map, and processes it according to the result of the comparison.
  • the process of processing is:
  • the point is black, that is, the non-shaded area.
  • the expanded mark pattern can be obtained, and the discontinuous area becomes continuous, which is more practical.
  • the hatched pattern of the simulated image 1 and the hatched pattern of the simulated image 2 are respectively expanded.
  • the shadow region is homogenous to the geological structure of a non-shaded region within a certain range around it, and the infrared characteristics are similar.
  • Fig. 10(a) and Fig. 11(a) there is a significant shadow distribution in the simulated tunnel 1 area and the simulated tunnel 2 area, so after detecting the shadow area, each independent shadow area and its adjacent non-shadow area are performed. Matching, complete the shadow area compensation operation, so as to achieve the purpose of removing the shadow.
  • the algorithm process is as follows:
  • ⁇ noshadow ⁇ p
  • mapping strategy is used to compensate the gray value of the shaded area:
  • F shadow (i,j) is the gray value of the shadow area before compensation
  • F' shadow (i,j) is the gray value of the shadow area after compensation
  • m shadow and ⁇ shadow are the mean values of the gray value of the shadow area.
  • the variance, m noshadow and ⁇ noshadow are the mean and variance of the gray values of the adjacent non-shadow areas
  • A is the compensation intensity coefficient
  • D is the compensation value required for the shadow area.
  • the shaded region is corresponding to the gray value compensation D of each point of the infrared remote sensing image to obtain the final corrected grayscale image.
  • the infrared remote sensing image of the simulated tunnel 1 area and the simulated tunnel 2 area after shadow compensation can be seen from the obtained infrared remote sensing image, and the shadow is eliminated to some extent. .
  • the specific process includes the following sub-steps:
  • the infrared remote sensing image after shadow compensation is used to detect the strip-shaped underground structure, and the thermal infrared remote sensing image is detected by pixel-by-pixel rolling segmentation according to the approximate position and direction information of the tunnel.
  • the length of each segment is 250 meters and the width is 30 meters (ie 25
  • the pixel points are 3 pixels wide, and the left and right traversals of the pixels are detected to detect whether there is a band-like underground structure, and finally the preliminary tunnel recognition result is obtained.
  • the same method is used to detect the underground structure in the direction perpendicular to the tunnel, and the false alarm rate is calculated.
  • the following comparison table is obtained. It can be seen that after the shadow is removed, the recognition rate is improved and the false alarm rate is also reduced.
  • the table is as follows:

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Geophysics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Geometry (AREA)
  • Image Processing (AREA)

Abstract

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

Description

一种基于太阳照射阴影补偿的带状地下结构探测方法 [技术领域]
本发明属于遥感技术、图像处理和模式识别的交叉领域，具体涉及一种基于太阳照射阴影补偿的带状地下结构探测方法。
[背景技术]
阴影是在自然界中普遍存在的一种物理现象，是由于光源被物体遮挡而产生的。图像中阴影的存在对计算机视觉领域的相关问题有不同的影响，这种影响有有利的也有不利的，比如在虚拟现实、3D游戏中为物体添加阴影，可以提高场景的真实感。但是更多的时候，图像中的阴影会对计算机视觉的相关问题产生不利的影响，如在遥感图像中，阴影的存在会影响后继的图像匹配、模式识别和地物的提取等多种遥感图像处理操作；在视频监控中，阴影和运动目标结合在一起，导致计算机对目标物体的提取和追踪出现错误。因此，有必要对图像中的阴影进行检测和分析，并根据需要进而消除或减弱阴影的影响。
目前，国内外很多专家学者对视频中的阴影去除进行了比较深入的研究，提出了很多有效的算法。在视频监控和车辆追踪等领域得到了广泛的应用，但是无法处理静态的图像。在静态图像中的阴影检测方法，一直是一个比较困难的问题，因为静态图像所包含的信息量更少，不能像在视频中一样利用帧间关系来进行阴影检测。同时，在遥感成像等基于静态图像的研究领域，在检测出阴影之后，还要尽可能的减弱或消除阴影的影响，恢复阴影中物体的本来面目，从而提高红外遥感图像的真实性。
目前国内外并未有使用红外成像技术探测带状地下结构的相关报 道，比如暗河，地下隧道等带状地下结构。红外遥感图像中阴影的存在增加了对带状地下结构探测的虚警率，降低了探测率。
发明内容
本发明提出了一种基于太阳照射阴影补偿的带状地下结构探测方法，结合数字高度模型(Digital Elevation Model，DEM)对单幅红外遥感图像进行阴影检测，并对阴影所在的区域的灰度进行了补偿，从而减弱了阴影对探测的影响。
本发明提供的一种基于太阳照射阴影补偿的带状地下结构探测方法，具体步骤如下：
(1)根据区域的地理位置得到区域的经纬度信息，然后通过图像的分辨率计算经纬度步长，获取指定区域的海拔高度信息，最后对所得到的海拔高度信息进行校正处理，具体过程包括以下子步骤：
(1.1)利用所检测的区域确定经纬度步骤
(1.1.1)经纬度信息的确定可以根据需要检测的位置和所需数据的分辨率来确定，待检测区域左上方点经纬度信息为A(x1，y1)，右下方点的经纬度信息为B(x2，y2)。
(1.1.2)根据仿真的红外遥感图像的设置参数可知红外遥感图像的分辨率为k米，因此需要将DEM地形数据的分辨率也设置为k米。地球上任意两点的实际距离D满足一定的关系，关系式如下：
C＝sin(y1)*sin(y2)*cos(x1-x2)+cos(y1)*cos(y2)
D＝R*arccos(C)*Pi/180
R＝6371.004
上式中(x1，y1)和(x2，y2)分别代表这两点的经纬度信息。
但是，我们只需求平行于赤道和垂直于赤道方向的两点距离，故平 行于赤道时的距离d1和经纬度的关系如下：
C＝sin(y1)2*cos(x1-x2)+cos(y1)2 (y1＝y2)
d1＝R*arccos(C)*Pi/180
R＝6371.004
lon＝k*(x2-x1)/d1
式中，d1为平行于赤道的两点距离，lon为经度步长。
垂直于赤道时的距离d2和经纬度的关系如下：
C＝cos(y1-y2) (x1＝x2)
d2＝R*arccos(C)*Pi/180
R＝6371.004
lat＝k*(y1-y2)/d2
式中，d2为垂直于赤道的两点距离，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大小的模板对海拔高度图进行均值滤波处理。
进行以上的步骤是由于所获取到的区域海拔数据并不是准确无误的，海拔高度信息会存在一定的误差，会给以后的检测带来不利的影响，故需要尽量使海拔高度信息准确。由于误差的出现具有随机性，有时海拔高度比实际的高，有时海拔高度比实际的低，所以利用其随机性消除这种误差，其原理主要是基于以下的公式：
Figure PCTCN2015072663-appb-000001
式中Hi+εi为读取的海拔高度，Hi为实际的海拔高度，εi为高度误差。从上式中可以看出，由于误差的随机性，平均误差
Figure PCTCN2015072663-appb-000002
最后得到比较准确的区域海拔高度信息数据。
(2)利用太阳方位角等图像信息对校正后的海拔高度图进行旋转处理，通过太阳高度角信息计算山体本体和山体投影的比例，确定阴影的大小和位置，再经过一系列的处理得基于DEM的图像阴影位置，具体过程主要包括以下子步骤：
(2.1)利用太阳方位角旋转海拔高度图步骤
利用产生仿真图像时设置的参数确定太阳方位角θ，然后对校正后的海拔高度图旋转太阳方位角θ，旋转可以看做是从一个坐标系变换到了另一个坐标系的仿射变换，为了减小旋转带来的误差，采用双线性插值的方法对图像进行旋转处理。仿射变换的计算公式如下：
Figure PCTCN2015072663-appb-000003
Figure PCTCN2015072663-appb-000004
最后得到旋转后的海拔高度图，超出的部分海拔高度设置为0。
上述中太阳方位角为竖立在地面上的直线的阴影和正南方的夹角。
(2.2)计算山体本体和山体投影的比例，确定阴影的大小和位置
其中阴影位置和太阳高度角也有关系，因为太阳高度角决定山体本体和山体投影的比例，这个比例直接决定阴影的大小，而太阳高度角的计算公式为：
Figure PCTCN2015072663-appb-000005
φ代表当地的地理纬度，
Figure PCTCN2015072663-appb-000006
代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+。
山体本体和山体阴影的关系为：
Figure PCTCN2015072663-appb-000007
由上式可以得到投影的范围大小，从而可以确定具体阴影的大小和位置。
(2.3)阴影的区域的检测标记步骤
(2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：
Figure PCTCN2015072663-appb-000008
式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，
Figure PCTCN2015072663-appb-000009
代表交换该符号前后的数值。
将旋转后得到的图进行左右倒置处理，是为了保证小的山丘在大的山丘的阴影部位的时候不被检测为非阴影区域，从而相当于太阳是从东边照射过来的。
(2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，为了对其中的每一个点进行判断，判断该点是否存在于阴影区，判断的公式如下：
H'(i,j)≠0 (i＝1,2,3...m、j＝1,2,3...n)
Figure PCTCN2015072663-appb-000010
式中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)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：
Figure PCTCN2015072663-appb-000011
D＝F'shadow(i,j)-Fshadow(i,j)
式中Fshadow(i,j)是补偿之前的阴影区域灰度值，F'shadow(i,j)是补偿之 后阴影区域灰度值，mshadow和σshadow是阴影区域灰度值的均值和方差，mnoshadow和σ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(x1，y1)＝(116.150049，40.296833)，右下方点的经纬度信息为
B(x2，y2)＝(116.194775，40.260787)。
(1.1.2)所仿真的红外遥感图像的分辨率为k米，需要将DEM地形数据的分辨率也设置为k米，地球上任意两点的实际距离d满足一定的关系，关系式如下：
C＝sin(y1)*sin(y2)*cos(x1-x2)+cos(y1)*cos(y2)
D＝R*arccos(C)*Pi/180
R＝6371.004
上式中(x1，y1)和(x2，y2)分别代表这两点的经纬度信息。
如图2所示，故平行于赤道时的距离d1和经纬度的关系如下：
C＝sin(y1)2*cos(x1-x2)+cos(y1)2(y1＝y2)
d1＝R*arccos(C)*Pi/180
R＝6371.004
lon＝k*(x2-x1)/d1
本实例中，d1＝3800m，lon＝0.0001177。
垂直于赤道时的距离d2和经纬度的关系如下：
C＝cos(y1-y2)(x1＝x2)
d2＝R*arccos(C)*Pi/180
R＝6371.004
lat＝k*(y1-y2)/d2
本实例中，d2＝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大小的模板对海拔高度图进行均值滤波处理。
进行以上的步骤是由于所获取到的区域海拔数据并不是准确无误的，海拔高度信息会存在一定的误差，会给以后的检测带来不利的影响，故需要尽量使海拔高度信息准确。由于误差的出现具有随机性，有时海拔高度比实际的高，有时海拔高度比实际的低，所以利用其随机性消除这种误差，其原理主要是基于以下的公式：
Figure PCTCN2015072663-appb-000012
式中Hi+εi为读取的海拔高度，Hi为实际的海拔高度，εi为高度误差。从上式中可以看出，由于误差的随机性，平均误差
Figure PCTCN2015072663-appb-000013
最后得到比较准确的区域海拔高度信息数据。本实例中，m为5。
(2)利用太阳方位角等图像信息对校正后的海拔高度图进行旋转处理，通过太阳高度角等图像信息计算山体本体和山体阴影的比例，再经过一系列的处理得基于DEM的图像阴影位置，具体过程主要包括以下子步骤：
(2.1)利用太阳方位角旋转海拔高度图步骤
利用产生仿真图像时设置的参数确定太阳方位角θ，太阳方位角如图6所示，然后对校正后的海拔高度图旋转太阳方位角θ。如图5所示，旋转可以看做是从一个坐标系变换到了另一个坐标系的仿射变换，为了减小旋转带来的误差，采用双线性插值的方法对图像进行旋转处理。仿射变换的计算公式如下：
Figure PCTCN2015072663-appb-000014
Figure PCTCN2015072663-appb-000015
最后得到旋转后的海拔高度图，超出的部分海拔高度设置为0。
上述中太阳方位角为竖立在地面上的直线的阴影和正南方的夹角。 本实例中，太阳方位角θ＝234.65。
(2.2)计算山体本体和山体投影的比例步骤
其中阴影位置和太阳高度角也有关系，因为太阳高度角决定山体本体和山体投影的比例，这个比例直接决定阴影的大小，太阳高度角如图7所示，太阳高度角的计算公式为：
Figure PCTCN2015072663-appb-000016
φ代表当地的地理纬度，
Figure PCTCN2015072663-appb-000017
代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+。本实例中，地理纬度和太阳直射点纬度在同一个半球，得β＝29。
山体本体和山体阴影的关系为：
Figure PCTCN2015072663-appb-000018
由上式可以得到投影的范围大小，从而可以确定具体阴影的大小和位置。本实例中，tanβ＝0.554
(2.3)阴影的区域的检测标记步骤
(2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：
Figure PCTCN2015072663-appb-000019
式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，
Figure PCTCN2015072663-appb-000020
代表交换该符号前后的数值。
将旋转后得到的图进行左右倒置处理，是为了保证小的山丘在大的 山丘的阴影部位的时候不被检测为非阴影区域，从而相当于太阳是从东边照射过来的。
(2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，为了对其中的每一个点进行判断，判断该点是否存在于阴影区，判断的公式如下：
H'(i,j)≠0 (i＝1,2,3...m、j＝1,2,3...n)
Figure PCTCN2015072663-appb-000021
式中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)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：
Figure PCTCN2015072663-appb-000022
D＝F'shadow(i,j)-Fshadow(i,j)
式中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。
从上述中得到需要补偿的灰度值后，遍历所有检测到的阴影区域，将阴影区域对应到红外遥感图像每一个点的灰度值补偿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)所示，通过对两幅图的对比可以看到对阴影补偿后的图进行带状地下结构检测结果比未补偿的检测结果要好。
通过对为去除阴影和去除阴影的红外遥感图像进行探测识别，得到了以下对比表，可以看出去除阴影后，识别率有所提高，同时虚警率也有所降低。表格如下：
Figure PCTCN2015072663-appb-000023
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

  1. 一种基于太阳照射阴影补偿的带状地下结构探测方法，其特征在于，所述方法包括如下步骤：
    (1)根据待检测区域的地理位置确定待检测区域的经纬度信息，然后根据红外遥感图像的分辨率计算经纬度步长，并根据待检测区域的经纬度信息和经纬度步长获取待检测区域的海拔高度信息，最后对所得到的海拔高度信息进行校正处理；
    (2)利用太阳方位角信息对校正后的海拔高度图进行旋转处理，通过太阳高度角信息计算山体本体和山体阴影的比例，确定海拔高度图中阴影的大小和位置，对旋转后的海拔高度图进行处理得到图像阴影位置的标记图；
    (3)对得到的图像阴影位置的标记图进行膨胀处理，并对膨胀后的阴影区域进行补偿；
    (4)对阴影区域补偿后的红外遥感图像进行初步的探测识别，然后利用位置的关联性进行二次的判断，从而探测得到带状地下结构。
  2. 如权利要求1所述的方法，其特征在于，所述步骤(1)中根据待检测区域的地理位置确定待检测区域的经纬度信息，然后根据红外遥感图像的分辨率计算经纬度步长，具体包括如下子步骤：
    (1.1.1)根据待检测区域的地理位置获得待检测区域左上方点经纬度信息为A(x1，y1)，以及右下方点的经纬度信息为B(x2，y2)；
    (1.1.2)获得红外遥感图像的分辨率为k米，分别根据下式求取经度步长lon和纬度步长lat，
    C＝sin(y1)2*cos(x1-x2)+cos(y1)2 (y1＝y2)
    d1＝R*arccos(C)*Pi/180
    R＝6371.004
    lon＝k*(x2-x1)/d1
    C＝cos(y1-y2) (x1＝x2)
    d2＝R*arccos(C)*Pi/180
    R＝6371.004   。
    lat＝k*(y1-y2)/d2
  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)中利用太阳方位角信息对校正后的海拔高度图进行旋转处理，具体包括：
    确定太阳方位角θ，对校正后的海拔高度图旋转太阳方位角θ，变换公式如下：
    Figure PCTCN2015072663-appb-100001
    Figure PCTCN2015072663-appb-100002
    对旋转后得到的海拔高度图，超出的部分海拔高度设置为0；
    其中，所述中太阳方位角为地面的竖立直线的阴影和正南方的夹角。
  6. 如权利要求1或2所述的方法，其特征在于，所述步骤(2)中通过太阳高度角信息计算山体本体和山体阴影的比例，确定海拔高度图中阴影的大小和位置，具体包括：
    根据公式
    Figure PCTCN2015072663-appb-100003
    计算太阳高度角，其中φ代表当地的地理纬度，
    Figure PCTCN2015072663-appb-100004
    代表太阳直射点地理纬度，(+/-)是所求地理纬度和太阳直射点的地理纬度是否在同一个半球：如果在同一半球符号为-，在南北半球符号为+；
    Figure PCTCN2015072663-appb-100005
    求得阴影的大小和位置。
  7. 如权利要求1或2所述的方法，其特征在于，所述步骤(2) 中对旋转后的海拔高度图进行处理得到图像阴影位置的标记图，具体包括：
    (2.3.1)将旋转后得到的海拔高度图左右倒置，具体地可以根据下述公式从第一列遍历到中间列对旋转后得到的图进行处理实现倒置，遍历的过程如下：
    Figure PCTCN2015072663-appb-100006
    式中H(:，i)代表第i列的所有行，n为海拔高度图的列数，
    Figure PCTCN2015072663-appb-100007
    代表交换该符号前后的数值；
    (2.3.2)对左右倒置后的海拔高度图进行从左到右、从上到下的全图遍历，判断每一个点是否存在于阴影区，判断公式如下：
    H'(i,j)≠0 (i＝1,2,3...m、j＝1,2,3...n)
    Figure PCTCN2015072663-appb-100008
     (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)在得到阴影区域和其邻近的非阴影区域之后，采用如下的映射策略对阴影区域的灰度值进行补偿：
    Figure PCTCN2015072663-appb-100009
    D＝F'shadow(i,j)-Fshadow(i,j)
    式中Fshadow(i,j)是补偿之前的阴影区域灰度值，F'shadow(i,j)是补偿之后阴影区域灰度值，mshadow和σshadow是阴影区域灰度值的均值和方差，mnoshadow和σnoshadow是邻近非阴影区域灰度值的均值和方差，A为补偿强度系数，D为阴影区域需要的补偿值；
    得到需要补偿的灰度值D后，遍历所有检测到的阴影区域，将阴影区域对应到红外遥感图像上每一个点的灰度值补偿D，得到最终补偿后的红外遥感图像。
  10. 如权利要求1或2所述的方法，其特征在于，所述步骤(4)具体包括：
    (4.1)初步判断带状地下结构位置
    利用阴影补偿后的红外遥感图像探测带状地下结构，依据隧道所在大致位置与方向信息对热红外遥感图像进行逐像素滚动分段探测识别，对像素点的左右遍历，探测是否存在带状地下结构，最后得到初步的隧道识别结果；同时在垂直于隧道的方向利用相同的方法探测地下结构，计算虚警率；
    (4.2)根据位置关联性进行二次判断
    对初步识别出的符合脉冲模式的分散段进行统计分析，依据预设长度上各个分散段出现的概率与位置分散程度进行判断可否忽视掉分离相邻的符合脉冲模式段的小空白：若根据判断可以忽视掉这些分散脉冲模式段中间的小空白，那么该预设长度段就确认为带状目标所在位置；若根据判断不可以忽视掉这些分散脉冲模式段中间的小空白，那么该预设长度段就确认为非带状目标所在位置。
PCT/CN2015/072663 2014-12-30 2015-02-10 一种基于太阳照射阴影补偿的带状地下结构探测方法 WO2016106950A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US15/106,686 US9582885B2 (en) 2014-12-30 2015-02-10 Zonal underground structure detection method based on sun shadow compensation

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2014108445785 2014-12-30
CN201410844578.5A CN104637073B (zh) 2014-12-30 2014-12-30 一种基于太阳照射阴影补偿的带状地下结构探测方法

Publications (1)

Publication Number Publication Date
WO2016106950A1 true WO2016106950A1 (zh) 2016-07-07

Family

ID=53215779

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2015/072663 WO2016106950A1 (zh) 2014-12-30 2015-02-10 一种基于太阳照射阴影补偿的带状地下结构探测方法

Country Status (3)

Country Link
US (1) US9582885B2 (zh)
CN (1) CN104637073B (zh)
WO (1) WO2016106950A1 (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109934788A (zh) * 2019-03-22 2019-06-25 鲁东大学 一种基于标准遥感图像的遥感图像缺失数据修复方法
CN111797679A (zh) * 2020-05-19 2020-10-20 中国地质大学（武汉） 一种遥感纹理信息处理方法、装置、终端及存储介质
CN112200858A (zh) * 2020-10-10 2021-01-08 长光卫星技术有限公司 基于高分辨率光学遥感图像的外浮顶油罐储量分析方法
CN112668441A (zh) * 2020-12-24 2021-04-16 中国电子科技集团公司第二十八研究所 一种结合先验知识的卫星遥感影像飞机目标识别方法
CN114323308A (zh) * 2022-01-06 2022-04-12 中国地质大学（北京） 一种遥感城市地表温度角度效应的校正方法
CN117671167A (zh) * 2023-10-19 2024-03-08 兰州交通大学 一种基于山体阴影分析的启发式dem综合方法
US12106024B2 (en) * 2022-01-27 2024-10-01 Huazhong University Of Science And Technology Geologically constrained infrared imaging detection method and system for urban deeply-buried strip-like passage

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104501959B (zh) * 2014-12-30 2016-08-17 华中科技大学 一种红外图谱关联智能探测方法及装置
CN105426881B (zh) * 2015-12-24 2017-04-12 华中科技大学 山体背景热场模型约束的地下热源昼间遥感探测定位方法
CN105844613A (zh) * 2016-03-14 2016-08-10 昆明理工大学 一种识别天文图像太阳黑子中半影纤维亮点的方法
CN107564017B (zh) * 2017-08-29 2020-01-10 南京信息工程大学 一种城市高分遥感影像阴影检测及分割方法
CN109671111B (zh) * 2018-11-09 2023-07-07 西安电子科技大学 基于可见光遥感图像的温度场调制方法
US11243329B2 (en) * 2019-03-01 2022-02-08 Nancy Kerr Del Grande Detecting subsurface objects and voids using thermal inertia
CN110210453B (zh) * 2019-06-14 2021-06-29 中国资源卫星应用中心 一种基于遥感图像特征的油罐存储量确定方法及系统
CN110969088B (zh) * 2019-11-01 2023-07-25 华东师范大学 一种基于显著性检测与深度孪生神经网络的遥感影像变化检测方法
CN111223172B (zh) * 2020-01-10 2023-06-30 成都中科合迅科技有限公司 一种雷达探测范围绘制和融合方法、设备及存储介质
CN111710016B (zh) * 2020-06-19 2023-07-11 北京世冠金洋科技发展有限公司 一种晨昏线图层获取方法及装置
CN112964643B (zh) * 2021-02-03 2022-04-19 福州大学 一种遥感影像可见光波段地形落影校正方法
CN113534851A (zh) * 2021-06-10 2021-10-22 深圳市爱深盈通信息技术有限公司 控制遮阳伞的方法
CN113870237B (zh) * 2021-10-09 2024-03-08 西北工业大学 一种基于水平扩散的复合材料图像阴影检测方法
CN115546648B (zh) * 2022-10-24 2023-05-23 南京师范大学 基于光线追踪的dem冲积扇提取方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006011832A (ja) * 2004-06-25 2006-01-12 Photo System:Kk 空撮画像の解析作業支援装置とそのための方法とプログラム。
CN1845088A (zh) * 2006-04-28 2006-10-11 上海大学 起伏地面上太阳直射光归一化分布图像的计算机生成方法
CN101030301A (zh) * 2007-03-29 2007-09-05 上海大学 起伏地面上遥感分布图像的计算机虚拟
CN101034477A (zh) * 2007-03-29 2007-09-12 上海大学 遥感数字图像上阴影消除和阴影中像元遥感值恢复方法
CN101034475A (zh) * 2007-03-29 2007-09-12 上海大学 无阴影卫星遥感正射数字图像的计算机生成方法
CN102622738A (zh) * 2012-03-08 2012-08-01 北京师范大学 一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008100308A1 (en) * 2006-12-19 2008-08-21 Weil Gary J Remote sensing of subsurface artifacts by use of visual and thermal imagery
CN101458328A (zh) * 2009-01-04 2009-06-17 上海大学 月球表面反射率遥感反演方法
JP5632173B2 (ja) * 2010-03-10 2014-11-26 一般財団法人宇宙システム開発利用推進機構 Ｓａｒデータ処理方法及びｓａｒデータ処理システム

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006011832A (ja) * 2004-06-25 2006-01-12 Photo System:Kk 空撮画像の解析作業支援装置とそのための方法とプログラム。
CN1845088A (zh) * 2006-04-28 2006-10-11 上海大学 起伏地面上太阳直射光归一化分布图像的计算机生成方法
CN101030301A (zh) * 2007-03-29 2007-09-05 上海大学 起伏地面上遥感分布图像的计算机虚拟
CN101034477A (zh) * 2007-03-29 2007-09-12 上海大学 遥感数字图像上阴影消除和阴影中像元遥感值恢复方法
CN101034475A (zh) * 2007-03-29 2007-09-12 上海大学 无阴影卫星遥感正射数字图像的计算机生成方法
CN102622738A (zh) * 2012-03-08 2012-08-01 北京师范大学 一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109934788A (zh) * 2019-03-22 2019-06-25 鲁东大学 一种基于标准遥感图像的遥感图像缺失数据修复方法
CN111797679A (zh) * 2020-05-19 2020-10-20 中国地质大学（武汉） 一种遥感纹理信息处理方法、装置、终端及存储介质
CN112200858A (zh) * 2020-10-10 2021-01-08 长光卫星技术有限公司 基于高分辨率光学遥感图像的外浮顶油罐储量分析方法
CN112668441A (zh) * 2020-12-24 2021-04-16 中国电子科技集团公司第二十八研究所 一种结合先验知识的卫星遥感影像飞机目标识别方法
CN112668441B (zh) * 2020-12-24 2022-09-23 中国电子科技集团公司第二十八研究所 一种结合先验知识的卫星遥感影像飞机目标识别方法
CN114323308A (zh) * 2022-01-06 2022-04-12 中国地质大学（北京） 一种遥感城市地表温度角度效应的校正方法
CN114323308B (zh) * 2022-01-06 2023-05-02 中国地质大学（北京） 一种遥感城市地表温度角度效应的校正方法
US12106024B2 (en) * 2022-01-27 2024-10-01 Huazhong University Of Science And Technology Geologically constrained infrared imaging detection method and system for urban deeply-buried strip-like passage
CN117671167A (zh) * 2023-10-19 2024-03-08 兰州交通大学 一种基于山体阴影分析的启发式dem综合方法
CN117671167B (zh) * 2023-10-19 2024-05-28 兰州交通大学 一种基于山体阴影分析的启发式dem综合方法

Also Published As

Publication number Publication date
CN104637073A (zh) 2015-05-20
US9582885B2 (en) 2017-02-28
CN104637073B (zh) 2017-09-15
US20160371841A1 (en) 2016-12-22

Similar Documents

Publication Publication Date Title
WO2016106950A1 (zh) 一种基于太阳照射阴影补偿的带状地下结构探测方法
Huang et al. An automatic change detection method for monitoring newly constructed building areas using time-series multi-view high-resolution optical satellite images
CN107093205B (zh) 一种基于无人机图像的三维空间建筑物窗户检测重建方法
CN111462200A (zh) 一种跨视频行人定位追踪方法、系统及设备
CN107610164B (zh) 一种基于多特征混合的高分四号影像配准方法
US10366501B2 (en) Method and apparatus for performing background image registration
Gance et al. Target Detection and Tracking of moving objects for characterizing landslide displacements from time-lapse terrestrial optical images
CN110569797B (zh) 地球静止轨道卫星影像山火检测方法、系统及其存储介质
US11132573B2 (en) Determining compass orientation of imagery
CN109741446B (zh) 一种三维数字地球动态生成精细海岸地形的方法
CN106096497B (zh) 一种针对多元遥感数据的房屋矢量化方法
CN103871039A (zh) 一种sar图像变化检测差异图生成方法
Gao et al. A general deep learning based framework for 3D reconstruction from multi-view stereo satellite images
Zhang et al. Lidar-guided stereo matching with a spatial consistency constraint
CN115937673B (zh) 一种基于移动终端照片的地理要素快速变化发现方法
Cheng et al. Semi-automatic road centerline extraction in high-resolution SAR images based on circular template matching
CN116879870A (zh) 一种适用于低线束3d激光雷达的动态障碍物去除方法
Wang Automatic extraction of building outline from high resolution aerial imagery
Zhou et al. Building occlusion detection from ghost images
Parmehr et al. Automatic registration of optical imagery with 3d lidar data using local combined mutual information
CN112489207B (zh) 一种空间约束的密集匹配点云平面基元提取方法
CN117367404A (zh) 基于动态场景下slam的视觉定位建图方法及系统
CN105631849A (zh) 多边形目标的变化检测方法及装置
CN114926524B (zh) 一种用于提高烟幕红外有效遮蔽面积测量精度的方法
CN106204596A (zh) 一种基于高斯拟合函数与模糊混合估计的全色波段遥感影像云检测方法

Legal Events

Date Code Title Description
WWE Wipo information: entry into national phase

Ref document number: 15106686

Country of ref document: US

121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 15874643

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 15874643

Country of ref document: EP

Kind code of ref document: A1