[GEE 6] Lab 6: Biến đổi ảnh và các chỉ số quang phổ

Lab 6: Biến đổi ảnh và các chỉ số quang phổ

Mục đích :  Mục đích của bài tập này là giúp bạn biết them về các chỉ số quang phổ, có thể được sử dụng để làm nổi bật đặc tính của các đối tượng mà chúng ta quan tâm trên ảnh viễn thám. Các bạn sẽ được giới thiệu các phương pháp để tạo các chỉ số về thực vật, nước, tuyết, bề mặt trống và bề mặt đã bị đốt cháy.

  1. Các chỉ số quang phổ

Các chỉ số quang phổ được xây dựng dựa trên phổ phản xạ khác nhau của các loại lớp phủ bề mặt. Các chỉ số được thiết kế để khai thác những sự khác biệt để làm nổi bật các loại đất cụ thể. Hãy tìm hiểu đồ thị về phản xạ phổ của các đối tượng bề mặt khác nhau trong đồ thị dưới đây.

Quan sát cho thấy các loại lớp phủ phản xạ khác nhau tại một hoặc nhiều bước sóng. Ví dụ đường cong phản xạ của thực vật (xanh) có hệ số phản xạ tương đối cao trong khoảng cận hồng ngoại (NIR), do năng lượng bức xạ tán xạ bởi các tế bào trong lá cây (Bowker et al. 1985). Bên cạnh đó. thực vật có phản xạ thấp trong dải màu đỏ, nơi năng lượng bức xạ được hấp thụ bởi chất diệp lục. Những quan sát này thúc đẩy việc xây dựng các chỉ số thực vật, ví dụ:

  1. NDVI.  

Chỉ số chuẩn hóa khác biệt về thực vật – The Normalized Difference Vegetation Index (NDVI) có lịch sử lâu đời trong viễn thám.  Công thức đặc trưng là

NDVI = (NIR – red) / (NIR + red)

Với kênh cận hồng ngoại và đỏ thể hiện phản xạ, tán xạ hoặc giá trị số cụ thể tại mỗi bước sóng cụ thể. Thực hiện chỉ số này trong GEE sẽ sử dụng phương pháp normalizedDifference(). Trước hết lấy một ảnh tại khu vực quan tâm bằng cách vẽ một điểm (Point), đặt tên điểm và nhập từ tập hợp ảnh Landsat8 TOA như là

Được gọi là landsat8 và sắp xếp tập hợp ảnh theo dữ liệu metadata về mây che phủ:

var image = ee.Image(landsat8

    .filterBounds(point)

    .filterDate(‘2015-06-01’, ‘2015-09-01’)

    .sort(‘CLOUD_COVER’)

    .first());

var trueColor = {bands: [‘B4’, ‘B3’, ‘B2’], min: 0, max: 0.3};

Map.addLayer(image, trueColor, ‘image’);

Công thức tính toán NDVI là dòng lệnh dưới đây:

var ndvi = image.normalizedDifference([‘B5’, ‘B4’]);

Hiển thị ảnh NDVI với bảng màu (có thể thay đổi)

var vegPalette = [‘red’, ‘blue’, ‘yellow’, ‘green’];

Map.addLayer(ndvi, {min: -1, max: 1, palette: vegPalette}, ‘NDVI’);

Sử dụng Inspector để kiểm tra các giá trị pixel trong các khu vực có thực vật và không có thực vật phủ.

  1. EVI.  

Chỉ số thực vật tăng cường (EVI) được thiết kế để giảm thiểu hiệu ứng bão hoà và các yếu tố nền của chỉ số NDVI (Huete et al. 2002). Vì nó không phải là một chỉ số khác biệt chuẩn, nó được tính toán theo biểu thức:

var evi = image.expression(

    ‘2.5 * ((NIR – RED) / (NIR + 6 * RED – 7.5 * BLUE + 1))’, {

      ‘NIR’: image.select(‘B5’),

      ‘RED’: image.select(‘B4’),

      ‘BLUE’: image.select(‘B2’)

});

Quan sát các kênh ảnh  được tham chiếu với sự giúp đỡ của một đối tượng sau khi được tính toán theo biểu thức image.expression().  

Hiển thị ảnh EVI:

Map.addLayer(evi, {min: -1, max: 1, palette: vegPalette}, ‘EVI’);

SO sánh kết quả của EVI và  NDVI. Bạn thấy có điểm khác biệt gi?

  1. NDWI.  

Chỉ số chuẩn hóa khác biệt về nước (NDWI) được phát triển bởi Gao (1996) là chỉ số về lượng thực vật và nước :

NDWI = (NIR – SWIR) / (NIR + SWIR)

Tính chỉ số NDWI trong GEE với công thức:

var ndwi = image.normalizedDifference([‘B5’, ‘B6’]);

Và hiển thị :

var waterPalette = [‘red’, ‘yellow’, ‘green’, ‘blue’];

Map.addLayer(ndwi, {min: -0.5, max: 1, palette: waterPalette}, ‘NDWI’);

Lưu ý rằng đây không phải là việc thực hiện chính xác NDWI, theo OLI spectral response, vì OLI không có kênh ảnh theo bước sóng đúng (1.26 𝛍m).

  1. NDWBI.  

Hai chỉ số NDWI khác nhau đã được phát minh một cách độc lập vào năm 1996. Để phân biệt, công thức NDWBI – chuẩn khác biệt về các bề  mặt nước là chỉ số được mô tả trong McFeeters (1996)

NDWBI = (green – NIR) / (green + NIR)

Giống như trên, thực hiện công thức NDWBI theo hàm normalizedDifference() và hiển thị kết quả:

var ndwi = image.normalizedDifference([‘B3’, ‘B5’]);

Map.addLayer(ndwi, {min: -1, max: 0.5, palette: waterPalette}, ‘NDWBI’);

So sánh kết quả NDWIvà NDWBI.  Bạn thấy có khác biệt gì ?

  1. NDBI.  

Chỉ số Chuẩn khác biệt về bề mặt trống Normalized Difference Bare Index (NDBI) được phát triển bởi  Zha et al. (2003)  để giúp đánh giá sự khác biệt ở khu vực đô thị:

NDBI = (SWIR – NIR) / (SWIR + NIR)

Lưu ý rằng chỉ số NDBI là chỉ số đảo của NDWI. Tính NDBI và hiển thị với bảng màu tùy chọn:

var ndbi = image.normalizedDifference([‘B6’, ‘B5’]);

var barePalette = waterPalette.slice().reverse();

Map.addLayer(ndbi, {min: -1, max: 0.5, palette: barePalette}, ‘NDBI’);

(Check this reference to demystify the palette reversal).

  1. BAI.  

Chỉ số cháy (BAI) được phát triển bởi Chuvieco et al. (2002) để hỗ trợ trong việc mô tả những khu vực bị cháy và đánh giá mức độ nghiêm trọng của đám cháy. Nó được dựa trên khoảng cách phản xạ phổ của than. Để kiểm tra chỉ số cháy, tải một hình ảnh từ đám cháy ở Sierra Nevadas năm 2013:

var burnImage = ee.Image(landsat8

    .filterBounds(ee.Geometry.Point(-120.083, 37.850))

    .filterDate(‘2013-08-17’, ‘2013-09-27’)

    .sort(‘CLOUD_COVER’)

    .first());

Map.addLayer(burnImage, trueColor, ‘burn image’);

Khi kiểm tra kỹ lưỡng theo hiển thị màu thực của bức ảnh này, bạn có thể chỉ ra các vị trí có đám cháy? Nếu không công thức BAI có thể giúp. Giống như EVI, chí số BAI được tính trên GEE như sau

:

var bai = burnImage.expression(

    ‘1.0 / ((0.1 – RED)**2 + (0.06 – NIR)**2)’, {

      ‘NIR’: burnImage.select(‘B5’),

      ‘RED’: burnImage.select(‘B4’),

});

Hiển thị kết quả. Khu vực bị cháy sẽ hiển thị ró ràng hơn trên ảnh kết quả BAI

var burnPalette = [‘green’, ‘blue’, ‘yellow’, ‘red’];

Map.addLayer(bai, {min: 0, max: 400, palette: burnPalette}, ‘BAI’);

  1. NBRT.  

Tỉ số chuẩn hóa nhiệt (NBRT) được phát triển dựa trên ý tưởng rằng bế mặt đất bị đốt cháy sẽ phản xạ thấp ở kênh cận hồng ngoại (ít thực vật), phản xạ cao ở kênh SWIR (tro), và độ sáng cao với nhiệt độ (Holden et al. 2005). Không giống như các chỉ số khác, một pixel có chỉ só NBRT thấp hơn có nghĩa là bị cháy nhiều hơn. Thực hiện NBRT với biểu thức:

var nbrt = burnImage.expression(

  ‘(NIR – 0.0001 * SWIR * Temp) / (NIR + 0.0001 * SWIR * Temp)’, {

    ‘NIR’: burnImage.select(‘B5’),

    ‘SWIR’: burnImage.select(‘B7’),

    ‘Temp’: burnImage.select(‘B11’)

});

Hiển thị kết quả và đảo ngược tý lệ :

Map.addLayer(nbrt, {min: 1, max: 0.9, palette: burnPalette}, ‘NBRT’);

Sự khác biệt trong chỉ số này, trước và sau khi cháy có thể được sử dụng để xác định các khu vực cháy nghiêm trọng (see van Wagtendonk et al. 2004).

  1. NDSI.  

Chỉ số chuẩn về tuyết (NDSI) được thiết kế để ước lượng số lượng pixel có tuyết phủ

 (Riggs et al. 1994):

NDSI = (green – SWIR) / (green + SWIR)

Trước hết, tìm một cảnh ảnh có tuyết để kiểm tra chí số:

var snowImage = ee.Image(landsat8

    .filterBounds(ee.Geometry.Point(-120.0421, 39.1002))

    .filterDate(‘2013-11-01’, ‘2014-05-01’)

    .sort(‘CLOUD_COVER’)

    .first());

Map.addLayer(snowImage, trueColor, ‘snow image’);

Tính và hiển thị kết quả NDSI trên GEE:

var ndsi = snowImage.normalizedDifference([‘B3’, ‘B6’]);

var snowPalette = [‘red’, ‘green’, ‘blue’, ‘white’];

Map.addLayer(ndsi, {min: -0.5, max: 0.5, palette: snowPalette}, ‘NDSI’);

  1. Chuyển đổi tuyến tính

Phép biến đổi tuyến tính là sự kết hợp tuyến tính của các giá trị pixel đầu vào. Đây có thể là kết quả của một loạt các chiến lược khác nhau, nhưng một chủ đề chung là pixel đang được coi là một dãy các giá trị của các kênh ảnh.

  1. Tasseled cap (TC).

Dựa trên những quan sát về đất nông nghiệp phản xạ trong khoảng phổ giữa cận hồng ngoại và đỏ, Kauth và Thomas (1976) đã phát minh ra cách luân chuyển theo form

p1 = RTp0

 Với p0 là vector px1 điểm ảnh gốc (một tập hợp các giá trị p được coi như một mảng các giá trị), p1 là pixel xoay và R là một cơ sở trực chuẩn của không gian mới (do đó RT là nghịch đảo của nó). Kauth và Thomas tìm thấy R bằng cách xác định trục đầu tiên của không gian chuyển đổi để song song với đường phản xạ của đất trong biểu đồ dưới đây, sau đó sử dụng phương pháp Gram-Schmidt để tìm các vectơ cơ sở khác.

Giả sử R là một biến. một chiều để tiến hành xoay trục trong GEE là với các mảng (array) này. Cụ thể, tạo một mảng của hệ số TC

var coefficients = ee.Array([

  [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],

  [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],

  [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],

  [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],

  [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],

  [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]

]);

Các hệ số ở trên là hệ số của vệ tinh ảnh Landsat TM sensor, tìm một ảnh ít mây. Trước hết tìm ‘landsat 5 toa’, sau đó Import ‘USGS Landsat 5 TOA Reflectance (Orthorectified)’.  Đặt tên ảnh mới nhập landsat5, sau đó lọc lựa chọn ảnh và sắp xếp theo mức độ mây che phủ:

var tcImage = ee.Image(landsat5

    .filterBounds(point)

    .filterDate(‘2008-06-01’, ‘2008-09-01’)

    .sort(‘CLOUD_COVER’)

    .first());

Để thực hiện phép nhân ma trận, lầu tiên chuyển đổi ảnh ban đầu từ một ảnh nhiều kênh phổ trở thành ảnh 1 kênh với các dãy giá trị, trong đó mỗi pixel chứa một mảng giá trị số:

var bands = [‘B1’, ‘B2’, ‘B3’, ‘B4’, ‘B5’, ‘B7’];

// Make an Array Image, with a 1-D Array per pixel.

var arrayImage1D = tcImage.select(bands).toArray();

// Make an Array Image with a 2-D Array per pixel, 6×1.

var arrayImage2D = arrayImage1D.toArray(1);

Thực hiện nhân ma trận, và chuyển đổi trở về ảnh đa kênh phổ:

var componentsImage = ee.Image(coefficients)

  .matrixMultiply(arrayImage2D)

  // Get rid of the extra dimensions.

  .arrayProject([0])

  // Get a multi-band image with TC-named bands.

  .arrayFlatten(

    [[‘brightness’, ‘greenness’, ‘wetness’, ‘fourth’, ‘fifth’, ‘sixth’]]);

Cuối cùng, hiển thị kết quả:

var vizParams = {

  bands: [‘brightness’, ‘greenness’, ‘wetness’],

  min: -0.1, max: [0.5, 0.1, 0.1]

};

Map.addLayer(componentsImage, vizParams, ‘TC components’);