图像辐射校正之基于BP的辐射校正代码实现

在无人机低空遥感的应用过程中，如何从无人机获取的遥感影像的DN值（Digital Number）反演出地表地物真实反射率是一个十分重要的问题。在理想情况下，太阳辐射出的能量在到达地球表面后，有一部分被地物所吸收，还有一部分经过地物反射后被无人机所接收，而这一部分与太阳辐射出来的总能量的比值就是地物的地表反射率。然而在实际的数据采集过程中，即使无人机的飞行高度比较低，由于存在环境因素、传输介质等影响因素，地物的辐射反射的能量也还是会受到大气微粒的吸收和散射作用，导致目标地物的辐射信号减弱，另一方面会受到背景环境的干扰，不能够反映地物目标真实的光谱反射率、光谱辐射亮度等物理特征。因此无人机传感器接受到的地物辐射信号在未经大气校正之前并不能准确的反映地物目标的真实辐射特征。虽然大气校正的研究已经取得了不错的进展，但是目前大多数的基于辐射传输模型的大气校正方法都是针对卫星传感器的，且用于无人机低空遥感影像数据的大气辐射校正方法都有着模型参数复杂、输入条件严苛、计算过程耗时较长的问题，不能较好的利用无人机传感器价格低廉、获取手段简单的特性。

因此采取BP神经网络的辐射校正方案能够较大提高辐射校正速度和效率。

辐射校正主要代码如下：（BP网络基于tensorflow框架实现）

d = 5  # 输入节点个数
    l = 1  # 输出节点个数
    q = 10  # 隐层个数,采用经验公式2d+1
    train_num = 480  # 训练数据个数
    test_num = 240  # 测试数据个数
    eta = 0.13  # 学习率
    error = 0.002  # 精度
    w1 = tf.Variable(tf.random_normal([d, q], stddev=1, seed=1))  # seed设定随机种子，保证每次初始化相同数据
    b1 = tf.Variable(tf.constant(0.0, shape=[q]))
    w2 = tf.Variable(tf.random_normal([q, l], stddev=1, seed=1))
    b2 = tf.Variable(tf.constant(0.0, shape=[l]))

    # 输入占位
    x = tf.placeholder(tf.float32, shape=(None, d))
    y_ = tf.placeholder(tf.float32, shape=(None, l))

    # 构建图：前向传播
    a = tf.nn.sigmoid(tf.matmul(x, w1) + b1)  # sigmoid激活函数
    y = tf.nn.sigmoid(tf.matmul(a, w2) + b2)
    mse = tf.reduce_mean(tf.square(y_ - y))  # 损失函数采用均方误差
    train_step = tf.train.AdamOptimizer(eta).minimize(mse)  # Adam算法
    a = excel()  # 返回整个函数的值
    excel_data = []
    excel_all_y_trues = []
    bandMax = 850
    elevationMax = 500
    DNMax = 65535
    solarAltitudeAngleMax = 90
    temperatureMax = 40
    for i in range(len(a)):
        # 进行归一化处理
        excel_data.append([a[i][0] / bandMax, a[i][1] / elevationMax, a[i][2] / DNMax, a[i][3] / solarAltitudeAngleMax,
                           a[i][4] / temperatureMax])
        excel_all_y_trues.append(a[i][10])
    # Define dataset
    data = np.array(excel_data)
    [m,n]=data.shape
    all_y_trues = np.array(excel_all_y_trues)
    all_y_trues=all_y_trues.reshape([m,1])





    #trainornot=0

    #if trainornot>0:
    #    train_y()
    #    message_askokcancel("提示", "BP网络模型训练完成！")
    band=450/850

    elevation = 100/500


    solarAltitudeAngle = 62/90

    temperature = 25/40

    img = cv2.imdecode(np.fromfile("./P45001420.tif",dtype=np.uint16),2)
    img=np.array(list(range(1,65535)))
    # img=reduce(operator.add, img)
    img = img.reshape((-1, 1))
    DNMax=max(img)
    img = img / 65535
    l = img.shape
    #DNN = 20000 / 65535
    #img = DNN * np.ones(l)*/
    band=np.ones(l)*band
    elevation = np.ones(l) * elevation
    solarAltitudeAngle = np.ones(l) * solarAltitudeAngle
    temperature =np.ones(l) * temperature
    tt=np.ones(l)*0.63
    testX=np.append(band,elevation,axis=1)
    testX = np.append(testX, img, axis=1)
    testX = np.append(testX, solarAltitudeAngle, axis=1)
    testX = np.append(testX, temperature, axis=1)

    testX=testX.astype(np.float32)
    with tf.Session() as sess:
        init_op = tf.global_variables_initializer()  # 初始化节点
        sess.run(init_op)

        STEPS = 0
        while True:
            sess.run(train_step, feed_dict={x: data, y_: all_y_trues})
            STEPS += 1
            train_mse = sess.run(mse, feed_dict={x: data, y_: all_y_trues})
            if STEPS % 10 == 0:  # 每训练100次，输出损失函数
                print("第 %d 次训练后,训练集损失函数为：%g" % (STEPS, train_mse))
            if train_mse < error:
                break
        print("总训练次数：", STEPS)
        a = tf.nn.sigmoid(tf.matmul(testX, w1) + b1)
        y = tf.nn.sigmoid(tf.matmul(a, w2) + b2)
        zz=sess.run(y)