Google Earth Engine笔记-计算时间序列hurst指数 GEE

3 计算时间序列hurst指数

理论介绍

这篇http://club.excelhome.net/thread-1414154-1-1.html文章对hurst计算介绍非常详细，讲得非常详细。

文章图片

文章图片

【Google Earth Engine笔记-计算时间序列hurst指数】hurst指数的意义

文章图片

数据准备

var wuhan = ee.FeatureCollection("users/yp7454982/wuhan"); Map.centerObject(wuhan,8); var startDate=ee.Date('2000-01-01'); var endDate=ee.Date('2020-01-01'); var dataset=ee.ImageCollection("MODIS/006/MOD13Q1") .filterDate(startDate,endDate) .filterBounds(wuhan) .select("EVI") .map(function (image){return image.clip(wuhan).divide(10000).toFloat()}) print(dataset)

选择武汉市2000年到2020年的MODIS的MOD13Q1的16天时间序列产品，选择其中的EVI波段，并乘上scale因子，转换为浮点类型，进行裁剪。最终形成了457张时间序列影像

文章图片

定义子序列长度列表 r 1 , r 2 , . . . , r m r_1,r_2,...,r_m r1?,r2?,...,rm?

var sub_list_dist=4; var sub_list_len=ee.List.sequence(2,dataset.size().divide(2).int(),sub_list_dist);

为计算方便，定义子序列长度间隔为4,由于 r m r_m rm?要等于 n / 2 n/2 n/2的整数部分，因此有如上定义子序列长度列表

hurst指数 R / S R/S R/S计算函数

function hurst(time_series,sub_len){ var RS_list=sub_len.map(function (dist){ var index_seq=ee.List.sequence(0,time_series.size(),dist) var sub_R_S=ee.List.sequence(0,index_seq.size().subtract(2),1).map(function (i){ var sub_series=ee.ImageCollection(time_series.toList(time_series.size()) .slice(index_seq.get(ee.Number(i)),index_seq.get(ee.Number(i).add(1)))) var sub_series_mean=sub_series.mean() var sub_x_subtract_mean=sub_series.map(function (image){return image.subtract(sub_series_mean)}) //Calculate the deviation var first = ee.List([ // Rename the first image ee.Image(0)]); // This is a function to pass to Iterate(). // As anomaly images are computed, add them to the list. var accumulate = function(image,list) { // Get the latest cumulative anomaly image from the end of the list with // get(-1).Since the type of the list argument to the function is unknown, // it needs to be cast to a List.Since the return type of get() is unknown, // cast it to Image. var previous = ee.Image(ee.List(list).get(-1)); // Add the current anomaly to make a new cumulative anomaly image. var added = image.add(previous); // Return the list with the cumulative anomaly inserted. return ee.List(list).add(added); }; // Create an ImageCollection of cumulative anomaly images by iterating. // Since the return type of iterate is unknown, it needs to be cast to a List. var cum_deviation = ee.ImageCollection(ee.List(sub_x_subtract_mean.iterate(accumulate, first)).removeAll(first)); var range=cum_deviation.max().subtract(cum_deviation.min()) var standard_dev=sub_x_subtract_mean.map(function (image){return image.pow(ee.Image(2))}) .sum().divide(sub_x_subtract_mean.size().subtract(1)).sqrt() var R_S= range.divide(standard_dev).rename('RS').toFloat(); return R_S }) var R_S_mean=ee.ImageCollection(sub_R_S).mean().set('dist',dist); return R_S_mean }) var R_Smean_series=ee.ImageCollection(RS_list) return R_Smean_series; }

算法核心思想：代码的所有实现均根据计算hurst指数理论分部往下写，其中要注意的是利用了两个map循环，其中第一个对dist进行map循环是对其中子序列长度进行遍历，每进行一次遍历，需要对原始时间序列进行划分为若干组，然后再对每个组的时间序列进行整合成新的ImageCollection，这里必须要将其转为list来进行slice选择，再进行类型转换，对每一个组内计算离差，累计离差，累计离差这里要用到iterate的迭代的思想，求其累加的和，然后求极差，最终求得每组的 R / S R/S R/S然后再对所有组形成的 R / S R/S R/S形成的ImageCollection求mean,返回当前dist下的 R / S R/S R/S最终返回了一个再每个子序列长度不同下，计算的 R / S R/S R/S

hurst指数计算——最小二乘拟合

var RS=hurst(dataset,sub_list_len) print(RS) var RS_list=RS.toList(RS.size()); //transfrom List we can loop in i in[1,RS.size()] var lnRS_ln_i=ee.ImageCollection(ee.List.sequence(0,RS.size().subtract(1),1).map(function (i){ var lnRS_i=ee.Image(RS_list.get(i)).log().rename('ln_RS'); var ln_i=ee.Image(ee.Number(sub_list_len.get(i)).log()).rename('ln_r_i').toFloat(); return lnRS_i.addBands(ln_i) })); print(lnRS_ln_i); var hurst=lnRS_ln_i.select(['ln_r_i','ln_RS']).reduce(ee.Reducer.linearFit()).select('scale'); print(hurst)

异常值检查与数据可视化与导出

var viz={ min:0, max:1, palette:[ '040274','040281','0502a3','0502b8','0502ce','0502e6', '0602ff','235cb1','307ef3','269db1','30c8e2','32d3ef', '3be285','3ff38f','86e26f','3ae237','b5e22e','d6e21f', 'fff705','ffd611','ffb613','ff8b13','ff6e08','ff500d', 'ff0000','de0101','c21301','a71001','911003' ] }; Map.addLayer(hurst,viz,'hurst') var chart=ui.Chart.image.histogram(hurst, wuhan, 250) print(chart) Export.image.toDrive({ image: hurst.clip(wuhan), description: "hurst_wuhan", scale: 250, region: wuhan.geometry(), maxPixels: 1e13 }); var legend = ui.Panel({ style: { position: 'bottom-left', padding: '8px 15px' } }); // Create legend title var legendTitle = ui.Label({ value: 'scale', style: { fontWeight: 'bold', fontSize: '18px', margin: '0 0 4px 0', padding: '0' } }); // Add the title to the panel legend.add(legendTitle); // create the legend image var lon = ee.Image.pixelLonLat().select('latitude'); var gradient = lon.multiply((viz.max-viz.min)/100.0).add(viz.min); var legendImage = gradient.visualize(viz); // create text on top of legend var panel = ui.Panel({ widgets: [ ui.Label(viz['max']) ], }); legend.add(panel); // create thumbnail from the image var thumbnail = ui.Thumbnail({ image: legendImage, params: {bbox:'0,0,10,100', dimensions:'10x200'}, style: {padding: '1px', position: 'bottom-center'} }); // add the thumbnail to the legend legend.add(thumbnail); // create text on top of legend var panel = ui.Panel({ widgets: [ ui.Label(viz['min']) ], }); legend.add(panel); Map.add(legend);