全自动Landsat影像温度反演软件开发

        许久没有更新遥感类软件开发了，都有点生疏了，这一次我带来了一个老的算法，新的东西， 为什么这么说呢，我们知道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
",t);
    mexPrintf("大气向上透过率参数: %f
",Lu);
    mexPrintf("大气向下透过率参数: %f
",Ld);
    mexPrintf("图像大小为: %d
",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
",BlaRad[i]);
//         mexPrintf("%f
",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个...
");
    
    // 检查输出变量数量是否正确，否则报错
    if (nlhs > 1)
        mexErrMsgTxt("输出变量个数不能超过1个...
");
    
    #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