其他分享
首页 > 其他分享> > 全自动Landsat影像温度反演软件开发

全自动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\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

 

标签:mexPrintf,软件开发,double,Landsat,反演,LST,imglength,aima
来源: https://www.cnblogs.com/wzp-749195/p/11000411.html