唔，最近在做作业的时候，一些实验内容涉及到了用ENVI处理遥感图像，然后自己手动操作软件一遍遍的输入各种参数神马的感觉挺无聊。。。。然后决定自己用python里面的opencv库写个脚本批处理图像反射率的计算试试~

核心步骤就是 遥感影像光谱辐射定标 →大气校正→计算反射率这三步了

①、遥感影像的光谱辐射定标

由遥感器的灵敏度特征引起的辐射畸变主要由其光学系统或光电转换系统的特征形成的,光电转换系统的灵敏性特征通常很重复,其校正一般是通过定期的地面测定值进行的。

遥感器光谱辐射定标时采用以下转换算式:

遥感器各波段偏移与增益值从论文找了找后，找到这么一张表~

那么这么个函数就能定标咯：

def computL(gain,Dn,bias):

return (gain*Dn+bias)

②、遥感影像的大气校正

任何一种依赖大气物理模型的大气校正方法都需要先进行遥感器的辐射校准。

公式是这个咯(Chavez P S，Jr. Image -Based Atmospheric Correction Revisited and Improved Photogrammetric Engineering and Remote Sensing, 1996,62,1025 -1036)

其中:Lhazel——大气层光谱辐射值;LI，min——遥感器每一波段最小光谱辐射值;LI,1%——反射率为1%的黑体辐射值。

关于LI，min和LI,1%的计算公式就省略了啊，感兴趣的同学可以自己去查查论文~

而计算Lhazel需要的参数可以从遥感图像的头文件中获得一部分，还有一部分是固定的参数~这些都藏在ENVI的背后，不过自己写脚本的时候找出他们还是废了一番功夫的。

计算Lhazel的代码如下：

#ESUN

ESUNI71=196.9

cos=math.cos(math.radians(90-41.3509605))

#

Lmini=-6.2

Lmax=293.7

#

Qcal=1

Qmax=255

LIMIN=Lmini+(Qcal*(Lmax-Lmini)/Qmax)

LI=(0.01*ESUNI71*cos*cos)/(math.pi*D*D)

Lhazel=LIMIN-LI

③、计算遥感影像的反射率

根据太阳辐射和大气传输原理与过程,TM/ETM+数据地面反射率反演的数学模型可综合表达为:

其中:ρ——地面相对反射率;D——日地天文单位距离;LsatI——传感器光谱辐射值,即大气顶层的辐射能量;LhazeI——大气层辐射值;ESUNl——大气顶层的太阳平均光谱辐射,即大气顶层太阳辐照度;SZ——太阳天顶角。

这里提一下其中两个参数的计算公式：

日地天文单位距离 D=1 -0.01674 cos(0.9856×(JD-4)×π/180);

(JD为遥感成像的儒略日(Julian Day)，计算公式为：

JD=K-32075+1461*(I+4800+(J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4

I、J、K分别为年、月、日

有了这些，最后就能直接算出来反射率啦，粗糙代码如下，因为是写着玩的，也没怎么处理：

不过需要注意的是，遥感图像进行计算跟输出的时候，需要使用uint16类型的数组来存储的(uint8长度不够啊。。)

一些参数涉及到浮点数计算，如果对处理结果有极高要求的话，最好使用专门的科学运算库(像我这种渣学校才不介意这些)

import cv2

import numpy as np

import math

img1=cv2.imread('F:\L71121040_04020030220_B10.TIF')

#图像格式转换

img10=cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)

#计算JD

I=2003

J=2

K=20

JD=K-32075+1461*(I+4800+ (J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4

#设置ESUNI值

ESUNI71=196.9

#计算日地距离D

D=1-0.01674*math.cos((0.9856*(JD-4)*math.pi/180))

#计算太阳天顶角

cos=math.cos(math.radians(90-41.3509605))

inter=(math.pi*D*D)/(ESUNI71*cos*cos)

#大气校正参数设置

Lmini=-6.2

Lmax=293.7

Qcal=1

Qmax=255

LIMIN=Lmini+(Qcal*(Lmax-Lmini)/Qmax)

LI=(0.01*ESUNI71*cos*cos)/(math.pi*D*D)

Lhazel=LIMIN-LI

def copy(img,new1):

new1= np.zeros(img.shape,dtype='uint16')

new1[:,:] = img[:,:]

def computL(gain,Dn,bias):

return (gain*Dn+bias)

if __name__ == '__main__':

print 'D=',D

print 'cosZS=',cos

print 'Lhazel=',Lhazel

#计算图像反射率

result=np.zeros(img.shape,dtype='uint16')

for i in range(0,img.shape(1)):

for j in range(0,img.shape(0)):

Lsat=computL(1.18070871,img10[i,j],-7.38070852)

result[i,j]=inter*(Lsat-Lhazel)*1000

#保存图像

cv2.imwrite("F:\\result.tif", result)

cv2.namedWindow("Image")

cv2.imshow("Image", result)

cv2.waitKey(0)