基于python的Hurst计算预测未来发展趋势（长时序栅格影像）

news2024/11/16 0:48:57

1.Hurst指数反映了时间序列长期记忆性的程度，即过去的信息对未来的影响程度。Hurst指数的取值范围为0到1之间，当Hurst指数等于0.5时，时间序列被认为是一种随机漫步，即具有随机性；当Hurst指数大于0.5时，时间序列显示出长期正相关性，表明趋势在未来可能持续；当Hurst指数小于0.5时，时间序列显示出长期负相关性，表明趋势在未来可能反转。

2.下面是一个例子，计算下面22幅影像的Hurst。
在这里插入图片描述
计算结果如下：

3.理论不再进一步讲解，直接上代码，有基础的可以自己调试，小白可以使用本人编好的exe程序，链接在下面。

import numpy as np
from osgeo import gdal
import os

def write(file_name, image, projection,geotransform,x_size,y_size):
    dtype = gdal.GDT_Float32
    # 数据格式
    driver = gdal.GetDriverByName('GTIFF')
    # 创建数据,设置文件路径及名称
    new_ds = driver.Create(file_name, x_size, y_size, 1, dtype)
    # 设置投影信息及6参数
    new_ds.SetGeoTransform(geotransform)
    new_ds.SetProjection(projection)
    # 将值写入new_ds中
    new_ds.GetRasterBand(1).WriteArray(image)
    # 把缓存数据写入磁盘
    new_ds.FlushCache()

    del new_ds

def Hurst(x):
    # x为numpy数组
    n = x.shape[0]
    t = np.zeros(n - 1)  # t为时间序列的差分
    for i in range(n - 1):
        t[i] = x[i + 1] - x[i]
    mt = np.zeros(n - 1)  # mt为均值序列,i为索引,i+1表示序列从1开始
    for i in range(n - 1):
        mt[i] = np.sum(t[0:i + 1]) / (i + 1)

    # Step3累积离差和极差,r为极差
    r = []
    for i in np.arange(1, n):  # i为tao
        cha = []
        for j in np.arange(1, i + 1):
            if i == 1:
                cha.append(t[j - 1] - mt[i - 1])
            if i > 1:
                if j == 1:
                    cha.append(t[j - 1] - mt[i - 1])
                if j > 1:
                    cha.append(cha[j - 2] + t[j - 1] - mt[i - 1])
        r.append(np.max(cha) - np.min(cha))
    s = []
    for i in np.arange(1, n):
        ss = []
        for j in np.arange(1, i + 1):
            ss.append((t[j - 1] - mt[i - 1]) ** 2)
        s.append(np.sqrt(np.sum(ss) / i))
    r = np.array(r)
    s = np.array(s)
    xdata = np.log(np.arange(2, n))
    ydata = np.log(r[1:] / s[1:])

    h, b = np.polyfit(xdata, ydata, 1)
    return h


def main(path1,result_path):
    filepaths = []
    for file in os.listdir(path1):
        filepath1 = os.path.join(path1, file)
        filepaths.append(filepath1)

    # 获取影像数量
    num_images = len(filepaths)
    # 读取影像数据
    img1 = gdal.Open(filepaths[0])
    # 获取影像的投影，高度和宽度
    transform1 = img1.GetGeoTransform()
    proj = img1.GetProjection()
    height1 = img1.RasterYSize
    width1 = img1.RasterXSize
    array1 = img1.ReadAsArray(0, 0, width1, height1)
    del img1

    # 读取所有影像
    for path1 in filepaths[1:]:
        if path1[-3:] == 'tif':

            img2 = gdal.Open(path1)
            array2 = img2.ReadAsArray(0, 0, width1, height1)
            array1 = np.vstack((array1, array2))
            del img2

    array1 = array1.reshape((num_images, height1, width1))
    # 输出矩阵，无值区用nan填充
    h_array = np.full([height1, width1], np.nan)

    # 只有有值的区域才进行计算
    c1 = np.isnan(array1)
    sum_array1 = np.sum(c1, axis=0)
    nan_positions = np.where(sum_array1 == num_images)
    positions = np.where(sum_array1<10)
    for i in range(len(positions[0])):
        x = positions[0][i]
        y = positions[1][i]
        hurst_list1 = array1[:, x, y]
        hurst_list1=hurst_list1[~np.isnan(hurst_list1)]
        h = Hurst(hurst_list1)
        h_array[x, y] = h

    write(result_path, h_array, proj, transform1, width1, height1)
if __name__ == "__main__":
    path1 = r"F:\rsei-2"
    result_path = r"F:\jieguo\h.tif"
    main(path1,result_path)