Google Earth Engine——利用Sentinel-1数据（db值）进行长时序均值加载

// 导入哨兵1号数据

var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')

.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))

.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))

.filter(ee.Filter.eq('instrumentMode', 'IW'))

.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))

.filterBounds(roi)

.filterDate('2014-01-01','2023-12-31')

.select(["VV","VH"])

//转化为DB值

function toNatural(img) {

return ee.Image(10.0).pow(img.select(0).divide(10.0));

}

//Function to convert to dB

function toDB(img) {

return ee.Image(img).log10().multiply(10.0);

}

//应用 SNAP 3.0 S1TBX 编码的精制李氏斑点滤波器：

///senbox-org/s1tbx/blob/master/s1tbx-op-sar-processing/src/main/java/org/esa/s1tbx/sar/gpf/filtering/SpeckleFilters/

//改编自 Guido Lemoine

function RefinedLee(img) {

// img 必须使用自然单位，即不是 dB！设置 3x3 内核

// 转换为自然单位！

var myimg = toNatural(img);

var weights3 = ee.List.repeat(ee.List.repeat(1,3),3);

var kernel3 = ee.Kernel.fixed(3,3, weights3, 1, 1, false);

var mean3 = myimg.reduceNeighborhood(ee.Reducer.mean(), kernel3);

var variance3 = myimg.reduceNeighborhood(ee.Reducer.variance(), kernel3);

// 使用 7x7 窗户内的 3x3 窗户样本来确定梯度和方向

var sample_weights = ee.List([[0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0], [0,1,0,1,0,1,0], [0,0,0,0,0,0,0], [0,1,0,1,0,1,0],[0,0,0,0,0,0,0]]);

var sample_kernel = ee.Kernel.fixed(7,7, sample_weights, 3,3, false);

// 计算采样窗口的平均值和方差，并存储为 9 个波段

var sample_mean = mean3.neighborhoodToBands(sample_kernel);

var sample_var = variance3.neighborhoodToBands(sample_kernel);

// 确定采样窗口的 4 个梯度

var gradients = sample_mean.select(1).subtract(sample_mean.select(7)).abs();

gradients = gradients.addBands(sample_mean.select(6).subtract(sample_mean.select(2)).abs());

gradients = gradients.addBands(sample_mean.select(3).subtract(sample_mean.select(5)).abs());

gradients = gradients.addBands(sample_mean.select(0).subtract(sample_mean.select(8)).abs());

// 并找出梯度带中的最大梯度

var max_gradient = gradients.reduce(ee.Reducer.max());

// 为最大梯度的带状像素创建掩码

var gradmask = gradients.eq(max_gradient);

// 重复的梯度蒙板带：每个梯度代表 2 个方向

gradmask = gradmask.addBands(gradmask);

// 确定 8 个方向

var directions = sample_mean.select(1).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(7))).multiply(1);

directions = directions.addBands(sample_mean.select(6).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(2))).multiply(2));

directions = directions.addBands(sample_mean.select(3).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(5))).multiply(3));

directions = directions.addBands(sample_mean.select(0).subtract(sample_mean.select(4)).gt(sample_mean.select(4).subtract(sample_mean.select(8))).multiply(4));

// The next 4 are the not() of the previous 4

directions = directions.addBands(directions.select(0).not().multiply(5));

directions = directions.addBands(directions.select(1).not().multiply(6));

directions = directions.addBands(directions.select(2).not().multiply(7));

directions = directions.addBands(directions.select(3).not().multiply(8));

// 屏蔽所有非 1-8 的值

directions = directions.updateMask(gradmask);

// 将堆栈 “折叠 ”为单个波段图像（由于屏蔽，每个像素在其方向波段中只有一个值（1-8），否则会被屏蔽）。

directions = directions.reduce(ee.Reducer.sum());

var sample_stats = sample_var.divide(sample_mean.multiply(sample_mean));

// 计算本地噪声方差

var sigmaV = sample_stats.toArray().arraySort().arraySlice(0,0,5).arrayReduce(ee.Reducer.mean(), [0]);

// 为定向统计设置 7*7 内核

var rect_weights = ee.List.repeat(ee.List.repeat(0,7),3).cat(ee.List.repeat(ee.List.repeat(1,7),4));

var diag_weights = ee.List([[1,0,0,0,0,0,0], [1,1,0,0,0,0,0], [1,1,1,0,0,0,0],

[1,1,1,1,0,0,0], [1,1,1,1,1,0,0], [1,1,1,1,1,1,0], [1,1,1,1,1,1,1]]);

var rect_kernel = ee.Kernel.fixed(7,7, rect_weights, 3, 3, false);

var diag_kernel = ee.Kernel.fixed(7,7, diag_weights, 3, 3, false);

// 使用原始内核创建均值和方差堆栈。屏蔽相关方向。

var dir_mean = myimg.reduceNeighborhood(ee.Reducer.mean(), rect_kernel).updateMask(directions.eq(1));

var dir_var = myimg.reduceNeighborhood(ee.Reducer.variance(), rect_kernel).updateMask(directions.eq(1));

dir_mean = dir_mean.addBands(myimg.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)));

dir_var = dir_var.addBands(myimg.reduceNeighborhood(ee.Reducer.variance(), diag_kernel).updateMask(directions.eq(2)));

// 并为旋转的内核添加波段

for (var i=1; i<4; i++) {

dir_mean = dir_mean.addBands(myimg.reduceNeighborhood(ee.Reducer.mean(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)));

dir_var = dir_var.addBands(myimg.reduceNeighborhood(ee.Reducer.variance(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)));

dir_mean = dir_mean.addBands(myimg.reduceNeighborhood(ee.Reducer.mean(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)));

dir_var = dir_var.addBands(myimg.reduceNeighborhood(ee.Reducer.variance(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)));

}

// 将堆栈 “折叠 ”为单个波段图像（由于屏蔽，每个像素在其方向波段中只有一个值，否则就会被屏蔽）。

dir_mean = dir_mean.reduce(ee.Reducer.sum());

dir_var = dir_var.reduce(ee.Reducer.sum());

// A 最后生成过滤值

var varX = dir_var.subtract(dir_mean.multiply(dir_mean).multiply(sigmaV)).divide(sigmaV.add(1.0));

var b = varX.divide(dir_var);

var result = dir_mean.add(b.multiply(myimg.subtract(dir_mean)));

//return(result);

return(img.addBands(ee.Image(toDB(result.arrayGet(0))).rename("filter")));

}

var collection = S1.map(RefinedLee);

var col = ee.ImageCollection(collection.select("filter"));

print(ui.Chart.image.series(col, roi, ee.Reducer.mean(), 20).setOptions({

title: 'TimeSeries analysis',

lineWidth: 1,

pointSize: 3 }));