时间动态规整算法（Dynamic Time Warping，DTW）是一种常用到的时间序列分析方法，常用于时间序列分类、模式发现。
卫星影像时间序列分类的动态时间规整介绍：https://medium.com/soilwatch/dynamic-time-warping-for-satellite-image-time-series-classification-872d9e54b8d
分类结果如下所示：
本文分享了外网上查到的 GEE平台上使用哨兵数据（Sentinel-2）实现时间加权或时间约束动态时间规整 (TW/TC-DTW) 的模块，以及例子。
代码来源：https://github.com/wouellette/ee-dynamic-time-warping
时间加权 DTW (TW-DTW) 方法取自：Maus, V., Câmara, G., Cartaxo, R., Sanchez, A., Ramos, F. M., & De Queiroz, G. R. (2016). A time-weighted dynamic time warping method for land-use and land-cover mapping. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9(8), 3729-3739.
时间约束的 DTW (TC-DTW) 方法来自：Csillik, O., Belgiu, M., Asner, G. P., & Kelly, M. (2019). Object-based time-constrained dynamic time warping classification of crops using Sentinel-2. Remote sensing, 11(10), 1257.
矢量距离 (VD-TW) 方法取自：Teke, Mustafa, and Yasemin Y. Çetin. “Multi-year vector dynamic time warping-based crop mapping.” Journal of Applied Remote Sensing 15.1 (2021): 016517.
要使用DTW模块，需要将其复制粘贴到您的代码编辑器环境中，或者使用下行公开可用的脚本：
var DTW = require('users/soilwatch/functions:dtw.js');
调用的时候，像下面的例子这么用。
var training_data_list = DTW.prepareSignatures(reference_signatures,
CLASS_NAME,
key,
BAND_NO,
PATTERNS_LEN,
band_names);
GEE平台上实现时间加权或时间约束动态时间规整 (TW/TC-DTW) 的模块
exports.DTWDist = function(patterns_arr, timeseries_col, options){
timeseries_col = ee.ImageCollection(timeseries_col);
patterns_arr = ee.Array(patterns_arr);
var patterns_no = options.patterns_no || patterns_arr.length().get([0]);
var band_no = options.band_no || patterns_arr.length().get([1]).subtract(1);
var timeseries_len = options.timeseries_len || timeseries_col.size();
var patterns_len = options.patterns_len || patterns_arr.length().get([2]);
var constraint_type = options.constraint_type || 'time-weighted';
var weight_type = options.weight_type || 'logistic';
var distance_type = options.distance_type || 'euclidean';
var beta = options.beta || 50;
var alpha = options.alpha || 0.1;
var cost_weight;
var dis_arr;
var dis_mat;
var _distCalc = function(img, j, k){
var wrap = function(n, previous2){
n = ee.Number(n);
previous2 = ee.List(previous2);
var x1 = img.select(n.subtract(1)).toArray().toInt16();
var y1 = patterns_arr.get(ee.List([ee.Number(k).subtract(1), n.subtract(1), j.subtract(1)]));
if (distance_type === 'euclidean') {
dis_arr = x1.subtract(y1).pow(2);
} else if (distance_type === 'angular') {
var img_prev = timeseries_col.filter(ee.Filter.lte('doy', img.get('doy')))
.limit(2, 'doy', false).sort('doy').first();
var x2 = img_prev.select(n.subtract(1)).toArray().toInt16();
var y2 = patterns_arr.get(ee.List([ee.Number(k).subtract(1), n.subtract(1), j.subtract(2)]));
dis_arr = x1.multiply(y1).add(x2.multiply(y2))
.divide(x1.pow(2).add(x2.pow(2)).sqrt().multiply(y1.pow(2).add(y2.pow(2)).sqrt()))
.acos();
}
return dis_arr.add(previous2)
}
return wrap
}
var matrix = ee.List.sequence(1, ee.Number(timeseries_len).subtract(1)).map(function(i){
var matrix_tmp = ee.List.sequence(1, ee.Number(patterns_len).subtract(1)).map(function(j){
return ee.List([i, j]);
});
return matrix_tmp;
});
matrix = ee.List(matrix.iterate(function(x, previous){
return ee.List(previous).cat(ee.List(x));
}, ee.List([])));
var dtw_image_list = ee.List.sequence(1, patterns_no).map(function(k){
if (constraint_type === 'time-constrained') {
var dt_list = ee.List(timeseries_col.toList(timeseries_col.size()).map(function(img) {
img = ee.Image(img);
var dt_list = ee.List.sequence(1, patterns_len).map(function(j) {
j = ee.Number(j);
var t1 = ee.Number(img.get('doy'));
var t2 = patterns_arr.get(ee.List([ee.Number(k).subtract(1), -1, j.subtract(1)]));
var time_arr = t1.subtract(t2).abs();
return ee.Feature(null, {
'dt': time_arr,
'i': t1,
'j': j
});
});
return dt_list;
}));
dt_list = dt_list.flatten();
dis_mat = timeseries_col.toList(timeseries_col.size()).map(function(img){
img = ee.Image(img);
var patterns_tmp = ee.FeatureCollection(dt_list)
.filter(ee.Filter.and(ee.Filter.eq('i', ee.Number(img.get('doy'))),
ee.Filter.lte('dt', beta)))
.aggregate_array('j');
var dis_list0 = patterns_tmp.map(function(j){
j = ee.Number(j);
var dis_sum = ee.Image(ee.List.sequence(1, ee.Number(band_no)).iterate(_distCalc(img, j, k),
ee.Image(0)));
if (distance_type === 'angular') {
dis_sum = dis_sum.multiply(t1.neq(timeseries_col.first().get('doy')));
}
return dis_sum.sqrt().set('j', j);
});
patterns_tmp = ee.FeatureCollection(dt_list)
.filter(ee.Filter.and(ee.Filter.eq('i', ee.Number(img.get('doy'))),
ee.Filter.gt('dt', beta)))
.aggregate_array('j');
var dis_list = patterns_tmp.map(function(j){
j = ee.Number(j);
return ee.Image(1e6).set('j', j);
});
dis_list = dis_list0.cat(dis_list);
dis_list = ee.ImageCollection(dis_list).sort('j').toList(dis_list.length());
return dis_list;
});
} else if (constraint_type === 'time-weighted') {
dis_mat = timeseries_col.toList(timeseries_col.size()).map(function(img){
img = ee.Image(img);
var dis_list = ee.List.sequence(1, patterns_len).map(function(j){
j = ee.Number(j);
var t1 = ee.Number(img.get('doy'));
var t2 = patterns_arr.get(ee.List([ee.Number(k).subtract(1), -1, j.subtract(1)]));
var time_arr = t1.subtract(t2).abs();
var dis_sum = ee.Image(ee.List.sequence(1, ee.Number(band_no)).iterate(_distCalc(img, j, k),
ee.Image(0)));
if (distance_type === 'angular') {
dis_sum = dis_sum.multiply(t1.neq(timeseries_col.first().get('doy')));
}
if (weight_type === 'logistic') {
cost_weight = ee.Image(1)
.divide(ee.Image(1).add(ee.Image(alpha)
.multiply(time_arr.subtract(beta))).exp());
} else if (weight_type === 'linear') {
cost_weight = ee.Image(alpha).multiply(time_arr).add(beta);
}
return dis_sum.sqrt().add(cost_weight);
});
return dis_list;
});
}
var D_mat = ee.List([]);
var dis = ee.List(dis_mat.get(0)).get(0);
var d_mat = D_mat.add(dis);
d_mat = ee.List(ee.List.sequence(1, ee.Number(patterns_len).subtract(1)).iterate(function(j, previous){
j = ee.Number(j);
previous = ee.List(previous);
var dis = ee.Image(previous.get(j.subtract(1))).add(ee.List(dis_mat.get(0)).get(j));
return previous.add(dis);
}, d_mat));
D_mat = D_mat.add(d_mat);
D_mat = ee.List(ee.List.sequence(1, ee.Number(timeseries_len).subtract(1)).iterate(function(i, previous){
i = ee.Number(i);
previous = ee.List(previous);
var dis = ee.Image(ee.List(previous.get(i.subtract(1))).get(0)).add(ee.List(dis_mat.get(i)).get(0));
return previous.add(ee.List(previous.get(i.subtract(1))).set(0, dis));
}, D_mat));
D_mat = ee.List(ee.List.sequence(0, matrix.length().subtract(1)).iterate(function(x, previous){
var i = ee.Number(ee.List(matrix.get(x)).get(0));
var j = ee.Number(ee.List(matrix.get(x)).get(1));
previous = ee.List(previous);
var dis = ee.Image(ee.List(previous.get(i.subtract(1))).get(j)).min(ee.List(previous.get(i)).get(j.subtract(1)))
.min(ee.List(previous.get(i.subtract(1))).get(j.subtract(1))).add(ee.List(dis_mat.get(i)).get(j));
return previous.set(i, ee.List(previous.get(i)).set(j, dis));
}, D_mat));
return ee.Image(ee.List(D_mat.get(-1)).get(-1))
.arrayProject([0])
.arrayFlatten([['DTW']]);
});
return ee.ImageCollection(dtw_image_list).min().toUint16();
};
exports.prepareSignatures = function(signatures, class_name, class_no, band_no, patterns_len, band_names){
var train_points = signatures.filter(ee.Filter.eq(class_name, class_no)).select(band_names);
var feature_list = train_points.toList(train_points.size());
var table_points_list = feature_list.map(function (feat){
return ee.List(band_names).map(function (band){return ee.Feature(feat).get(band)});
});
return ee.Array(table_points_list).reshape([-1, band_no+1, patterns_len]).toInt16();
};
根据上述模块实现TW-DTW（时间加权的动态时间规整）分类活跃耕地/弃耕地的例子：
var palettes = require('users/gena/packages:palettes');
var wrapper = require('users/adugnagirma/gee_s1_ard:wrapper');
var S2Masks = require('users/soilwatch/soilErosionApp:s2_masks.js');
var composites = require('users/soilwatch/soilErosionApp:composites.js');
var DTW = require('users/soilwatch/functions:dtw.js');
var CLASS_NAME = 'lc_class';
var AGG_INTERVAL = 30;
var TIMESERIES_LEN = 6;
var PATTERNS_LEN = 6;
var CLASS_NO = 7;
var S2_BAND_LIST = ['B2', 'B3', 'B11', 'B12', 'ndvi'];
var S1_BAND_LIST = ['VV', 'VH'];
var BAND_NO = S1_BAND_LIST.concat(S2_BAND_LIST).length;
var DOY_BAND = 'doy';
var BETA = 50;
var ALPHA = 0.1;
var not_water = ee.Image("JRC/GSW1_2/GlobalSurfaceWater").select('max_extent').eq(0);
var dem = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2").select("DSM");
dem = dem.mosaic().setDefaultProjection(dem.first().select(0).projection());
var slope = ee.Terrain.slope(dem);
var dem_mask = dem.lt(3600);
var slope_mask = slope.lt(30);
var crop_mask = dem_mask.and(slope_mask);
var addNDVI = function(img){
return img.addBands(img.normalizedDifference(['B8','B4']).multiply(10000).toInt16().rename('ndvi'));
};
var adm0_name = 'Sudan';
var adm1_name = 'Sennar';
var adm2_name = 'Sennar';
var year_dict = {
'2019': 'COPERNICUS/S2',
'2018': 'COPERNICUS/S2',
'2017': 'COPERNICUS/S2'
};
var county = ee.Feature(
ee.FeatureCollection(ee.FeatureCollection("FAO/GAUL/2015/level2")
.filter(ee.Filter.and(ee.Filter.equals('ADM0_NAME', adm0_name),
ee.Filter.equals('ADM2_NAME', adm2_name)))
).first());
Map.centerObject(county.geometry());
Map.layers().reset([ui.Map.Layer(county, {}, adm2_name)]);
var DTWClassification = function(year, collection_type){
var date_range = ee.Dictionary({'start': year + '-07-01', 'end': year + '-12-30'});
var s2_cl = S2Masks.loadImageCollection(collection_type, date_range, county.geometry());
var masked_collection = s2_cl
.filterDate(date_range.get('start'), date_range.get('end'))
.map(S2Masks.addCloudShadowMask(not_water, 1e4))
.map(S2Masks.applyCloudShadowMask)
.map(addNDVI);
var time_intervals = composites.extractTimeRanges(date_range.get('start'), date_range.get('end'), 30);
var s2_stack = composites.harmonizedTS(masked_collection, S2_BAND_LIST, time_intervals, {agg_type: 'geomedian'});
var parameter = {
START_DATE: date_range.get('start'),
STOP_DATE: date_range.get('end'),
POLARIZATION:'VVVH',
ORBIT : 'DESCENDING',
GEOMETRY: county.geometry(),
APPLY_ADDITIONAL_BORDER_NOISE_CORRECTION: true,
APPLY_SPECKLE_FILTERING: true,
SPECKLE_FILTER_FRAMEWORK: 'MULTI',
SPECKLE_FILTER: 'LEE',
SPECKLE_FILTER_KERNEL_SIZE: 9,
SPECKLE_FILTER_NR_OF_IMAGES: 10,
APPLY_TERRAIN_FLATTENING: true,
DEM: dem,
TERRAIN_FLATTENING_MODEL: 'VOLUME',
TERRAIN_FLATTENING_ADDITIONAL_LAYOVER_SHADOW_BUFFER: 0,
FORMAT : 'DB',
CLIP_TO_ROI: false,
SAVE_ASSETS: false
}
var s1_ts = wrapper.s1_preproc(parameter)[1]
.map(function(image){return image.multiply(1e4).toInt16()
.set({'system:time_start': image.get('system:time_start')})});
var s1_stack = composites.harmonizedTS(s1_ts, S1_BAND_LIST, time_intervals, {agg_type: 'geomedian'});
var filter = ee.Filter.equals({
leftField: 'system:time_start',
rightField: 'system:time_start'
});
var simpleJoin = ee.Join.inner();
var innerJoin = ee.ImageCollection(simpleJoin.apply(s1_stack, s2_stack, filter))
var joined = innerJoin.map(function(feature) {
return ee.Image.cat(feature.get('primary'), feature.get('secondary'));
});
joined = joined.map(function(image){
var currentDate = ee.Date(image.get('system:time_start'));
var meanImage = joined.filterDate(currentDate.advance(-AGG_INTERVAL-1, 'day'),
currentDate.advance(AGG_INTERVAL+1, 'day')).mean();
var ddiff = currentDate.difference(ee.Date(ee.String(date_range.get('start')))
.format('YYYY').cat('-01-01'),
'day');
return meanImage.where(image, image).unmask(0)
.addBands(ee.Image(ddiff).rename('doy').toInt16())
.set({'doy': ddiff.toInt16()})
.copyProperties(image, ['system:time_start']);
}).sort('system:time_start');
var s1s2_stack = ee.Image(joined.iterate(function(image, previous){return ee.Image(previous).addBands(image)}, ee.Image([])))
.select(ee.List(S1_BAND_LIST.concat(S2_BAND_LIST)).add(DOY_BAND).map(function(band){return ee.String(band).cat('.*')}));
var band_names = s1s2_stack.bandNames();
var reference_signatures = reference_signatures || s1s2_stack
.sampleRegions({
collection: signatures,
properties: [CLASS_NAME],
scale : 10,
geometries: true
});
var dtw_min_dist = function(key, val){
key = ee.Number.parse(key);
var training_data_list = DTW.prepareSignatures(reference_signatures,
CLASS_NAME,
key,
BAND_NO,
PATTERNS_LEN,
band_names);
return ee.ImageCollection(DTW.DTWDist(training_data_list,
joined.select('[^'+DOY_BAND+'].*'),
{patterns_no: val,
band_no: BAND_NO,
timeseries_len: TIMESERIES_LEN,
patterns_len: PATTERNS_LEN,
constraint_type: 'time-weighted',
beta: BETA,
alpha: ALPHA
})
).min()
.rename('dtw')
.addBands(ee.Image(key).toByte().rename('band'));
};
var dtw_image_list = reference_signatures_agg.map(dtw_min_dist);
var array = ee.ImageCollection(dtw_image_list.values()).toArray();
var axes = {image:0, band:1};
var sort = array.arraySlice(axes.band, 0, 1);
var sorted = array.arraySort(sort);
var values = sorted.arraySlice(axes.image, 0, 1);
var min = values.arrayProject([axes.band]).arrayFlatten([['dtw', 'band']]);
var dtw_score = min.select(0).rename('score_' + year);
var dtw_class = min.select(1).rename('classification_' + year);
return [dtw_class.addBands(dtw_score), reference_signatures, s1s2_stack];
};
var signatures = ee.FeatureCollection('users/soilwatch/Sudan/SennarSignatures')
.filterBounds(county.geometry());
var reference_signatures_agg = signatures.aggregate_histogram(CLASS_NAME);
print('Number of reference signatures per land cover class:')
print(reference_signatures_agg);
var classification_names = ['built-up',
'water',
'trees',
'bare soil',
'active cropland',
'rangelands',
'wetland',
'abandoned/long-term fallow cropland (> 3 years)',
'short-term fallow cropland (
];
var classification_palette = ['#d63000',
'blue',
'#4d8b22',
'grey',
'ebff65',
'#98ff00',
'purple',
'#00ce6f',
'#FFA33F'
];
var score_palette = palettes.colorbrewer.RdYlGn[9].reverse();
var dtw_outputs = DTWClassification('2020', 'COPERNICUS/S2');
var dtw = dtw_outputs[0];
var reference_signatures = dtw_outputs[1];
var s1s2_stack = dtw_outputs[2];
var s1s2_list = ee.List([s1s2_stack]);
Object.keys(year_dict).forEach(function(i) {
var dtw_outputs = DTWClassification(i, year_dict[i])
dtw = dtw.addBands(dtw_outputs[0]);
s1s2_list = s1s2_list.add(dtw_outputs[2]);
});
var imageVisParam = {bands: ["ndvi_5", "ndvi_3", "ndvi_1"],
gamma: 1,
max: 7000,
min: 1000,
opacity: 1
};
Map.addLayer(s1s2_stack.clip(county.geometry()), imageVisParam, 'ndvi stack 2020');
imageVisParam['bands'] = ["VV_5", "VV_3", "VV_1"];
Map.addLayer(s1s2_stack.clip(county.geometry()), imageVisParam, 'VV stack 2020');
var dtw = ee.Image('users/soilwatch/Sudan/dtw_sennar_s1s2_2017_2020').clip(county.geometry());
var dtw_score = dtw.select('score_2020');
var dtw_class = dtw.select('classification_2020');
var vito_lulc_crop = ee.ImageCollection("ESA/WorldCover/v100").filterBounds(county.geometry()).mosaic().eq(40);
Map.addLayer(dtw_score.clip(county.geometry()),
{palette: score_palette, min: 0, max: 20000},
'DTW dissimilarity score (30 days signature, 30 days images)');
Map.addLayer(dtw_class.updateMask(crop_mask).clip(county.geometry()),
{palette: classification_palette.slice(0, classification_palette.length-2), min: 1, max: CLASS_NO},
'DTW classification (30 days signature, 30 days images)');
var dtw_class_strat = dtw_class.where(dtw_class.eq(6)
.and(dtw.select('classification_2019').eq(6))
.and(dtw.select('classification_2018').eq(6))
.and(dtw.select('classification_2017').eq(6))
.and(vito_lulc_crop),
ee.Image(8))
.where(dtw_class.eq(6)
.and(dtw.select('classification_2019').eq(5)
.or(dtw.select('classification_2018').eq(5))
.or(dtw.select('classification_2017').eq(5))
).and(vito_lulc_crop),
ee.Image(9));
Map.addLayer(dtw_class_strat.updateMask(crop_mask).clip(county.geometry()),
{palette: classification_palette, min: 1, max: CLASS_NO+2},
'DTW classification: active/abandoned cropland distinction');
var area_histogram = ee.Dictionary(dtw_class_strat.updateMask(crop_mask).reduceRegion(
{reducer: ee.Reducer.frequencyHistogram(),
geometry: county.geometry(),
scale: 100,
maxPixels:1e13,
tileScale: 4
}).get('classification_2020')).values();
var area_list = ee.List([]);
for (var i = 0; i CLASS_NO+1; i++) {
area_list = area_list.add(ee.Feature(null, {
'area': area_histogram.get(i),
'class': classification_names[i]
}));
}
var options = {title: 'Land Cover Classes Distribution',
colors: classification_palette,
sliceVisibilityThreshold: 0
};
var area_chart = ui.Chart.feature.byFeature({
features: ee.FeatureCollection(area_list),
xProperty: 'class',
yProperties: ['area']
})
.setChartType('PieChart')
.setOptions(options);
print(area_chart);
var video_args = {
'dimensions': 2500,
'region': county.geometry(),
'framesPerSecond': 1,
'crs': 'EPSG:4326',
'min': 1,
'max': 9,
'palette': classification_palette
}
var generateGIF = function(img, band_name){
print(band_name+' GIF')
print(img.select(band_name).getVideoThumbURL(video_args));
}
generateGIF(ee.ImageCollection(ee.List([dtw.select('classification_2020').rename('classification')])
.add(dtw_class_strat.rename('classification'))), 'classification');
Export.video.toDrive({
collection: ee.ImageCollection(s1s2_list.map(function(img){return ee.Image(img).divide(100).toByte()}))
.select(["ndvi_5", "ndvi_3", "ndvi_1"]),
description:'temporal ndvi rgb',
region: county.geometry(),
scale: 100,
crs:'EPSG:4326',
maxPixels: 1e13
});
var setPointProperties = function(f){
var class_val = f.get("lc_class");
var mapDisplayColors = ee.List(classification_palette);
return f.set({style: {color: mapDisplayColors.get(ee.Number(class_val).subtract(1))}})
}
var styled_td = signatures.map(setPointProperties);
Map.addLayer(styled_td.style({styleProperty: "style"}), {}, 'crop type sample points');
Export.image.toAsset({
image: dtw.clip(county.geometry()),
description:'dtw_sennar_s1s2_2017_2020',
assetId: 'users/soilwatch/Sudan/dtw_sennar_s1s2_2017_2020',
region: county.geometry(),
crs: 'EPSG:4326',
scale: 10,
maxPixels:1e13
});
var legend = ui.Panel({
style: {
position: 'bottom-left',
padding: '8px 15px'
}
});
var legendTitle = ui.Label({
value: 'Legend',
style: {
fontWeight: 'bold',
fontSize: '18px',
margin: '0 0 4px 0',
padding: '0'
}
});
legend.add(legendTitle);
var makeRow = function(color, name) {
var colorBox = ui.Label({
style: {
backgroundColor: color,
padding: '8px',
margin: '0 0 4px 0'
}
});
var description = ui.Label({
value: name,
style: {margin: '0 0 4px 6px'}
});
return ui.Panel({
widgets: [colorBox, description],
layout: ui.Panel.Layout.Flow('horizontal')
});
};
for (var i = 0; i CLASS_NO+1; i++) {
legend.add(makeRow(classification_palette[i], classification_names[i]));
}
Map.add(legend);
Original: https://blog.csdn.net/qq_35591253/article/details/125710307
Author: 养乐多
Title: GEE：DTW（Dynamic Time Warping）动态时间规整，Sentinel-2 时间序列分类
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/666952/
转载文章受原作者版权保护。转载请注明原作者出处！