CẢNH BÁO: Các quy trình được mô tả trong hướng dẫn này để điều chỉnh dữ liệu Landsat ETM+ và OLI đã lỗi thời và không được khuyến nghị hoặc không cần thiết khi làm việc với dữ liệu hệ số phản xạ bề mặt Landsat Collection 2. Để biết thêm thông tin, hãy xem mục Câu hỏi thường gặp về việc chuẩn hoá Landsat. Dữ liệu của Tập hợp 1 không còn được dùng nữa và mã trong hướng dẫn này sẽ không còn hoạt động hoặc được cập nhật.
Hướng dẫn này liên quan đến việc điều chỉnh độ phản xạ bề mặt của Landsat ETM+ cho phù hợp với độ phản xạ bề mặt của Landsat OLI. Thư viện này cung cấp:
- Hàm biến đổi phổ
- Các hàm để tạo dữ liệu sẵn sàng phân tích
- Ví dụ về hình ảnh trực quan hoá chuỗi thời gian
Đây là hướng dẫn toàn diện về cách điều chỉnh và trực quan hoá chuỗi thời gian khu vực gồm dữ liệu Landsat trong hơn 35 năm. Bạn có thể áp dụng ngay hướng dẫn này cho(các) khu vực mà bạn quan tâm.
Xin lưu ý rằng các hệ số để chuyển đổi ETM+ sang OLI cũng áp dụng cho TM. Do đó, trong hướng dẫn này, ETM+ được dùng thay cho TM.
Giới thiệu về Landsat
Landsat là một chương trình chụp ảnh vệ tinh đã thu thập hình ảnh Trái đất có độ phân giải trung bình từ năm 1972. Là chương trình quan sát Trái Đất từ không gian dài nhất, Landsat cung cấp một bản ghi thời gian có giá trị để xác định các xu hướng không gian và thời gian trong sự thay đổi của cảnh quan. Dữ liệu của công cụ Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) và Operational Land Imager (OLI) được dùng trong hướng dẫn này. Chúng có mối quan hệ chặt chẽ và tương đối dễ dàng hợp nhất thành một chuỗi thời gian nhất quán, tạo ra một bản ghi liên tục từ năm 1984 đến nay, với tần suất 16 ngày/cảm biến và độ phân giải không gian là 30 mét. Thiết bị Máy quét đa phổ mở rộng bản ghi Landsat từ năm 1972, nhưng không được dùng trong hướng dẫn này. Dữ liệu của cảm biến này khá khác biệt, gây khó khăn cho việc tích hợp với các cảm biến sau này. Hãy xem Savage và cộng sự (2018) và Vogeler và cộng sự (2018) để biết ví dụ về việc điều chỉnh trên tất cả các cảm biến.
Lý do cần hài hoà
Roy và cộng sự (2016) cho thấy có những điểm khác biệt nhỏ nhưng có khả năng đáng kể giữa đặc điểm quang phổ của Landsat ETM+ và OLI, tuỳ thuộc vào ứng dụng. Lý do bạn nên điều chỉnh tập dữ liệu bao gồm: tạo chuỗi thời gian dài trải dài trên Landsat TM, ETM+ và OLI, tạo các thành phần gần ngày trong năm để giảm ảnh hưởng của các quan sát bị thiếu do khoảng trống SLC-off của ETM+ và việc che phủ đám mây/bóng, hoặc tăng tần suất quan sát trong chuỗi thời gian. Vui lòng xem bản thảo được liên kết ở trên để biết thêm thông tin.
Hướng dẫn
Hàm
Sau đây là một loạt các hàm cần thiết để điều chỉnh ETM+ thành OLI và tạo dữ liệu sẵn sàng để phân tích. Dữ liệu này sẽ được dùng trong phần ứng dụng của hướng dẫn này để trực quan hoá lịch sử quang phổ-thời gian của một pixel.
Harmonization
Việc điều chỉnh được thực hiện thông qua phép biến đổi tuyến tính của không gian quang phổ ETM+ thành không gian quang phổ OLI theo các hệ số được trình bày trong Bảng 2 của Roy và cộng sự (2016) về các hệ số hồi quy OLS.
Các hệ số tương ứng với băng tần được xác định trong từ điển sau đây với các hằng số hình ảnh độ dốc (slopes) và đường cắt (itcps). Xin lưu ý rằng các giá trị giao điểm trên trục y được nhân với 10.000 để khớp với tỷ lệ của dữ liệu độ phản xạ bề mặt Landsat của USGS.
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
Tên dải ETM+ và OLI cho cùng một cửa sổ phản xạ quang phổ không bằng nhau và cần được chuẩn hoá. Các hàm sau đây chỉ chọn các dải phản xạ và dải pixel_qa từ mỗi tập dữ liệu, đồng thời đổi tên chúng theo dải bước sóng mà chúng đại diện.
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
Cuối cùng, hãy xác định hàm biến đổi, hàm này sẽ áp dụng mô hình tuyến tính cho dữ liệu ETM+, truyền kiểu dữ liệu dưới dạng Int16 để nhất quán với OLI và đính kèm lại dải pixel_qa để sử dụng sau này trong việc che phủ đám mây và bóng.
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
Che mây và bóng
Dữ liệu sẵn sàng để phân tích phải được che phủ bằng mây và bóng mây. Hàm sau đây sử dụng CFmask (Zhu và cộng sự, 2015)
Băng tần pixel_qa được đưa vào mỗi hình ảnh độ phản xạ bề mặt Landsat USGS để đặt các pixel được xác định là đám mây và bóng mây thành giá trị rỗng.
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
Tính toán chỉ số quang phổ
Ứng dụng sắp ra mắt này sử dụng chỉ số quang phổ tỷ lệ cháy được chuẩn hoá (NBR) để biểu thị lịch sử quang phổ của một điểm ảnh rừng bị ảnh hưởng bởi cháy rừng. NBR được sử dụng vì Cohen và cộng sự (2018) nhận thấy rằng trong số 13 chỉ số/băng tần quang phổ, NBR có tỷ lệ tín hiệu trên nhiễu lớn nhất liên quan đến tình trạng xáo trộn rừng (tín hiệu) trên khắp Hoa Kỳ. Giá trị này được tính bằng hàm sau.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
Trong ứng dụng của riêng mình, bạn có thể chọn sử dụng một chỉ số quang phổ khác. Đây là nơi bạn có thể thay đổi hoặc thêm một chỉ mục khác.
Kết hợp các hàm
Xác định một hàm bao bọc cho từng tập dữ liệu để hợp nhất tất cả các hàm ở trên cho thuận tiện khi áp dụng chúng vào các tập hợp hình ảnh tương ứng.
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
Trong ứng dụng của mình, bạn có thể quyết định đưa vào hoặc loại trừ các hàm. Thay đổi các hàm này nếu cần.
Ví dụ về chuỗi thời gian
Các hàm bao bọc prepOli và prepEtm ở trên có thể được ánh xạ trên các tập hợp hệ số phản xạ bề mặt Landsat để tạo dữ liệu sẵn sàng phân tích trên nhiều cảm biến nhằm trực quan hoá trình tự quang phổ của một pixel hoặc vùng pixel. Trong ví dụ này, bạn sẽ tạo một chuỗi thời gian hơn 35 năm và hiển thị nhật ký quang phổ cho một pixel duy nhất. Pixel cụ thể này liên quan đến lịch sử gần đây của một mảng rừng cây lá kim trưởng thành ở Tây Bắc Thái Bình Dương (Hình 1), nơi đã trải qua một số biến động vào những năm 1980 và 1990, cũng như một vụ cháy lớn vào năm 2011.
Hình 1. Vị trí và ký tự trang web cho ví dụ về khu vực mà bạn quan tâm. Rừng lá kim trưởng thành ở Tây Bắc Thái Bình Dương trên sườn phía bắc của Núi Hood, tiểu bang Oregon, Hoa Kỳ. Hình ảnh do: Google Earth, Cục Lâm nghiệp Hoa Kỳ, Landsat và Copernicus cung cấp.
Xác định khu vực quan tâm (AOI)
Kết quả của ứng dụng này là một chuỗi thời gian gồm các dữ liệu quan sát Landsat cho một pixel. Đối tượng ee.Geometry.Point được dùng để xác định vị trí của pixel.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
Trong ứng dụng của mình, bạn có thể chọn một pixel khác. Bạn có thể thực hiện việc này bằng cách thay thế toạ độ kinh độ và vĩ độ ở trên. Ngoài ra, bạn có thể tóm tắt nhật ký quang phổ của một nhóm pixel bằng cách sử dụng các định nghĩa đối tượng ee.Geometry khác, chẳng hạn như ee.Geometry.Polygon() và ee.Geometry.Rectangle().
Vui lòng xem phần Hình học trong Hướng dẫn cho nhà phát triển để biết thêm thông tin.
Truy xuất các bộ sưu tập dữ liệu từ cảm biến Landsat
Nhận bộ sưu tập độ phản xạ bề mặt Landsat USGS cho OLI, ETM+ và TM.
Hãy truy cập vào các đường liên kết để tìm hiểu thêm về từng tập dữ liệu.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
Xác định bộ lọc bộ sưu tập hình ảnh
Xác định bộ lọc để giới hạn các tập hợp hình ảnh theo ranh giới không gian của khu vực quan tâm, mùa quang hợp cao điểm và chất lượng.
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
Bạn có thể thay đổi các tiêu chí lọc này theo ý muốn. Xác định các thuộc tính khác để lọc theo trong thẻ "Thuộc tính hình ảnh" của các trang mô tả dữ liệu được liên kết ở trên.
Chuẩn bị bộ sưu tập
Hợp nhất các tập hợp, áp dụng bộ lọc và ánh xạ hàm prepImg trên tất cả các hình ảnh. Kết quả của đoạn mã sau đây là một ee.ImageCollection duy nhất chứa hình ảnh từ các bộ sưu tập cảm biến OLI, ETM+ và TM đáp ứng tiêu chí lọc và được xử lý thành NBR sẵn sàng phân tích.
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
Tạo biểu đồ chuỗi thời gian hiển thị tất cả các điểm quan sát
Bạn có thể đánh giá định tính sự hài hoà giữa các cảm biến bằng cách vẽ tất cả các kết quả quan sát dưới dạng biểu đồ phân tán, trong đó cảm biến được phân biệt bằng màu sắc. Hàm ui.Chart.feature.groups của Earth Engine cung cấp tiện ích này, nhưng trước tiên, bạn phải thêm giá trị NBR quan sát được vào từng hình ảnh dưới dạng một thuộc tính. Ánh xạ một hàm giảm khu vực trên tập hợp hình ảnh để tính toán giá trị trung bình của tất cả các pixel giao nhau với AOI và đặt kết quả làm thuộc tính hình ảnh.
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
Trong ứng dụng của riêng bạn, nếu AOI là một khu vực rộng lớn, bạn có thể cân nhắc tăng scale và chỉ định các đối số bestEffort, maxPixels, tileScale cho hàm reduceRegion để đảm bảo rằng thao tác không vượt quá giới hạn về số lượng pixel tối đa, bộ nhớ hoặc thời gian chờ. Bạn cũng có thể thay thế hàm giảm trung vị bằng số liệu thống kê mà bạn muốn. Hãy xem phần Số liệu thống kê về một vùng hình ảnh trong Hướng dẫn cho nhà phát triển để biết thêm thông tin.
Hàm ui.Chart.feature.groups hiện có thể chấp nhận bộ sưu tập này.
Hàm này mong đợi một tập hợp và tên của các thuộc tính đối tượng tập hợp để ánh xạ đến biểu đồ. Ở đây, thuộc tính "system:time_start" đóng vai trò là biến trục x và "NBR" đóng vai trò là biến trục y. Xin lưu ý rằng thuộc tính NBR được đặt tên bằng hàm reduceRegion ở trên dựa trên tên dải tần trong hình ảnh đang được giảm. Nếu bạn sử dụng một dải tần khác (không phải "NBR"), hãy thay đổi đối số tên thuộc tính trục y cho phù hợp. Cuối cùng, hãy đặt biến nhóm (loạt) thành "SATELLITE" (VỆ TINH), trong đó các màu riêng biệt được chỉ định.
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
Một biểu đồ tương tự như Hình 2 sẽ xuất hiện trong bảng điều khiển sau một khoảng thời gian xử lý. Lưu ý một số điều sau:
- Mỗi năm có nhiều lần quan sát.
- Chuỗi thời gian này bao gồm 3 cảm biến khác nhau.
- Không có độ lệch cảm biến nào có thể nhận biết được.
- Tần suất quan sát thay đổi hằng năm.
- Có một số biến động trong năm, nhất quán qua các năm và tương đối chặt chẽ. Thử nghiệm này với NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) cho thấy sự biến đổi phản hồi trong và giữa các năm lớn hơn nhiều.
- NBR vẫn ở mức cao trong phần lớn chuỗi thời gian với một số biến động nhỏ.
- NBR giảm đáng kể và nhanh chóng do cháy rừng, thể hiện rõ trong
- Tuy nhiên, xin lưu ý rằng vụ cháy (Hồ Dollar) xảy ra vào tháng 9 năm 2011. Lỗi về năm phát hiện là do phạm vi ngày tổng hợp hằng năm là từ tháng 7 đến tháng 8; những thay đổi xảy ra sau khoảng thời gian này sẽ không được ghi nhận cho đến khi có được quan sát khác không phải là giá trị rỗng, có thể là vào năm sau hoặc muộn hơn.
- Quá trình phục hồi thảm thực vật (tăng phản hồi NBR) bắt đầu sau 2 năm kể từ khi phát hiện thấy sự sụt giảm lớn về NBR.
Hình 2. Biểu đồ chuỗi thời gian phản hồi quang phổ cho một pixel duy nhất với dữ liệu đại diện từ Landsat TM, ETM+ và OLI. Hình ảnh TM và ETM+ được điều chỉnh thành OLI bằng phép biến đổi tuyến tính. Hình ảnh được chụp từ tháng 7 đến tháng 8 và được lọc để đảm bảo chất lượng cao.
Tạo biểu đồ chuỗi thời gian cho thấy giá trị trung vị hằng năm
Để đơn giản hoá và loại bỏ nhiễu khỏi chuỗi thời gian, bạn có thể áp dụng phương pháp giảm số lượng quan sát trong năm. Ở đây, giá trị trung vị được sử dụng.
Bước đầu tiên là nhóm ảnh theo năm. Thêm một thuộc tính "year" mới vào mỗi hình ảnh bằng cách ánh xạ qua chế độ cài đặt bộ sưu tập "year" từ ee.Image.Date của mỗi hình ảnh.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
Bạn có thể dùng thuộc tính "year" (năm) mới để nhóm các hình ảnh trong cùng một năm bằng một thao tác kết hợp. Thao tác kết hợp giữa một tập hợp năm riêng biệt của hình ảnh và tập hợp hoàn chỉnh sẽ cung cấp danh sách cùng năm, từ đó có thể thực hiện các thao tác giảm trung bình. Tạo một tập hợp con từ bộ sưu tập đầy đủ thành một tập hợp gồm các đại diện riêng biệt theo năm.
var distinctYearCol = col.distinct('year');
Xác định một bộ lọc và một phép kết hợp để xác định những năm trùng khớp giữa tập hợp năm riêng biệt (distinctYearCol) và tập hợp đầy đủ (col). Các kết quả trùng khớp sẽ được lưu vào một thuộc tính có tên là "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Áp dụng thao tác kết hợp và chuyển đổi FeatureCollection kết quả thành ImageCollection.
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
Thay đổi ee.Reducer ở trên nếu bạn muốn dùng một số liệu thống kê khác ngoài số trung vị để biểu thị bộ sưu tập hình ảnh trong một năm nhất định.
Cuối cùng, hãy tạo một biểu đồ hiển thị các giá trị NBR trung bình từ các nhóm quan sát mùa hè trong năm. Thay đổi số liệu thống kê về mức giảm của khu vực theo ý muốn.
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
Một biểu đồ tương tự như Hình 3 sẽ xuất hiện trong bảng điều khiển sau một khoảng thời gian xử lý.
Hình 3. Biểu đồ chuỗi thời gian phản hồi quang phổ trung bình vào mùa hè cho một pixel duy nhất có dữ liệu đại diện từ Landsat TM, ETM+ và OLI. Hình ảnh TM và ETM+ được điều chỉnh thành OLI bằng một phép biến đổi tuyến tính. Hình ảnh được chụp từ tháng 7 đến tháng 8 và được lọc để đảm bảo chất lượng cao.
Thông tin khác
Các hàm biến đổi thay thế
Bảng 2 của Roy và cộng sự (2016) cung cấp các hệ số hồi quy OLS và RMA để chuyển đổi độ phản xạ bề mặt ETM+ thành độ phản xạ bề mặt OLI và ngược lại. Hướng dẫn trên chỉ minh hoạ quá trình chuyển đổi ETM+ sang OLI bằng hồi quy OLS. Bạn có thể tìm thấy các hàm cho tất cả lựa chọn dịch trong tập lệnh Trình chỉnh sửa mã hoặc tập lệnh nguồn GitHub.
Vụ cháy ở hồ Dollar
Vụ cháy được đề cập trong hướng dẫn này là vụ cháy Dollar Lake, xảy ra vào giữa tháng 9 năm 2011 ở sườn phía bắc của núi Hood, Oregon, Hoa Kỳ. Hãy truy cập vào bài đăng này trên blog và xem Varner và cộng sự (2015) để biết thêm thông tin về vấn đề này.