UAV-BASED THERMAL ANOMALY DETECTION FOR DISTRIBUTED HEATING NETWORKS

District Heating Systems (DHS) distribute heat in terms of hot water or steam. Loss of media (water or steam), and thus energy, is expensive and has a negative impact on the environment. It is therefore of great interest to develop techniques to detect and localize potential leakages fast and cost-effectively. To avoid interference with the operating process of a DHS, airborne thermography comes into place. The use of Unmanned Aerial Vehicles (UAV) as a flexible and low-cost platform equipped with a Thermal InfraRed (TIR) camera is a promising alternative to a conventional manned flight. This paper describes a method for automatic thermal anomaly detection of a DHS. Thermal data acquisition using a UAV is followed by photogrammetric processing of the TIR images. In this way, a thermal orthophoto is produced. The next step is an identification of anomalies by means of image analysis. We apply the Laplacian of Gaussian (LoG) blob detector to find high temperature regions in areas of interest of a thermal orthomosaic. This area of interest is defined around the DHS position in the images as defined in a Geographic Information System. Finally, segmentation and classification are employed to reduce false alarms and localize thermal anomalies. An experimental evaluation using real-world data is presented, showing that the developed method deliverers promising results.


INTRODUCTION
Conventionally operated District Heating Systems (DHS) use the principle of cogeneration and distribute heat which is generated as waste heat during the generation of electricity. For the integration of renewable heat, district heating networks are often characterized by a high degree of flexibility, especially when low-temperature (or lowenergy) heat sources are to be made use of. While these networks have a high technological relevance, they come with high costs for construction, maintenance and repair (Sernhed and Jönsson, 2017).
Today's maintenance strategy is based on a statistical assessment of the state of damage of the networks. Relevant pipeline sections are gradually taken out of operation, emptied and inspected for possible errors, which is slow and costly. This problem can be remedied by airborne Thermal InfraRed (TIR) imaging of heating networks, which avoids interference with the operating process. This technology makes it possible to detect temperature differences of a surface. Unusually high radiation and temperature differences (anomalies) give indications for emerging water and thus damage to the district heating system (Ljungberg et al., 1988;Axelsson, 1988). Currently, airborne thermal flights usually are carried out over a large area with airplanes or helicopters. However, their application is limited due to cost considerations and high planning effort, in particular for small-scale networks.
Particularly in the last years, a rapid development in the Unmanned Aerial Vehicles (UAV) technology equipped with different sensors could be observed (Nex et al. 2014;Colomina and Molina et al. 2014). In our context, the use of UAV as a flexible and low-cost platform equipped with a TIR camera is a promising alternative to airplanes and helicopters. Using a UAV, it is possible to acquire data of a region within a comparatively short period of time at a reasonable cost, in order to detect possible leakages.
In light of the foregoing, this paper aims to develop and test an economical, automated, measurement system for the management of DHS based on thermal imaging. The focus is to define UAV operation steps for data acquisition and then to identify thermal anomalies related to a DHS by means of automatic image analysis methods.

RELATED WORK
Monitoring techniques for DHS evolved over the years. One of the common methods is to measure the water loss between inlet and outlet. Despite being effective, those methods are only capable of providing an approximate location for a leakage. Methods for large-scale monitoring by manned flight thermography provide an alternative, these were studied e.g. by Ljungberg et al. (1988) and Axelsson (1988). Friman et al. (2014) show a method for automatic analysis of TIR images to localize leakages in pipes of a DHS. The TIR images were captured from manned aircraft at 800 meters altitude, which leads to 24 cm ground sampling distance (GSD). Leakages are located by selecting the upper percentile from the distribution of the heat radiation for a pipe, while pipe localization is done by projecting Geographic Information System (GIS) data of the DHS into the acquired images. In order to reduce the false alarm rate related to abnormal temperatures around buildings, for example from warm chimneys, vents and entrances, buildings are segmented by a watershed transformation applied to the orthophoto followed by a classification of the resulting segments. It is shown that the Adaboost classifier achieves an accuracy of 97% on the training data and 86% on the test data. Berg et al. (2014) also address the problem of reducing the number of false alarms among automatically detected potential leakages in district heating networks. The leakages are detected in images captured by an airborne TIR camera, a region-based shape representation of the detection and proximity to buildings. Feature selection was done by calculating the Mahalanobis distance between class means for each feature. The final feature set consists of eight features. The classification problem was formulated as binary classification, where the target was to determine if the detection is a leakage or not. While different types of linear and nonlinear classifiers were employed in the task, the Random Forest classifier achieved the best results. Zhong et al. (2019) used a saliency analysis method for DHS pipe leakage detection using remotely sensed infrared imagery and optical imagery captured by manned flights and UAVs with a min. GSD of about 20 cm, and GIS data. A saliency map is an image that emphasises the object with respect to its surroundings. Itti et al. (1998a) established one of the first saliency computation models. In this sense pipeline leakages can be seen as salient in infrared imagery as the leaked material creates a local high-temperature area. Zhong et al. (2019) develop a DHS leakage detection framework characterized by local and global saliency analysis (LGSA). The local saliency map shows image regions distinct with respect to their local neighbourhoods, where intensity and orientation features are adopted in the analysis. On the other hand, the global saliency map represents saliency of an image region by calculation of the seldomness of features over the entire scene. Global saliency in these contents enhance homogenous regions inside DHS leakages detected by local saliency. In the last step, the infrared saliency map is created by fusing the local and the global saliency map using the dot product. The infrared saliency map is an index map, where a pixel with a high index has a high probability to be part of a high-temperature anomaly area. The authors make the assumption that they have at hand GIS data for the location of the DHS. In the absence of such data, they use the road network as a substitute for their location, as according to the authors these pipes can typically be found alongside roads. At this point optical images come into place to extract road areas, which are subsequently used as complementary information to the DHS GIS data. Finally, instead of using thresholding techniques to detect potential leakages, maximum entropy segmentation is done on the final fused saliency map. The validation of the proposed method was carried out on data captured in mid and high altitude in Gävle in Sweden and Datong in China. The field experiments revealed that the proposed LGSA method successfully detected the leakages of the underground DHS pipelines with a detection accuracy of 90%, and it outperformed the existing methods in false alarm reduction. However, in the absence of DHS GIS data, the success rests on the mentioned assumption that pipes can be found alongside roads.
In the work mentioned above, all flights, except those reported in Zhong et al. (2019), were carried out using an manned aircraft at a flying height of hundred meters and more, employing a cooled IR camera,. Of course, flying at higher altitude provides larger area coverage, but the result is a lower GSD. In addition, employing a manned aircraft needs complex preparations and is rather costly. This paper deals with temperature anomaly detection in thermal orthomosaics. The main contribution is the usage of a blob detector, followed by segmentation, to detect abnormal temperature variations related to a DHS from UAV TIR orthophotos, generated using rigorous photogrammetric techniques. As the TIR camera is roughly calibrated we are able to convert the grey values of the thermal image into temperatures. To reduce false alarms, we use GIS data, incl. pipe locations and height data, and we do not rely on the assumptions introduced in Zhong et al. (2019). In our investigation, all flights were done at an altitude below 50 meters using a UAV equipped with both, a TIR and an optical camera, giving us a much higher GSD (5.2 cm TIR, 1 cm optical) and the flexibility of easily choosing when and where to do the flights

System overview
Before describing the developed anomaly detection method, we first provide an overview of the automatic leakage detection system we are in the process of developing, as shown in Figure 1. The system consists of three steps: photogrammetric processing of TIR and optical images, anomaly detection followed by classification. Photogrammetric processing was done employing rigorous photogrammetric techniques (see Sledz et al., 2018 for details). In the following, we concentrate on anomaly detection, while a discussion with respect to anomaly classification is beyond the scope of this paper. Nevertheless, we provide a short discussion in the last chapter.

Hot spot detection
Based on an orthomosaic generated by rigorous photogrammetric processing, the aim of the anomaly detection method is to determine areas with abnormal temperatures on the ground that are produced by leakages of a DHS. This is done by detecting hot spots in the thermal orthomosaic. Hot spots are brighter than a surrounding in the thermal images, so they can be regarded as blobs. A hot spot is defined as an anomaly, if the temperature difference to the surroundings is higher than a predefined threshold.
As a pre-processing step, the search space in the orthomosaic is limited to only a buffer around the pipe locations based on GIS information usually available from the DHS operator. The buffer size needs to be selected with respect to the geometric accuracy of the GIS data and the process of orthoprojection. The GIS is rasterised to the same ground resolution as the orthomosaic. The outcome is a binary raster, where pixels with the values "1" represent the position around the pipe. The next step is to detect hotspots in the orthomosaics. For this task, we use the Laplacian of Gaussian (LoG) blob detector in scale space (Lindeberg 1998). Each blob is represented by its centre coordinates ( , ) and the corresponding scale (σ) (see upper left images in Figure 3 and Figure 4) that is matched with the scale of the Laplacian. As we have observed that an elliptical blob representation leads to a better anomaly localisation, we subsequently determine the major and minor axis of an ellipse fitted to the blob area. The axes are calculated from the eigenvalues and eigenvectors of the area, derived from the weighted Hessian matrix (H) of the smoothed image ( ) as described at equation (1). In case one of the axis of the ellipse is too small (we use a couple of pixels as threshold), the scale of the blob (σ) is used as the length of the axis instead of the axis calculated from H.
where ( , ) is the Gaussian filter and is an image smoothed by the Gaussian filter with blob scale σ. and are the size of the Gaussian filter (we use an interval of [−3 , 3 ].)

Clustering and anomaly detection
In the next step, anomalies are detected based on the found hot spots by subdividing the blob areas into segments with higher and lower temperatures; those with higher temperatures are then considered to be anomalies. After determining minor and major axis of the ellipse, the ellipse is enlarged by a factor of 1.5. Such the factor was chosen to provide enough pixels from the surroundings for the following clustering. Temperature values of pixels in the enlarged ellipse are clustered by k-means with k = 3. The upper right images in Figure 3 and Figure 4 show two green ellipses: the inner one is the calculated ellipse, while the outer one is the area in which temperature values are clustered. The idea of using three clusters is to ensure a better separation of higher temperatures. The outcome is three segments that are characterized by their mean temperature [ , , ℎ ℎ ]. In the bottom left images in Figure 3 and Figure 4 the results of clustering are presented, where blue, yellow and red colours represent the clusters with low, mid and high temperature. One drawback of the usage of three clusters instead of two is that pixels that indeed belong to a hot spot may be clustered into the middle cluster. To overcome this problem, fusion between and ℎ ℎ is done in the following manner: check if is close (less than a threshold ℎ ℎ ) to ℎ ℎ and far enough (larger than another threshold ℎ ) from . In case, these conditions are satisfied, the clusters of and ℎ ℎ are merged (see bottom right image in Figure 3). Figure 4 shows the examples where such a condition is not met; as the result the hot spot area continues to have three different segments. The final decision -is the investigated hot spot an anomaly or not? -is taken by checking, if the temperature difference between the region with the highest temperature ( ℎ ℎ ) and its surroundings ( or ) is higher than a predefined threshold( ℎ ). All steps to decide if a hotspot is an anomaly are shown in the decision tree in Figure 2. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition)

EXPERIMENTS AND RESULTS
In our experiments, we used a DJI M200 UAV, a vertical take-off and landing (VTOL) quadcopter equipped with GNSS receiver, Inertial Measurement Unit (IMU) and barometer. The camera in use was the DJI Zenmuse XT2. Its gimbal-stabilized, dual-sensor design rigidly combines a FLIR radiometric thermal imager and a CMOS optical camera to capture thermal and visual image data while in flight. The uncooled FLIR thermal sensor has a resolution of 640x512 pixels, while the optical camera comes with a 12MP sensor. The thermal camera in use has a frame rate of 9 HZ. The DJI Zenmuse XT2 comes with two MicroSD cards for optical and thermal data storage.
A total number of 11 flights were carried out in different regions in Hannover, Germany. Flight planning was based on the TIR camera, due to lower GSD. Based on a pixel size of 17 µm and a focal length of 13 mm, the flying height was set to 40 meters above ground, leading to a ground resolution of 5.2 cm for the TIR images (the corresponding GSD of the optical camera was 0.95 cm). This flying height was chosen as a compromise between obtaining a small GSD and a flying height above the objects in the scene. Another factor is the legal restriction not to exceed 50 meters flying height, because of the UAV operational requirements in Germany. 80% overlap in both directions for the TIR images was chosen for flight strips to ensure enough tie points in the TIR images, which exhibit rather low contrast. With respect to environmental requirements, all flights were done at ambient temperatures below 5 ˚C, with the first morning sun light and wind speed lower than 12 m/s.
Rigorous photogrammetric processing of the images was carried out using Agisoft Metashape including GNSS measurements for the position of the image projection centres. Table 1 shows the average photogrammetric statistics of all the flights. Due to finer GSD and the better image quality, the optical block produces more tie points (sparse point cloud) and also more points for the DSM (dense point cloud). In addition, the optical camera has a wider field of view and thus covers a slightly larger area.
Tie point matching was done within one group of images (thermal or optical) only. However, observations for camera position were treated differently for the optical and the thermal block. For the optical block, GNSS positions were used as observation for the projection centres and a bundle adjustment was computed. Due to the fact that both cameras are rigidly fixed in the gimbal, the complete estimated optical camera exterior orientation was then subsequently used as observations in a second bundle adjustment to compute the exterior orientation of the thermal camera, followed by dense matching to generate a digital surface model (DSM).
The final product is a thermographic orthomosaic of the region of interest. Because of the low accuracy of the lowcost GNSS in use, a misalignment between the image block and the GIS data exists. The misalignment was corrected manually. Examples of the produced orthomosaics for three test sites (with detected anomalies) are shown in Figure 10, Figure 11 and Figure 12. Area coverage (m 2 ) 12092 9335 During our experiments we concluded, that ℎ ℎ = 1˚C and ℎ = 2˚C provide good cluster fusion results. Due to many reasons, among them variations in water temperature of the DHS, pipe depth installation, pipe type (width, insulation type, etc.) and the material between the pipe and the surface, not all pipes produce enough heat signature to be observed by the TIR sensor. Based on literature research, our initial assumption was, that a leakage should produce a temperature difference of 5˚C and more (FFI personal communication) with respect to its surroundings. Due to the fact, that no actual leakage was present in the captured data, ℎ was set to 2˚C. Around 200 such anomalies were then detected in the images. Next, labels were manually assigned to the detected anomalies:  Table 2 summarises the number of detected anomalies for each label. From these numbers it becomes clear that the amount of anomalies interesting for leakage detection is indeed very small.
Not surprisingly, with respect to anomalies related to the DHS, we thus detected a large number of false alarms. Nevertheless, many of these false alarms are easy to detect and eliminate. In particular, anomalies related to 3D objects such as buildings, street lamps, cars and people can be identified using the GIS building layer and by inspecting a normalised DSM derived from the DSM used for orthoprojection, e.g. by one of the methods studied by Sithole and Vosselman (2004). In addition, moving objects (cars, people) can also be found by comparing flights carried out on different days.
In the remaining results, two type of anomalies were observed. The first type refers to anomalies caused by a hot area; those are relevant in our context. The second type of anomaly is in close proximity to very cold objects (see Figure 9 for an example). In this case, the temperature difference of the area and its surroundings did not exceed ℎ , however the proximity to a cold object was the cause for trigger for the anomaly. The same reason mostly applies to anomalies related to buildings. As stated by Friman et al. (2014) and as observed in our findings, building roofs indeed appear to be much colder, especially in early hours. Here, absolute temperatures need to be investigated to distinguish these cases from the sought anomalies related to the DHS.
In order to further identify false alarms and thus reduce the results to those meaningful for leakage detection, a classification of the remaining anomalies is the next step we will carry out.

CONCLUSION AND DISCUSSION
This paper showed an automatic anomaly localisation procedure for a district heating system (DHS) based on blob detection with the utilisation of GIS data. Such methods are a valuable tool for DHS monitoring. In principle, UAV-based thermography thus provides a fast response for thermal anomaly localisation with relatively low cost compared to manned flight thermography. The method described in this paper is one of the important steps to achieve a fully automated approach for DHS monitoring.
In future work, we intend to further integrate the thermal and the optical images in various ways. From the photogrammetric point of view, an additional step could be defined: thermal and optical block registration. As shown in Table 1, the optical block achieves a better accuracy than the thermal one. A closer integration of the two sets of images will thus lead to a higher accuracy of the thermal orthomosaic. Thermal and optical block registration could be based on optical and thermal Ground Control Points (GCP). In such way, both blocks, thermal and optical are georeferenced to a same datum. However, such a step will require more complex flight planning in terms of GCP placement and measurement. On the other hand, registration between the thermal and optical point cloud provides a fully automatic way. As was proposed by Maset et al. (2017), the Iterative Closest Point (ICP) algorithm is a suitable solution. Finally, a common bundle adjustment between thermal and optical images, which requires multi-modal image matching (e.g. Mehltretter and Heipke, 2018), is worth an investigation. The outcome is a more consistent dataset a more detailed DSM and a geometrically more accurate thermal orthomosaic.
The final step is an automatic classification of the detected anomalies in order to reduce the number of the false alarms. Recently, convolutional neural networks (CNN) have emerged as a powerful tool for many classification tasks. Despite their huge advantage, the major drawback is the requirement for large amounts of training data. Alternatives exist in more conventional approaches such as support vector machines and random forest classifiers, see also Berg et al. (2014). One of the advantages of a closer integration between optical and thermal images is the fact that both image types can be simultaneously used for the classification. In addition, a refined DSM will also help in obtaining accurate results. Finally, with respect to the potentially small and unbalanced anomaly dataset (see Table 2), all images of the flight can be taken advantage of in searching for anomalies (as opposed to only the orthomosaic) by back-projecting an anomaly found in the orthomosaic to all related images using the accurate DSM and the image orientation parameters. .