基于大气所铁塔照片的能见度识别

利用照片识别能见度是一个非常有意思的话题，今天我就来分享一下最近试验成功的基于大气所铁塔照片识别能见度的方案。

大气所气象铁塔

从大气所官网（https://iap.cas.cn）页面的下方，可以进入铁塔俯瞰图的主页：http://view.iap.ac.cn:8080/imageview

下面是主页里对铁塔的简介：

中国科学院大气物理研究所325米气象观测塔（简称气象塔）建成于1979年8月，主要服务于首都北京的空气污染研究和大气边界层、大气湍流扩散等研究。在世界众多气象塔中，大气所气象塔具有独特的优势，其高度是世界第三，塔上的观测平台分布较细致。气象塔位于39°58′N，116°22′E，海拔高度49米。在气象塔280米高度南北两个方向设置的两部高清相机24小时运行，每半小时自动拍照，实时监测周边天气状况和地表状况。您现在看到的是半小时内更新的奥运区实时图，左侧为东北向，右侧为西南向照片。

铁塔相机的拍摄视角是静止的，大型遮挡物较少，可以远眺很远的物体。因此它可以为我们做能见度识别提供位置稳定且景观丰富的照片。鉴于西南方向镜头容易受强光影响出现过度曝光的情况，因此我们选用东北方向的镜头做能见度识别的工作。

景观测距与分离

如果要通过照片来计算能见度，一种很直观的方案就是把照片中的各种不同距离的景观标注出来，然后通过判断最远的可辨识的景观来确定能见度。因此我们的第一个工作就是对景观进行测距，得益于 Google Earth 等工具的帮助，让我们的测距工作变得容易很多。

我们先根据公开的坐标信息在 Google Earth 上标注好铁塔的位置，然后通过 Google Earth 上的卫星图片与照片中的建筑物做人工对照和地图标注。

在完成地图标注以后，我们可以使用 Google Earth 的测距工具，测量铁塔与各个景观之间的距离。

完成定距之后，我们还需要借助一些绘图工具，将不同景观的选框坐标记录下来，最后保存为一个 json 配置文件：

{
    "华庭嘉园": {
        "x1": 896,
        "y1": 2285,
        "x2": 3154,
        "y2": 3246,
        "distance": 1.27
    },
    "水立方": {
        "x1": 2420,
        "y1": 1886,
        "x2": 3524,
        "y2": 2191,
        "distance": 2.18
    },
    "盘古大观": {
        "x1": 1691,
        "y1": 1079,
        "x2": 2385,
        "y2": 2104,
        "distance": 1.71
    },
    ... // 篇幅限制省略部分内容
}

完成配置之后，我们就可以根据的选框信息，对照片中的各个景观进行切割分离，参考代码：

import cv2

def crop_image(x1, y1, x2, y2, image_path, output_path):
    img = cv2.imread(image_path)
    cropped_img = img[y1:y2, x1:x2]
    cv2.imwrite(output_path, cropped_img)

景观辨识度分析

根据前述的思路，我们将不同景观做分离的目的是判断不同距离景观的辨识度，从而计算最终能见度。那么现在的关键问题是如何来判断景观的辨识度。众所周知，清晰照片的色彩丰富度较高，色彩分布较分散；而模糊照片的色彩丰富度较低，色彩分布较集中。鉴于这种思路，加之一些实验的佐证，我们提出下面一种辨识度评分方案。

其中 是指 Recognizability（辨识度）； 是指 Standard Deviation（标准差）； 是指 Spectral Width（谱宽）。在这里， 标准差是指将图片以灰度图方式读取的色彩值（取值范围在0-255之间）的标准差， 谱宽是指将图片以灰度图方式读取后的色彩值中，最大值与最小值之间的差值。

该公式可以兼顾色彩丰富度和离散度， 值越大则辨识度越高，越低则辨识度越低。下图显示的是不同能见度条件下奥林匹克塔所计算的 值及其对应的景观照片。

可以看到， 值在 200 以下基本可以被视为难以辨识，我们可以根据自己对辨识度的容忍度来设置一个阈值，来判定景观结果是“可见”还是“不可见”。在这里我将阈值设置为 250。

参考代码实现：

import matplotlib.pyplot as plt
from skimage.color import rgb2gray

def is_clear(image_fp, threshold=250):
    # 读取图像并转换为灰度图像
    image = rgb2gray(plt.imread(image_fp))
    # 将图像的数据类型转换为uint8，并将其范围缩放到0-255
    image = (image * 255).astype("uint8")
    std = image.std()
    span = image.max() - image.min()

    score = std * span
    if score > threshold:
        return True
    else:
        return False

计算全图能见度

当我们可以判断单个景观是“可见”还是“不可见”以后，判断全图能见度就变得简单了，只需要把那些被判定为“可见”的景观的保留下来，然后获取其中最远的景观距离即可。参考代码：

import os
import shutil

CONFIG_FP = os.path.join(os.path.dirname(__file__), "config.json")

def get_config(config_fp=CONFIG_FP):
    with open(config_fp) as f:
        config = json.load(f)

    return config

def analysis_visibility(image_fp, threshold=250, tmpdir="./tmp", keep_tmp=False):
    config = get_config()
    visibles = []
    for key, value in config.items():
        x1, y1, x2, y2 = value["x1"], value["y1"], value["x2"], value["y2"]
        distance = value["distance"]
        savefp = os.path.join(tmpdir, f"{key}.jpg")
        os.makedirs(tmpdir, exist_ok=True)
        crop_image(x1, y1, x2, y2, image_fp, savefp)
        if is_clear(savefp, threshold):
            visibles.append(distance)

    if not keep_tmp:
        shutil.rmtree(tmpdir)

    return round(max(visibles), 0)

至此，我们就可以使用 analysis_visibility 傻瓜式地获取铁塔照片的能见度值了。下面我们用几张实战照片来看一下识别的效果：

本项目代码已在 github 上开源共享，仓库地址：https://github.com/caiyunapp/tower-eye/tree/main