LandTrendr算法在GEE上实现过程中注意事项

摘要：目前LandTrendr实现算法比较多了，最权威的还是俄勒冈州立大学的LT-GEE Guide。但是该算法很长，~1500行，如果只是普通实现LandTrendr，需要自己精简、优化。这个过程中，我们发现了一些容易出错地方，归纳如下：（1）LandTrendr结果的提取与使用；（2）影像值的反转与是否乘以-1，这关系到是提取gain和loss。

1.LandTrendr结果的提取与使用

没有看到以下解释之前，我对GEE中LandTrendr结果完全是懵的，不知道内在结果具体是什么。后面就知道是Band1的LandTrendr，Band2的RMSE，Band3不一定有，如果没有多的图层导入。本研究只用了NDVI做时序检测，所以只有Band1和Band2。对于Band1的LandTrendr，如图2所示，里面含有4个list，分别对应：时间年、NDVI原值、NDVI拟合值和是否是断点（1或0）。

图1 LT-GEE 的结果本质上是每个像素的列表的列表，这些列表描述分割，并可选地提供拟合的年度光谱数据 (FTV)。输出以至少包含 2 个波段的 GEE 图像形式提供，一个波段包含年度分割信息，另一个波段包含分割拟合的 RMSE。此外，如果 LT-GEE 的输入图像集合包含多个波段，则第一个波段之后的每个波段都将表示为光谱拟合的年度系列 (FTV)。

[
  [1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, ...] // Year list 
  [ 811,  809,  821,  813,  836,  834,  833,  818,  826,  820,  765,  827,  775, ...] // Source list
  [ 827,  825,  823,  821,  819,  817,  814,  812,  810,  808,  806,  804,  802, ...] // Fitted list
  [   1,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    0,    1, ...] // Is Vertex list
]

图2 LT-GEE的LandTrendr结果

    它包含 4 行，列数与时间序列中给定像素的年度观测值一样多。

• Row 1 is the observation year.
• Row 2 is the observation value corresponding to the year in row 1, it is equal to the first band in the input collection.
• Row 3 is the observation value corresponding to the year in row 1, fitted to line segments defined by breakpoint vertices identified in segmentation.
• Row 4 is a Boolean value indicating whether an observation was identified as a vertex.

     代码中是使用extractDisturbance函数来提取LandTrendr结果，其中第一步是提取断点，即得到扰动发生年图层，这里使用的是ee.Image.arraySlice(axis, start, end, step)函数，axis中0代表row方向，1代表column方向。所以这里lt.arraySlice(0, 3, 4)，起始index是0，就是行方向从第3行开始，到第4行结束，对应提取的就是Row 4。再利用掩膜功能arrayMask，1值（断点）保留，0值（非断点）去除。

// 选择只包含变化的顶点，将其分配到变量vertexMask get the vertex - yes(1)/no(0) dimension
//lt = lt.select('LandTrendr')
var vertexMask = lt.arraySlice(0, 3, 4); 
// 通过arrayMask方法获取该遮罩的值，并返回一个只包含变化的顶点的图像对象convert the 0's to masked
var vertices = lt.arrayMask(vertexMask); 

    得到断点后，需要对断点时间片段进行处理。代码中还是利用ee.Image.arraySlice()函数，不过此时axis轴变为1，即对一行值进行处理。left和right错开一位。

// 对时序片段进行处理，construct segment start and end point years and index values

  var left = vertices.arraySlice(1, 0, -1);    
// slice out the vertices as the start of segments，只保留左侧变化点V1,V2,V3
//-1表示去除最后一个点
  var right = vertices.arraySlice(1, 1, null); 
// slice out the vertices as the end of segments，只保留右侧变化点V2,V3,V4
//null表示一直到结束
  
  var startYear = left.arraySlice(0, 0, 1);    
// 左片段年，get year dimension of LT data from the segment start vertices，
  var startVal = left.arraySlice(0, 2, 3);     
// 左片段拟合值，get spectral index dimension of LT data from the segment start vertices
  var endYear = right.arraySlice(0, 0, 1);     
// 右片段年，get year dimension of LT data from the segment end vertices 
  var endVal = right.arraySlice(0, 2, 3);      
//右片段拟合值， get spectral index dimension of LT data from the segment end vertices
  
  var dur = endYear.subtract(startYear);       
// subtract the segment start year from the segment end year to calculate the duration of segments 
  var mag = endVal.subtract(startVal);         
// substract the segment start index value from the segment end index value to calculate the delta of segments 

图3 右片段减去左片段的示意图，这样就可得到3时间段（Y12，Y23，Y34）持续时间和变化程度。

    通过cat函数制作片段属性的图像。

  var distImg = ee.Image.cat([startYear.add(1), mag, dur, startVal])
                        .toArray(0); 
//制作片段属性的图像,concatenate segment start year, delta, duration, and starting spectral index value to an array 

2.影像值的反转与是否乘以-1，这关系到是提取gain和loss

     俄勒冈州立大学的LT-GEE Guide 中提到对于Image Collection的进行分割前的2个重要步骤：去云和光谱值方向（确保△VI为正时代表植被loss）。第一个步骤很好理解，但是第二个步骤其实并不是必备的，这就导致使用LandTrendr结果时需要乘以-1或1才能正确提取gain或loss。

    以我们研究中的NDVI为例。我们利用的最优合成年际影像直接进行分割，因此，植被loss就是负的△NDVI，植被gain是正的△NDVI。这里mag.multiply(-1)就是得到gain，mag.multiply(1)就是得到loss。

// sort the segments in the disturbance attribute image delta by spectral index change delta  
  var distImgSorted = distImg.arraySort(mag.multiply(-1));  // flip the delta around so that the greatest delta segment is first in order
                                                            //arraySort is ascending order, namely 1,2,3,4,5,6···
  // slice out the first (greatest) delta
  var tempDistImg = distImgSorted.arraySlice(1, 0, 1); //get the first segment in the sorted array