2023-02-05

1 Lee滤波去噪

  • Lee滤波是选择一定长度的窗口作为局部区域,利用图像局部统计特性进行SAR图像斑点滤波的典型方法,其是基于完全发育的斑点噪声模型(图像上斑点处于均匀区域,且斑点噪声与图像信号不相关)。

2 Refined Lee滤波去噪

  • Refined Lee滤波器是在Lee滤波器的基础上,改进的滤波算法,与Lee滤波器最主要的区别是采用了非正方形的滤波窗口。
  • 由于Lee滤波器存在缺陷,即对靠近边缘或点目标的同质区域像素滤波不够充分,之后,Lee又提出了一种基于边缘检测的自适应滤波算法,通过重新定义中心像素的邻域来提高估计的准确性。通常使用7×7的滑动窗口。
  • Refined Lee滤波器的算法实现主要包括两步:模板窗口选择和局部统计滤波
    • 模板窗口选择:将7×7的滑窗分为九个子区间,区间之间有重叠,每个子区间大小为3×3。
    • 局部统计滤波:在边缘方向窗口中,利用滤波器进行滤波。计算各子窗的均值。用这个均值构造一个3×3的矩阵M,来估计局域窗中边缘的方向,将3×3梯度模板应用到均值矩阵,梯度绝对值最大的方向被认为是边缘的方向。这里只需要用水平、垂直、45度和135度四个方向的梯度模板,相反方向互为相反数。

3 具体实现步骤与代码

  • 导入Python库,设置网络代理
import ee
import geemap
import os
geemap.set_proxy(port=19180)
Map = geemap.Map(center=(30, 120), zoom=4)
  • 设置研究区域
region =  ee.Geometry.Polygon(
        [[[98.95782793608701, 37.37657466863545],
          [98.95782793608701, 37.083535663646416],
          [99.50165117827451, 37.083535663646416],
          [99.50165117827451, 37.37657466863545]]])
# 定义从分贝转换的函数
def toNatural(img):
  return ee.Image(10.0).pow(img.divide(10.0))
# 定义转换成分贝的函数
def toDB(img):
  return ee.Image(img).log10().multiply(10.0)
# 定义RefinedLee主函数
def RefinedLee(img):
# 设置3x3内核
  img2 = img
  weights3 = ee.List.repeat(ee.List.repeat(1,3),3)
  kernel3 = ee.Kernel.fixed(3,3, weights3, 1, 1, False)
  mean3 = img2.reduceNeighborhood(ee.Reducer.mean(), kernel3)
  variance3 = img2.reduceNeighborhood(ee.Reducer.variance(), kernel3)
# 使用7x7窗口内的3x3窗口的样本来确定梯度和方向
  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]])
  sample_kernel = ee.Kernel.fixed(7,7, sample_weights, 3,3, False)
# 计算采样窗口的均值和方差,存储为9个波段
  sample_mean = mean3.neighborhoodToBands(sample_kernel); 
  sample_var = variance3.neighborhoodToBands(sample_kernel)
# 确定采样窗口的4个梯度
  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())
# 找到梯度带中的最大梯度
  max_gradient = gradients.reduce(ee.Reducer.max())
# 为最大梯度波段像素创建一个掩码
  gradmask = gradients.eq(max_gradient)
# 重复的渐变带: 每个渐变代表2个方向
  gradmask = gradmask.addBands(gradmask)
# 定义八个方向
  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))
# 接下来的4个是前面4个中的相反—Not()
  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()) 
#pal = ['ffffff','ff0000','ffff00', '00ff00', '00ffff', '0000ff', 'ff00ff', '000000']
#Map.addLayer(directions.reduce(ee.Reducer.sum()), {min:1, max:8, palette: pal}, 'Directions', False)
  sample_stats = sample_var.divide(sample_mean.multiply(sample_mean))
# 计算局部噪声方差
  sigmaV = sample_stats.toArray().arraySort().arraySlice(0,0,5).arrayReduce(ee.Reducer.mean(), [0])
# 设置7 * 7内核进行方向统计
  rect_weights = ee.List.repeat(ee.List.repeat(0,7),3).cat(ee.List.repeat(ee.List.repeat(1,7),4))
  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]])
  rect_kernel = ee.Kernel.fixed(7,7, rect_weights, 3, 3, False)
  diag_kernel = ee.Kernel.fixed(7,7, diag_weights, 3, 3, False)
# 使用原始内核创建均值和方差堆栈,用相关方向掩膜。
  dir_mean = img2.reduceNeighborhood(ee.Reducer.mean(), rect_kernel).updateMask(directions.eq(1))
  dir_var = img2.reduceNeighborhood(ee.Reducer.variance(), rect_kernel).updateMask(directions.eq(1))
  dir_mean = dir_mean.addBands(img2.reduceNeighborhood(ee.Reducer.mean(), diag_kernel).updateMask(directions.eq(2)))
  dir_var = dir_var.addBands(img2.reduceNeighborhood(ee.Reducer.variance(), diag_kernel).updateMask(directions.eq(2)))
# 然后加上旋转的内核
  for i in range(1,4):
    dir_mean = dir_mean.addBands(img2.reduceNeighborhood(ee.Reducer.mean(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)))
    dir_var = dir_var.addBands(img2.reduceNeighborhood(ee.Reducer.variance(), rect_kernel.rotate(i)).updateMask(directions.eq(2*i+1)))
    dir_mean = dir_mean.addBands(img2.reduceNeighborhood(ee.Reducer.mean(), diag_kernel.rotate(i)).updateMask(directions.eq(2*i+2)))
    dir_var = dir_var.addBands(img2.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())
# 最终生成过滤值
  varX = dir_var.subtract(dir_mean.multiply(dir_mean).multiply(sigmaV)).divide(sigmaV.add(1.0))
  b = varX.divide(dir_var)
  result = dir_mean.add(b.multiply(img2.subtract(dir_mean)))
  result = result.arrayGet(0).set('system:time_start', img.get('system:time_start'))
  return ee.Image(result)
  • S1雷达数据导入与筛选
S1 =  ee.ImageCollection('COPERNICUS/S1_GRD')\
  .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'))\
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))\
  .filter(ee.Filter.eq('instrumentMode', 'IW')).filterDate('2021-01-01', '2021-12-31')\
  .filterBounds(region).select(['VV','VH'])\
  .map(toNatural).map(RefinedLee).map(toDB)\
  .median().clip(region)
  • VV、VH波段与波段组合结果显示
sarParams_VV = {'bands': ['VV'],'min':-25, 'max':0}
sarParams_VH = {'bands': ['VH'],'min':-25, 'max':0}
sarParams_VVVH = {'bands': ['VV','VH'],'min':-25, 'max':0}
Map.addLayer(S1, sarParams_VV, 'VV')
Map.addLayer(S1, sarParams_VH, 'VH')
Map.addLayer(S1, sarParams_VVVH, 'VVVH')
Map
  • 0
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值