Estimating the Volume of Oil Tanks Based on High-Resolution Remote Sensing Images

The purpose of this study is to obtain oil tank volumes from high-resolution satellite imagery to meet the need to measure oil tank volume globally. A preprocessed remote sensing HSV image is used to extract the shadow of the oil tank by Otsu thresholding, shadow area thresholding, and morphological closing. The oil tank shadow is crescent-shaped. Hence, a median method based on sub-pixel subdivision positioning is used to calculate the shadow length of the oil tank and then determine its height with high precision. The top of the tank and its radius in the image are identified using the Hough transform. The final tank volume is calculated using its height and radius. A high-resolution Gaofen 2 optical remote sensing image is used to evaluate the proposed method. The actual height and volume of the tank we tested were 21.8 m and 109,532 m3. The experimental results show that the mean absolute error of the height of the tank calculated by the median method is 0.238 m, the relative error is within 1.15%, and the RMES is 0.23. The result is better than the previous work. The absolute error between the calculated and the actual tank volumes ranges between 416 and 3050 m3, and the relative error ranges between 0.38% and 2.78%. These results indicate that the proposed method can calculate the volume of oil tanks with high precision and sufficient accuracy for practical applications.


Introduction
For many years now, the status and production of oil have been highly valued.In an oil strategy, how to estimate the current amount of global oil reserves and predict changes in production are key for determining oil procurement and transportation measures.They are also important factors in economic development.If the volume of an oil tank can be calculated, it is possible to estimate the storage capacity of an oil storage base.In recent years, with the development and improvement of satellite remote sensing technology, high-resolution image information from any part of the world can be obtained easily [1].Since 1989, in aerial photogrammetry, researchers have long used shadow information to estimate building height [2].However, most of the buildings considered in the literature are houses, and their shadows are approximately rectangular.The shape of a tank shadow in a remote sensing image is similar to a crescent shape.There are, hence, few methods to calculate the actual length of a tank shadow.
The HSV color space was first proposed by Smith in 1978.In the HSV color space, hue, saturation, and value are used to express the color of the image [3] since the geography of tank areas, deep blue water, and dark green vegetation areas cannot be observed in oil tank images.This study hence detects shadows in the HSV color space.Shadows have low brightness, high saturation, and increased hue in the HSV color space.To take advantage of this characteristic, a hue-to-value ratio image can be obtained, thereby enhancing the low luminance and hue values of the shadow region [4].Using this image, the method proposed in this paper performs a series of segmentation and optimization algorithms to obtain the tank shadow.
The sub-pixel edge detection method can improve the accuracy of edge detection [5].The proposed method uses the idea of sub-pixel subdivision positioning to improve the accuracy of shadow length calculation.The tank height is calculated based on the geometric relationship between the length of the shadow and the actual height of the tank.Finally, the tank volume is calculated using the tank radius obtained by the Hough transform algorithm.
The rest of this paper is organized as follows: Section 2 presents a summary of work related to shadow extraction and building height calculation.Section 3 details the proposed approach for calculating tank volume based on high-resolution optical remote sensing.The image preprocessing, oil tank shadow extraction, and oil tank volume calculation are described in detail.The experimental results are presented and discussed in Section 4. Finally, the paper concludes in Section 5 with some recommendations for further work.

Remote Sensing Image Shadow Detection
Many approaches for shadow detection have been reported in the literature; they are often categorized into property-based and model-based methods [6].However, the model-based algorithms have strong disadvantages [7].The core of the algorithm is based on the geometry of the occlusion object, digital surface model data and solar illumination orientation, sensors and other related parameters.However, in the practical research, this information is difficult to obtain and has limitations.
Property-based methods take into account shadow properties that are generally directly deduced from the image data information.This includes both radiometric attributes with spectral features and textural attributes through spatial features.Some examples of the methods in this class are simple histogram thresholding and invariant color models [8].Histogram thresholding is the simplest shadow detection tool, but is also the most popular in the literature.For example, Nagao et al. [9] used a linear combination of RGB and NIR channels.Then, the search for a threshold in the histogram is based on various assumptions about the histogram such as a bimodal distribution [10], an underlying Gaussian mixture model [11], or the number of peaks and valleys in the histogram [12].The threshold can also be chosen arbitrarily by visual inspection [13].The main drawback of these methods is the misclassification of sunlit dark objects (such as water) as shadows as well as the misclassification of shadowed bright objects as sunlit regions.
The shadow chromaticity expressed in a 2D image is exploited in invariant color model methods.Shadows have high hue value, high saturation in blue-violet channels, and low intensity.Polidorio et al. [14] subtracted the luminance I component from the saturation component S in the HSI color space, separating the shaded and non-shaded regions.Tsai [4] developed a spectral ratio between hue and intensity, and showed that his method outperforms those of Polidorio et al. [14] and Huang et al. [15].In contrast, Chung et al. [16] improved the spectral ratio of Tsai and showed that their method has better accuracy.Luo et al. [17] based their method on the properties of each component of the HSI color space combined with the maximum inter-class variance threshold method for the separation of shadow regions and non-shadow regions.

Shadow-Based Building Height Calculation
In the past 20 years, the technique of inferring the height of buildings based on the shadows of buildings in remote sensing images has become relatively mature [2].Cheng et al. established a simple calculation model of building shadows and their heights, and obtained a mean-square error of 3.69 m over 42 building height measurements [18].The average height of 42 buildings in [18] is about 20 m, so the result is not so good.Hartl and Cheng [19] extended the work of Cheng and Thiel.They tested 77 buildings and obtained a root-mean-square error of 6.19 m.Massalabi [20] provided a formula for the determination of building heights based on shadow length, solar elevation, solar azimuth, satellite elevation, and satellite azimuth to make calculation easier.Izadi [21] and Wang [22] respectively deduced formulas for building height based on the sine-cosine theorem using shadow length, solar elevation, satellite elevation, and the difference between solar azimuth and satellite azimuth.To determine building heights based only on single satellite images, Qi [23] presented a new method named the corner-shadow-length ratio (CSLR) method.Cheng and Thiel [18] and Turker and Sumer [24] used the building shadows to calculate the height of buildings that collapsed in an earthquake.More recently, Qi [25] estimated building height from Google Earth images by first calculating the ratio of building height to the shadow length of known buildings, and thereafter utilizing the identified shadow-length ratio to obtain the heights of other buildings with unknown heights.Liasis and Stavrou [26] developed a custom filter for enhancing shadows and reducing the spectral heterogeneity of the regions of interest (ROIs) to form an optimized contour model for estimating building height using a shadow length and a solar elevation angle.Liu [27].vectorized the extracted shadows and used parallel lines to measure the actual length of the building's shadows, thereby calculating the height of the building.Zhang [28] extracted the shadow area by multi-band spectral difference, calculated the shadow length using the pixel number method, and then found the height of the building.

Overview of the Approach
A flowchart of the proposed approach is illustrated in Figure 1.Image preprocessing is performed on the original Gaofen 2 satellite image data.The shadow area of the oil tank is obtained by converting the image to HSV color space, using the maximum inter-class variance, a threshold, and morphological operators.The actual length of the tank shadow is calculated using a median method based on sub-pixel subdivision positioning.The height of the tank is determined based on the geometric relationship between the length of the shadow and the height of the tank.Then, the Hough transform is used to identify the top of the tank and calculate its radius.Finally, the volume of the tank is obtained according to the formula for the volume of a cylinder.

Image Preprocessing
To ensure the accuracy of shadow extraction, the original remote sensing data image need to be preprocessed.The main processes in the proposed method include radiation correction, geometric correction, image fusion, and image cropping.
Radiation correction is divided into two key steps: radiometric calibration and atmospheric correction.Radiation calibration is to establish a quantitative relationship between the brightness value of the image pixel and the actual radiance value.Atmospheric correction can eliminate the influence of atmospheric and illumination factors on the reflection of ground objects, so physical parameters, such as reflectivity, emissivity, and surface temperature of the object, can be obtained.The original remote sensing image data of Gaofen 2 has been processed by on-orbit radiation calibration [29][30][31][32].In the proposed method, atmospheric correction is performed using the radiation transfer model, which is more accurate.We performed general geometric correction.Pan-sharpening fusion is employed to fuse the images.As an example, the 0.81-m panchromatic image is shown in Figure 2a

Image Preprocessing
To ensure the accuracy of shadow extraction, the original remote sensing data image need to be preprocessed.The main processes in the proposed method include radiation correction, geometric correction, image fusion, and image cropping.
Radiation correction is divided into two key steps: radiometric calibration and atmospheric correction.Radiation calibration is to establish a quantitative relationship between the brightness value of the image pixel and the actual radiance value.Atmospheric correction can eliminate the influence of atmospheric and illumination factors on the reflection of ground objects, so physical parameters, such as reflectivity, emissivity, and surface temperature of the object, can be obtained.The original remote sensing image data of Gaofen 2 has been processed by on-orbit radiation calibration [29][30][31][32].In the proposed method, atmospheric correction is performed using the radiation transfer model, which is more accurate.We performed general geometric correction.Pan-sharpening fusion is employed to fuse the images.As an example, the 0.81-m panchromatic image is shown in Figure 2a and the 3.2-m multispectral image is shown in Figure 2b.The merged image is shown in Figure 3.The remote sensing images were cropped to extract the areas of interest for the experiments in this paper.The cropped images are shown in Figure 4.

Oil Tank Shadow Extraction
In the HSV color space, hue, saturation, and value are used to express the color of the image.In shadow areas, the hue values are relatively consistent and high in value, whereas non-shadow areas have relatively low hue values, as shown in Figure 5.In contrast, in Figure 6, the saturation of the shaded area and the saturation of the surrounding ground are difficult to distinguish.The value of the shadow area is smaller than that of the non-shadow area.We use the (H + 1)/(V + 1) ratio image to enhance the low luminance and hue values of the shadow region, where H is the hue channel image and V is the value channel image.

Oil Tank Shadow Extraction
In the HSV color space, hue, saturation, and value are used to express the color of the image.In shadow areas, the hue values are relatively consistent and high in value, whereas non-shadow areas have relatively low hue values, as shown in Figure 5.In contrast, in Figure 6, the saturation of the shaded area and the saturation of the surrounding ground are difficult to distinguish.The value of the shadow area is smaller than that of the non-shadow area.We use the ( ) ( ) H+1 / V+1 ratio image to enhance the low luminance and hue values of the shadow region, where H is the hue channel image and V is the value channel image.After preprocessing the satellite remote sensing image, the image of the experimental region in the RGB color space is obtained by cropping and then converted into an HSV image, as shown in Figure 7a.The ratio image is obtained, as shown in Figure 7b.Then, the ratio image is thresholded using Otsu's method (maximum inter-class variance method) to obtain the image shown in Figure 7c.These images mainly consist of the tank shadows and other small areas of speckle noise.The average area of the connected components is calculated and used as a threshold to obtain Figure 7d.Finally, a morphological closing operation is used to obtain the shadow of the tank.The structing elements (SE) we used is a 3 × 3 matrix.The values in the matrix are all 1.The result as shown in Figure 7e.If the line segment parallel to the sun's rays is obtained, the distance between two points can be obtained using the Euclidean distance as follows: p, q  −   − where (s, t) and (x, y) are the coordinates of two points and coordinates are in pixel units.S is the building shadow length, which is calculated as: and k is the pixel size in meters.Figure 9 is the interpretation diagram of shadow length calculation.The detailed steps for calculating the shadow length are as follows: First, obtain the coordinates of the upper and lower edges respectively.The boundary coordinates of the image tank shadow are extracted by the findContours() function in OpenCV.The key parameters of the function are: CV_RETR_EXTERNAL (only the outer contour is detected), and CV_CHAIN_APPROX_NONE (all the points on the contour are stored).The storage rule of the contour point is to store the point with the smallest y value first.If there is more than one point with the smallest y value, the point with the smallest value of x is taken as the first position of the storage list.Besides, all the contour points are sorted in a counterclockwise order.The right and left endpoints are the point with the largest x value and the point with the largest y value, respectively.Now we have stored the coordinates of all the contours and have identified the coordinates of the left and right endpoints.Starting from the first bit of the list, the y value of the coordinate is compared to ymax.If y is less than ymax, the coordinate is classified as the upper edges.If more than one y value is equal to ymax, the coordinate of the last ymax is the left end point.The coordinates of the remaining ymax are classified as the upper edges.Starting from the left end point, the x value of the remaining coordinates is compared with xmax.If x is less than xmax, the coordinate is classified as the lower edges.If more than one x value is equal to xmax, the coordinate of the last xmax is the right end point.The coordinates of the remaining xmax are classified as the lower edges.At this time, the remaining coordinates in the contour If the line segment parallel to the sun's rays is obtained, the distance between two points can be obtained using the Euclidean distance as follows: where (s, t) and (x, y) are the coordinates of two points and coordinates are in pixel units.S is the building shadow length, which is calculated as: and k is the pixel size in meters.Figure 9 is the interpretation diagram of shadow length calculation.The detailed steps for calculating the shadow length are as follows: First, obtain the coordinates of the upper and lower edges respectively.The boundary coordinates of the image tank shadow are extracted by the findContours() function in OpenCV.
The key parameters of the function are: CV_RETR_EXTERNAL (only the outer contour is detected), and CV_CHAIN_APPROX_NONE (all the points on the contour are stored).The storage rule of the contour point is to store the point with the smallest y value first.If there is more than one point with the smallest y value, the point with the smallest value of x is taken as the first position of the storage list.Besides, all the contour points are sorted in a counterclockwise order.The right and left endpoints are the point with the largest x value and the point with the largest y value, respectively.Now we have stored the coordinates of all the contours and have identified the coordinates of the left and right endpoints.Starting from the first bit of the list, the y value of the coordinate is compared to y max .If y is less than y max , the coordinate is classified as the upper edges.If more than one y value is equal to y max , the coordinate of the last y max is the left end point.The coordinates of the remaining y max are classified as the upper edges.Starting from the left end point, the x value of the remaining coordinates is compared with x max .If x is less than x max , the coordinate is classified as the lower edges.If more than one x value is equal to x max , the coordinate of the last x max is the right end point.The coordinates of the remaining x max are classified as the lower edges.At this time, the remaining coordinates in the contour coordinate point set are classified as the upper edges.The coordinate values of each point are the pixel's center-point coordinate values.10) and traverse the upper edge of the shadow to find angle γ.Select the coordinates of the center point of the upper edge of the shadow such that the absolute value of the difference between angle γ and supplementary angle θ of the solar azimuth is less than 0.5°.
According to sub-pixel subdivision positioning [33], subdivide the upper edge of the selected pixel.Find the point such that | − | ≤ 0.05.The distance between the selected point and the center point of the pixel at the  Convert the obtained shadow edge pixel center-point coordinate values into floating point values.Take a pixel center point on the lower edge of the shadow (e.g., point O shown in Figure 10) and traverse the upper edge of the shadow to find angle γ.Select the coordinates of the center point of the upper edge of the shadow such that the absolute value of the difference between angle γ and supplementary angle θ of the solar azimuth is less than 0.5 • .
positional relationships between the height of a building and its shadow on a remote sensing image can be divided into three cases: i) the satellite azimuth differs from the solar azimuth by more than 180°; ii) the satellite azimuth is equal to the sun azimuth; and iii) the azimuth of the satellite is not equal to the azimuth of the sun and the difference is less than 180°.Since our experimental data fall into the first case, we only cover this case in detail.
When imaging high-altitude remote sensing satellite images, if the satellite azimuth differs from the solar azimuth by more than 180°, i.e., the satellite and sun are on the opposite sides of the building, then all the shadows of the building can be observed on the satellite image, as shown in Figure 11.The geometry for this case is shown in Figure 12, where α is the solar elevation angle, h is the height of the building, and S is the length of the shadow projected by the building on the ground.According to sub-pixel subdivision positioning [33], subdivide the upper edge of the selected pixel.Find the point such that |γ − θ| ≤ 0.05.The distance between the selected point and the center point of the pixel at the lower edge of the shadow is calculated according to Equation (1).The bottom pixel of the shadow may indicate the edge pixel of the tank, so the center point of the pixel at the bottom edge of the shadow is used.
The above steps are used to calculate the length of all such segments in the tank shadow.The median, maximum, and average of these values are each used to calculate the actual length of the shadow according to Equation (2).

Oil Tank Height Calculation
The formation of building shadows in high-resolution satellite remote sensing images is related to the size of the building itself, the solar elevation angle, solar azimuth, satellite azimuth, and satellite elevation angle.In addition, the experiments in this study consider the oil tanks of a national oil reserve base where the terrain area does not have large fluctuations, does not affect the shadow of the tank on the ground, and the tank is perpendicular to the ground surface.There is also a certain interval between the tanks.In the selection of remote sensing satellite images, the shadow of the tanks completely falls on the spaced ground.Under these conditions, according to the traditional method, the height information of the tank can be obtained by using the positional relationships among the building, the satellite azimuth, the solar azimuth, and the solar elevation angle.The traditional positional relationships between the height of a building and its shadow on a remote sensing image can be divided into three cases: (i) the satellite azimuth differs from the solar azimuth by more than 180 • ; (ii) the satellite azimuth is equal to the sun azimuth; and (iii) the azimuth of the satellite is not equal to the azimuth of the sun and the difference is less than 180 • .Since our experimental data fall into the first case, we only cover this case in detail.
When imaging high-altitude remote sensing satellite images, if the satellite azimuth differs from the solar azimuth by more than 180 • , i.e., the satellite and sun are on the opposite sides of the building, then all the shadows of the building can be observed on the satellite image, as shown in Figure 11.
The geometry for this case is shown in Figure 12, where α is the solar elevation angle, h is the height of the building, and S is the length of the shadow projected by the building on the ground.In this case, when the remote sensing satellite takes an image, the geometric relationship between the shadow length of the building and the height of the building itself is given by: ℎ   . (3)

Oil Tank Top Radius Calculation
The top of the tank is assumed to be circular in the two-dimensional plane.When the tank height is known, the radius of the top of the tank is required to calculate the volume.In the proposed method, the Hough transform is used to detect the oil tank and identify its top radius.The Hough transform is suitable for detecting targets with specific shapes.It uses point-line duality to determine a specific shape in the image by a voting mechanism.The analytical geometry equation of a circle is: where ( ) c ,c is the center of the circle and r is the radius of the circle.At any point ( ) , x y on the circle in a two-dimensional image space, there will be a conical surface corresponding to it in the parameter space C1C2R.It can also be said that the points on the same circle in the two-dimensional image space are mapped to the parameter space, and a cluster of intersecting conical surfaces is In this case, when the remote sensing satellite takes an image, the geometric relationship between the shadow length of the building and the height of the building itself is given by: h = S × tan α. (3)

Oil Tank Top Radius Calculation
The top of the tank is assumed to be circular in the two-dimensional plane.When the tank height is known, the radius of the top of the tank is required to calculate the volume.In the proposed method, the Hough transform is used to detect the oil tank and identify its top radius.The Hough transform is suitable for detecting targets with specific shapes.It uses point-line duality to determine a specific shape in the image by a voting mechanism.The analytical geometry equation of a circle is: where (c 1 , c 2 ) is the center of the circle and r is the radius of the circle.At any point (x, y) on the circle in a two-dimensional image space, there will be a conical surface corresponding to it in the parameter space C 1 C 2 R. It can also be said that the points on the same circle in the two-dimensional image space are mapped to the parameter space, and a cluster of intersecting conical surfaces is formed.As shown in Figure 13, a cluster of conical surfaces intersects at point (c 1 , c 2 , r) in parameter space.The specific algorithm for detecting circles using the Hough transform is as follows: Solve for the gradient on a circular point in the image, and form a set of edge points according to the threshold processing.Initialize a voting matrix , , 0 A c c r = in parameter space C1C2R.
According to the set of edge points, take all possible values of 1 2 c ,c to obtain the corresponding value of γ.Accumulate matrix , , A c c r according to the obtained ( ) c ,c ,r , that is, ( ) ( ) , , , , 1 A c c r A c c r = + .Local extrema are detected in matrix ( ) , , A c c r to obtain the parameters ( ) c ,c ,r .

Experimental Data
The experimental data used in this paper are data from China's independent Gaofen 2 satellite.The Gaofen 2 satellite was successfully launched on August 19, 2014, and it has a spatial resolution of 0.81 m at the sub-satellite point.The imaging spectrum includes four multi-spectral spectral segments and one panchromatic segment.The detailed technical specifications and parameters of Gaofen 2 are shown in Table 1.
The remote sensing images were taken at 10:51 h on April 16, 2015, and consist of full-color image data at a resolution of 0.81 m and blue, green, red, and near-infrared multi-band image data at a resolution of 3.20 m.In addition, the solar elevation angle is 57.96°, the solar azimuth is 150.66°, and the satellite azimuth is 310.48°.The satellite image was preprocessed to form a multi-spectral fullcolor fused image, eliminating the effects of solar azimuth and satellite azimuth, and the solar azimuth and satellite azimuth differ by 159.82°, which is close to 180°.Hence, Equation (3) was used to calculate the height of the tank.We calculated the height of 24 oil tanks.The actual height of the tanks in the experimental area is 21.8 m [34].

Experimental Data
The experimental data used in this paper are data from China's independent Gaofen 2 satellite.The Gaofen 2 satellite was successfully launched on August 19, 2014, and it has a spatial resolution of 0.81 m at the sub-satellite point.The imaging spectrum includes four multi-spectral spectral segments and one panchromatic segment.The detailed technical specifications and parameters of Gaofen 2 are shown in Table 1.The remote sensing images were taken at 10:51 h on April 16, 2015, and consist of full-color image data at a resolution of 0.81 m and blue, green, red, and near-infrared multi-band image data at a resolution of 3.20 m.In addition, the solar elevation angle is 57.96 • , the solar azimuth is 150.66 • , and the satellite azimuth is 310.48 • .The satellite image was preprocessed to form a multi-spectral full-color fused image, eliminating the effects of solar azimuth and satellite azimuth, and the solar azimuth and satellite azimuth differ by 159.82 • , which is close to 180 • .Hence, Equation (3) was used to calculate the height of the tank.We calculated the height of 24 oil tanks.The actual height of the tanks in the experimental area is 21.8 m [34].

Evaluation Metrics of Extracted Shadows
The evaluation metrics proposed by Shufelt (accuracy and precision) were used to evaluate the accuracy of the detected shadow area.Recall CF and shadow detection accuracy CA are calculated as evaluation metrics, respectively [35], as follows: Here, TP is the number of shadow pixels that are correctly detected, FP is the number of unshaded pixels detected as shadows, FN is the number of shadow pixels that are detected as non-shadow pixels.

Evaluation Metrics of Oil Tank Height, Radius, and Volume
We used the absolute value and relative error to evaluate the accuracy of the estimated oil tank height, radius, and volume.The absolute error is the absolute value of the calculated height and the actual height, and the relative error is the ratio of the absolute error to the actual height.

Shadow Extraction
Images were preprocessed on the ENVI (The Environment for Visualizing Images) 5.3 software platform, which is a product of Exelis Visual Information Solutions, USA.On the Eclipse Neon 3.2 platform, which is managed by the Eclipse Foundation, Python 3.6.3and OpenCV 3.3.1 were used to implement the tank shadow extraction.To verify the effectiveness of the algorithm, the traditional histogram threshold and morphological methods were used to extract the oil tank shadow in the experimental images.These two algorithms are compared with the proposed algorithm, as shown in Figure 14.
The results are shown in Table 2.The length calculation of the oil tank shadow is carried out according to the traditional pixel number-based method, the mean method based on direct manual measurement, and the methods proposed in this paper, which include the maximum method, median method, and average method, all based on the idea of sub-pixel subdivision positioning.The specific processes of the pixel number method and direct measurement method are given in Appendix A. The shaded areas are numbered as shown in Figure 15.The results are shown in Table 2.The length calculation of the oil tank shadow is carried out according to the traditional pixel number-based method, the mean method based on direct manual measurement, and the methods proposed in this paper, which include the maximum method, median method, and average method, all based on the idea of sub-pixel subdivision positioning.The specific processes of the pixel number method and direct measurement method are given in Appendix A. The shaded areas are numbered as shown in Figure 15.The results are shown in Table 3.The results are shown in Table 3. Figure 16 compares the height of the tanks estimated by the five methods.In addition, Figure 17 compares the mean of the absolute value error, and Figure 18 compares the mean of the relative error for these methods.The root mean square error (RMES) of the pixel number method, direct measurement method, median method, maximum method, and mean method are 5.58, 0.53, 0.23, 3.63, and 0.90, respectively.18 compares the mean of the relative error for these methods.The root mean square error (RMES) of the pixel number method, direct measurement method, median method, maximum method, and mean method are 5.58, 0.53, 0.23, 3.63, and 0.90, respectively.
The results of oil tank height calculation are shown in Table 4.The results of oil tank height calculation are shown in Table 4.The proposed method uses the "HoughCircles()" function in OpenCV to identify and detect oil tanks in remote sensing images.After many experiments, when the function parameter "param1" is 132.55,"param2" is 20, "minRadius" is 48 or 49, and "maxRadius" is 51, the recognition rate of the tank sin the image is the best, as shown in Figure 19.The proposed method uses the "HoughCircles()" function in OpenCV to identify and detect oil tanks in remote sensing images.After many experiments, when the function parameter "param1" is 132.55,"param2" is 20, "minRadius" is 48 or 49, and "maxRadius" is 51, the recognition rate of the tank sin the image is the best, as shown in Figure 19.As Figure 19 shows, the numbers of tanks identified are 23 and 22 for minimum radius parameter values of 48 and 49, respectively.In this case, we selected the average of the number of tanks with the top radius correctly identified in the experiment as 22, as shown in Table 5.As Figure 19 shows, the numbers of tanks identified are 23 and 22 for minimum radius parameter values of 48 and 49, respectively.In this case, we selected the average of the number of tanks with the top radius correctly identified in the experiment as 22, as shown in Table 5.

Tank Volume
The results of tank volume estimation for the 22 identified oil tanks are shown in Table 6.In addition, Figure 20 shows the absolute error of the calculated volume and Figure 21 shows the relative error of the calculated volume.

Discussion
This paper used spectral ratio techniques in the HSV color space to detect tank shadows.This method was compared with classical algorithms, a histogram threshold algorithm, and a mathematical morphological algorithm.All three algorithms can detect the shadow of the tank.However, traditional histogram threshold algorithms and mathematical morphological algorithms erroneously detect dark-colored parts of the surface that are similar to shadows.The outline of the oil tank shadow extracted by the proposed algorithm is basically consistent with the original image area and there is less interference.An evaluation of the accuracy of the detected shadow area used the method proposed by Shufelt.Table 2 shows that the accuracy and precision of the shadow detection algorithm proposed in this article are above 96% and better than the comparison methods.The other two algorithms also have high accuracy, but their precision is low.The oil tank shadow extracted by the algorithm in this paper is better that those of other methods, which lays a good foundation for the calculation of the subsequent tank height.The accuracy of the shadow length calculation determines the accuracy of the tank height calculation.The shape of tank shadows in remote sensing images are similar to crescent shapes, in contrast to the shadows of houses, which resemble rectangles.In this paper, according to the idea of sub-pixel subdivision positioning, the lengths of all the line segments that meet the conditions in the shadow region are obtained, and the median, maximum, and mean values of the line segment length are obtained separately.We also used the classic pixel number method and direct measurement of the shadow length using ENVI software to compare the results.As can be seen from Figure 16, the heights of each tank obtained from the averaging and median methods are very close to the actual heights of the tanks.Figures 17 and 18 show that the median method has the smallest absolute and relative errors.The mean absolute error of the height of the tank calculated by the median method is within 0.25 m, and the relative error is within 1.15%, which is better than the result in [36].The smallest absolute error of calculated tank height is 0.5 m and the smallest relative error is 2.85% in [36].Additionally, The overall accuracy of our approach is compared against previous works on image-based building height estimation in Table 7.The "-" in Table 7 means no relevant data.We used the Hough transform to find the tank radius.Once the height and radius are obtained, we can calculate the volume of the tank.The absolute error of the calculated volume of the tank ranges from 416-3050 m 3 and the relative error ranges from 0.38%-2.78%.The estimation error of most tank volumes is within 3000 m 3 .Therefore, the need to estimate the actual volumes of large tanks is met by the proposed method.

Conclusions
In the method proposed in this article, we converted remote sensing images into HSV images, and the ratio image was obtained to emphasize the shadows.The ratio image was thresholded using the maximum inter-class variance of the levels, then thresholded according to area and processed using morphological operators to obtain the shadow of the oil tank.The accuracy of the proposed shadow extraction was verified by comparing it with the traditional method.The calculation of the length of the tank shadow was carried out according to the traditional pixel number method, manual direct measurement method, and proposed method (using the maximum, median, and average methods).The heights of the tanks calculated by these five methods were compared with the actual heights of the tanks.The final results showed that the median method proposed in this paper outperformed the other methods.Combined with the Hough transform to detect the radius of the top of the tank, the volume of the tank can be estimated with sufficient accuracy.
In the future, we will further study the nature of shadows in other color spaces and combine the texture features of shadows to extract tank shadows from a wider range of images.In this study, the Hough transform did not identify all oil tanks.In the future, the Hough transform and machine-learning-related algorithms should be further combined to improve the accuracy of the tank identification rate and estimation of the radius of the tank.

Figure 1 .
Figure 1.Flowchart of the proposed approach.

3. 4 . 23 3. 4 .
Oil Tank Volume Calculation3.4.1.Shadow Length CalculationMost of the shadows used in the literature to calculate building height are the shadows of rectangular buildings.In contrast, the shadows of the tanks extracted in Section 3.3 are similar to a crescent, as shown in Figure8.Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of Oil Tank Volume Calculation 3.4.1.Shadow Length CalculationMost of the shadows used in the literature to calculate building height are the shadows of rectangular buildings.In contrast, the shadows of the tanks extracted in Section 3.3 are similar to a crescent, as shown in Figure8.

Figure 8 .
Figure 8. Shadow of a tank.If the line segment parallel to the sun's rays is obtained, the distance between two points can be obtained using the Euclidean distance as follows:

Figure 8 .
Figure 8. Shadow of a tank.

Figure 9 .
Figure 9. Interpretation diagram of shadow length calculation.Figure 9. Interpretation diagram of shadow length calculation.

Figure 9 .
Figure 9. Interpretation diagram of shadow length calculation.Figure 9. Interpretation diagram of shadow length calculation.

Figure 10 .
Figure 10.Geometry used in the algorithm.

Figure 10 .
Figure 10.Geometry used in the algorithm.

Figure 11 .
Figure 11.Image in which the satellite and sun are on the opposite sides of the building.Figure 11.Image in which the satellite and sun are on the opposite sides of the building.

Figure 11 . 23 Figure 12 .
Figure 11.Image in which the satellite and sun are on the opposite sides of the building.Figure 11.Image in which the satellite and sun are on the opposite sides of the building.

Figure 12 .
Figure 12.Geometry of the case shown in Figure 10.

Figure 13 .
Figure 13.Conical surfaces in parameter space.The specific algorithm for detecting circles using the Hough transform is as follows: Solve for the gradient on a circular point in the image, and form a set of edge points according to the threshold processing.Initialize a voting matrix A(c 1 , c 2 , r) = 0 in parameter space C 1 C 2 R. According to the set of edge points, take all possible values of c 1 , c 2 to obtain the corresponding value of γ.Accumulate matrix A(c 1 , c 2 , r) according to the obtained (c 1 , c 2 , r), that is, A(c 1 , c 2 , r) = A(c 1 , c 2 , r) + 1. Local extrema are detected in matrix A(c 1 , c 2 , r) to obtain the parameters (c 1 , c 2 , r).

Figure 16
Figure 16  compares the height of the tanks estimated by the five methods.In addition, Figure17compares the mean of the absolute value error, and Figure18compares the mean of the relative error for these methods.The root mean square error (RMES) of the pixel number method, direct measurement method, median method, maximum method, and mean method are 5.58, 0.53, 0.23, 3.63, and 0.90, respectively.The results of oil tank height calculation are shown in Table4.

Figure 16 .
Figure 16.Comparison of tank heights obtained by various methods.

Figure 16 .
Figure 16.Comparison of tank heights obtained by various methods.

Figure 17 .
Figure 17.Mean of the absolute error of the tank height calculated by different methods.

Figure 17 . 23 Figure 18 .
Figure 17.Mean of the absolute error of the tank height calculated by different methods.

Figure 18 .
Figure 18.Mean of the relative error of the tank height calculated by different methods.

Figure 20 .
Figure 20.Absolute error of the estimated volumes of the tanks.

Figure 20 .
Figure 20.Absolute error of the estimated volumes of the tanks.

Figure 21 .
Figure 21.Relative errors of the estimated volumes of the oil tanks.

Figure 21 .
Figure 21.Relative errors of the estimated volumes of the oil tanks.

23 Figure A2 .
Figure A2.Direct manual measurement of shadow length.

Figure A2 .
Figure A2.Direct manual measurement of shadow length.
Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 23 coordinate point set are classified as the upper edges.The coordinate values of each point are the pixel's center-point coordinate values.Convert the obtained shadow edge pixel center-point coordinate values into floating point values.Take a pixel center point on the lower edge of the shadow (e.g., point O shown in Figure

Table 2 .
Comparison of the shadow detection results.

Table 2 .
Comparison of the shadow detection results.

Table 4 .
Tank heights obtained by various methods.

Table 6 .
Calculated volumes of the oil tanks.

Table 7 .
Algorithm performance against previous works.