許久沒有更新遙感類軟件開發了,都有點生疏了,這一次我帶來了一個老的算法,新的東西, 為什么這么說呢,我們知道Landat8、Landsat5等影像,單個影像去做溫度反演,並沒有什么太大的難度,
但是呢,如果遇到大批量、多源的數據怎么辦呢,如果一景景去調參、去設計模型,那就是在太浪費時間了,我看過市面上面所有的同類溫度反演軟件,幾乎沒有人做到全自動Landat系列的溫度反演,
但是規划類的業務應用需求,往往需要大尺度范圍的高精度溫度反演結果,這就有點尷尬了,源於這個需求,我開發了市面上第一個全自動批量化Landsat8溫度反演軟件,在此拋磚引玉,相信山外有山,
定有其他高手能做出更優秀的全自動溫度反演軟件,如需交流,qq:1044625113,下面我講一下開發流程與算法模型基礎:
1.算法模型基礎
算法層面,我不講太多,大家可以參考一些覃志豪老師、國外的一些Professor寫的經典論文,歸根到底就是兩類:單通道和多通道。其中針對Landsat8數據源,覃志豪老師已經明確了,單通道更為合適,這是由於該
數據源有一個熱紅外通道出現了問題(不知道這樣表述是否明確)。我這里全自動的溫度反演算法,就是采用經典的覃志豪單通道算法,為了提高自動化程度,一些大氣參數,我采用了NASA開發的網站進行查詢獲取,比如大氣
向上輻射參數、大氣透過率等。這里我采用了一個技巧,大氣參數需要手動獲取,我采用爬蟲技術自動獲取大氣參數,這里就“曲線救國”的方式得到了自動化的大氣參數,這樣算法層面就可以自動化了,如果大家有其他算法,需要大氣參數的,
當然可以參考我的思路。
對於完整的溫度反演,我這里核心部分采用c++實現,部分代碼如下:
#include "mex.h" #include <cmath> using namespace std; /*地表溫度反演 * t: 大氣參數 * Lu: 大氣參數 * Ld: 大氣參數 * aima: 地表比輻射率 * band10_rad: 輻射定標后第10波段 */ void LST_Landsat8(double t, double Lu, double Ld, double *aima, double *band10_rad, int imglength, double *LST) { // double *LST = new double[imglength]; double *BlaRad = new double[imglength]; mexPrintf("大氣透過率參數: %f\n",t); mexPrintf("大氣向上透過率參數: %f\n",Lu); mexPrintf("大氣向下透過率參數: %f\n",Ld); mexPrintf("圖像大小為: %d\n",imglength); for (int i = 0; i < imglength; i++) { BlaRad[i] = (band10_rad[i] - Lu - t*(1 - aima[i]) * Ld) / (t*aima[i]); LST[i] = 1321.08 / log(774.89 / BlaRad[i] + 1) - 273; //mexPrintf("%f\n",BlaRad[i]); // mexPrintf("%f\n",LST[i]); } delete[] BlaRad; //釋放內存 //delete[] aima; //delete[] band10_rad; // return LST; } // c++入口主函數 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) //編譯代碼 { // 檢查輸入變量數量是否正確,否則報錯 if (nrhs != 5) mexErrMsgTxt("變量個數應當為5個...\n"); // 檢查輸出變量數量是否正確,否則報錯 if (nlhs > 1) mexErrMsgTxt("輸出變量個數不能超過1個...\n"); #define t prhs[0] // matlab傳過來的參數列表 #define Lu prhs[1] #define Ld prhs[2] #define aima prhs[3] #define b10_rad prhs[4] // #define LST_band plhs[0] // 傳回matlab的反演溫度波段 int M = mxGetM(aima); int N = mxGetN(aima); int imglength = M*N; //printf("%d", imglength); plhs[0] = mxCreateDoubleMatrix(M, N, mxREAL); double *t1; // 得到輸入參數指針變量 double *Lu1; double *Ld1; double *aima1; double *b10_rad1; double *lst_band; t1 = mxGetPr(t); // 得到輸入參數指針變量 Lu1 = mxGetPr(Lu); Ld1 = mxGetPr(Ld); aima1 = mxGetPr(aima); b10_rad1 = mxGetPr(b10_rad); /*mexPrintf("%f", *t1); mexPrintf("%f", *Lu1); mexPrintf("%f", *Ld1);*/ // lst_band = mxGetPr(LST_band); lst_band = mxGetPr(plhs[0]); // 執行c++溫度反演主函數 // 特別注意一下,c++混合編程時,輸出變量應當作為函數的輸入 // lst_band = LST_Landsat8(t1[0], Lu1[0], Ld1[0], aima1, b10_rad1, imglength, lst_band); LST_Landsat8(t1[0], Lu1[0], Ld1[0], aima1, b10_rad1, imglength, lst_band); }
上面是黑體輻射部分,之所以采用c++實現,是c++效率比較高,比matlab大概要快3倍左右。
全自動反演軟件:
2.實驗結果
實驗算法,我采用湖南省、重慶市夏季、冬季等L8影像,部分反演結果如下所示:
這是渲染后的結果,可以看到結果還是比較好的。
3.繼續改進
這里面我采用爬蟲技術曲線救國實現自動化批量溫度反演,后面我會繼續采用全自動的爬蟲技術獲取遙感影像,免去數據下載的痛苦,實現真正的一鍵大尺度溫度反演。交流qq:1044625113