MATLAB与Python：长时间遥感图像变异系数计算

变异系数简介：

变异系数是衡量观测序列值变异程度的一个统计量，可以很好地反映空间数据在时间序列上变化的差异程度，评价数据时间序列的稳定性。计算公式为：

计算指标的均值：

计算指标的标准差：

计算变异系数：

Vj 越大，数据波动越大。

应用于长时间序列遥感影像：就是把多幅影像叠加在一起，逐像元构建一组时间序列，然后计算它的变异系数。

从x年到x+n-1年共n年栅格数据

一般情况下，将区域CV划分五个等级，来表征区域数据稳定情况：

MATLAB代码实现：

[a,R]=geotiffread(&#39;E:\test\pre\2001year_mean.tif&#39;);%先导入某个图像的投影信息，为后续图像输出做准确
info=geotiffinfo(&#39;E:\test\pre\2001year_mean.tif&#39;);
[m,n]=size(a);
years=5; %表示有多少年份需要做回归
data=zeros(m*n,years);
k=1;
for year=2001:2005 %起始年份
    file=[&#39;E:\test\pre\&#39;,int2str(year),&#39;year_mean.tif&#39;];%注意自己的名字形式，这里使用的名字是2001year_mean.tif，根据这个可修改
    bz=importdata(file);
    bz=reshape(bz,m*n,1);
    data(:,k)=bz;
    k=k+1;
end
    bya=zeros(m,n);
for i=1:length(data)
    bz=data(i,:);
    bz=bz&#39;;
    if min(bz)&gt;0 %此处根据数据实际最小值范围更改
        %计算变异系数
    A=mean(bz) %求每列平均值
    S=std(bz,1)  %求每列方差
    V=S./A %变异系数

    bya(i)=V;
    end
end
name1=&#39;E:\test\pre\变异系数.tif&#39;;
geotiffwrite(name1,bya,R,&#39;GeoKeyDirectoryTag&#39;,info.GeoTIFFTags.GeoKeyDirectoryTag);

python代码实现

from osgeo import gdal
import numpy as np
import copy
# 计算变异系数
images = [
    r&#39;E:\test\pre\2001year_mean.tif&#39;,
    r&#39;E:\test\pre\2002year_mean.tif&#39;,
    r&#39;E:\test\pre\2003year_mean.tif&#39;,
    r&#39;E:\test\pre\2004year_mean.tif&#39;,
    r&#39;E:\test\pre\2005year_mean.tif&#39;,
]
def coefficient_of_variation(data): # 变异系数
    mean = np.mean(data) # 计算平均值
    std = np.std(data, ddof=0) # 计算标准差
    cv = std/mean
    return cv
# 栅格图像组计算变异系数
def CV(images, outpath):
    images_pixels = [] # 存放多个图像像元矩阵的空数组
    for image in images:
        tif = str(image)
        open_tif = gdal.Open(tif) # 打开栅格图像
        band = open_tif.ReadAsArray() # 获取波段的矩阵
        images_pixels.append(band) # 把图像的矩阵加入到数组中
    CV = copy.deepcopy(images_pixels[0]) # 获取一个矩阵作为要写入的模板
    # 读取像元，并计算变异系数
    for i in range(len(CV)):
        for j in range(len(CV[1])):
            CV_data = [] # 存放ij坐标下像元值的数组，以计算变异系数
            for px in range(len(images)): # 遍历多个图像下ij坐标的像元值
                CV_data.append(images_pixels[px][i][j]) # 同一坐标的多点加入数组，以计算变异系数
            CV_value = coefficient_of_variation(CV_data) # 变异系数计算
            CV[i][j] = CV_value # 写入该坐标下的变异系数
    # 保存栅格图像
    gtiff_driver = gdal.GetDriverByName(&#39;GTiff&#39;)
    out_tif = gtiff_driver.Create(outpath, CV.shape[1], CV.shape[0], 1, gdal.GDT_Float32)
    # 将数据坐标投影设置为原始坐标投影
    out_tif.SetProjection(open_tif.GetProjection())
    out_tif.SetGeoTransform(open_tif.GetGeoTransform())
    out_band = out_tif.GetRasterBand(1)
    out_band.WriteArray(CV)
    out_band.FlushCache()
    print(&#39;栅格图像组变异系数计算完成&#39;)
CV(images, r&#39;E:\test\pre\cv.tif&#39;)