唔，最近在做作業的時候，一些實驗內容涉及到了用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);

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)