【Google Earth Engine线性拟合的显著性问题】

Google Earth Engine在进行统计分析时的不足

在遥感大数据时代，我们面临着海量的遥感数据——长时序、多平台、多传感器、多分辨率。GEE作为遥感云计算平台给研究人员在收集、处理数据上提供了很大的便利。在进行时序分析时，GEE也提供了诸如相关性、线性拟合等常用的方法。然而，GEE提供的线性拟合函数并不提供结果的显著性，因此在分析时存在一些麻烦。关于这个问题，也有很多同事私信我，本文在这里提供一个使用GEE进行显著性分析的思路，希望和大家一起讨论学习。

显著性检验方法（以小样本为例）

本人研究的时间范围为2000-2020年，时间分辨率为一年，因此每个pixel上有21个样本，在pixel尺度上开展线性拟合则需要使用t检验方式（t检验适用于样本数小于30）。t检验具体公式可以参照https://blog.csdn.net/semine_shen/article/details/80061594，本文使用GEE （JS）实现。

线性拟合第一步：增加年份波段

本人是做2000-2020年度的时间序列拟合，因此需要给每年的image增加一个年份波段用于进行线性拟合。本人自写了两个函数，AddTime2000to2020是增加年份属性，addyear则为增加名为year的波段

function AddTime2000to2020(IC){
    var timeList = ee.List([946656000000,
                            978278400000,
                            1009814400000,
                            1041350400000,
                            1072886400000,
                            1104508800000,
                            1136044800000,
                            1167580800000,
                            1199116800000,
                            1230739200000,
                            1262275200000,
                            1293811200000,
                            1325347200000,
                            1356969600000,
                            1388505600000,
                            1420041600000,
                            1451577600000,
                            1483200000000,
                            1514736000000,
                            1546272000000,
                            1577808000000])
    var idx = ee.Number(0)
    var iniDict = {
      index:idx,
      ic:ee.List([])
    }
    
    var Dp = ee.Dictionary(ee.ImageCollection(IC).iterate(setDate,iniDict));
    return Dp;
    
    function setDate(image,dict){
      var img = ee.Image(image)
      dict = ee.Dictionary(dict)
      
      var index = ee.Number(dict.get("index"))
      var starttime = ee.Number(timeList.get(index)).add(ee.Number(28800000))
      img = img.set("system:time_start",starttime)
      
      var new_ic = ee.List(dict.get("ic")).add(img)
      var new_idx = index.add(1)
      var new_dict = {
        index:new_idx,
        ic:new_ic
      }
      return new_dict
    }
}

function addyear(image){
   var image_date = ee.Date(image.date());
   var image_doy = ee.Number.parse(image_date.format('YYYY'));
   image = image.set("year",image_doy)
   return image.addBands(ee.Image(image_doy).rename('year').toInt())  
}

其中，AddTime2000to2020中timeList中的值为2000-2020年每年1月1日对应的距离1977年的毫秒数，因此可以通过修改这些毫秒数来给任意时间段的年度序列增加年份信息，本人写了一个python代码来获取期望时间段的数据，下面展示获取1985-2019年每年的年份信息的代码，只需要更改range内的时间范围即可：

import time, datetime

from numpy import long

timeStr = "-01-01 00:00:00"
for i in range(1985,2000):
    now = str(i) + timeStr
    time1 = datetime.datetime.strptime(now, "%Y-%m-%d %H:%M:%S")
    secondsFrom1970 = time.mktime(time1.timetuple())
    print(long(secondsFrom1970 * 1000))

只需使用下面一行代码便可以给指定的ImageCollection中的每年的Image增加一个对应的年份波段，只需修改下面代码中的evitgs为自己的ImageCollection即可：

evitgs = ee.ImageCollection(ee.List(ee.Dictionary(AddTime2000to2020(evitgs)).get("ic")))

线性拟合第二步： 线性拟合函数（以SenSlope为例）

senslope函数X在前，Y在后，select可以分别获取斜率和截距：

 var sen = evitgs.select(["year","evitgs"]).reduce(ee.Reducer.sensSlope()).select("slope")
  var offset = evitgs.select(["year","evitgs"]).reduce(ee.Reducer.sensSlope()).select("offset")

线性拟合第三步： 显著性分析（以t检验为例）

//x平均值
var x_ = evitgs.select("year").reduce(ee.Reducer.mean())
//y平均值
var y_ = evitgs.select("evitgs").reduce(ee.Reducer.mean())
//x平方
var x2 = evitgs.select("year").map(function(img){
      return img.subtract(x_).multiply(img.subtract(x_))
}).reduce(ee.Reducer.sum())
//y平方
var y2 = evitgs.select("evitgs").map(function(img){
      return img.subtract(y_).multiply(img.subtract(y_))
}).reduce(ee.Reducer.sum())
// (x-x_mean) * (y-y_mean)
var xy = evitgs.map(function(img){
            var x = img.select("year")
            var y = img.select("evitgs")
            var xy = x.subtract(x_).multiply(y.subtract(y_))
            return xy
}).reduce(ee.Reducer.sum())

var bottom2 = evitgs.map(function(img){
      var temp = img.select("evitgs").subtract(offset).subtract(sen.multiply(img.select("year")))
      return temp.multiply(temp)
 }).reduce(ee.Reducer.sum()).divide(n.subtract(2)).sqrt()
 
var t_linear = sen.multiply(x2.sqrt()).divide(bottom2)

我做的是双尾t检验，就对结果取了平均值，然后根据样本数计算自由度，通过查表来得到p<0.05时对应的t绝对值对应的阈值

t_linear = t_linear.abs()
t_linear.gt(1.729)

这样就可以得到一个0-1值的mask，再使用updateMask即可得到p<0.05对应的范围