Giới thiệu dự án

Tốc độ đô thị hóa nhanh chóng tại các vùng kinh tế trọng điểm phía Bắc, đặc biệt là tỉnh Quảng Ninh, đặt ra bài toán cấp bách trong công tác quản lý tài nguyên đất đai. Thị xã Đông Triều – cửa ngõ phía Tây tỉnh Quảng Ninh – đang trải qua quá trình chuyển dịch cơ cấu kinh tế mạnh mẽ từ nông nghiệp sang công nghiệp, dịch vụ và phát triển đô thị. Quá trình này thúc đẩy sự gia tăng diện tích bề mặt không thấm (đất xây dựng) cùng với sự xuất hiện của các diện tích đất trống do giải phóng mặt bằng, chuyển mục đích sử dụng đất hoặc các dự án quy hoạch dở dang.

Phương pháp đo đạc địa chính truyền thống đòi hỏi chi phí nhân lực lớn, chu kỳ cập nhật kéo dài từ 3–5 năm và gặp nhiều trở ngại khi theo dõi biến động diện tích quy mô lớn. Ứng dụng công nghệ viễn thám đa phổ (Multispectral Remote Sensing) kết hợp hệ thống thông tin địa lý (GIS) cung cấp giải pháp giám sát không gian diện rộng với tần suất chụp lặp định kỳ, tối ưu hóa ngân sách và nâng cao độ chính xác trong công tác quản lý Nhà nước về đất đai.

+-----------------------------------------------------------------------------------+
|                              VẤN ĐỀ VÀ THÁCH THỨC                                  |
| - Đô thị hóa nhanh làm xáo trộn lớp phủ bề mặt (đất xây dựng vs đất trống)        |
| - Hiện tượng pixel hỗn hợp (mixed pixels) gây nhầm lẫn phổ giữa đất trống & bê tông|
| - Phương pháp khảo sát thực địa tốn kém chi phí, dữ liệu cập nhật chậm trễ         |
+---------------------------------------------------------+-------------------------+
                                                          |
                                                          v
+---------------------------------------------------------+-------------------------+
|                                GIẢI PHÁP ĐỀ XUẤT                                  |
| - Khai thác ảnh vệ tinh Landsat 8/9 (OLI/TIRS) độ phân giải 30m                   |
| - Áp dụng phương pháp phân tích chỉ số phổ: NDBI, UI, EBBI                         |
| - Kết hợp phân loại có kiểm định (Supervised Classification) và lọc nhiễu 3x3      |
| - Đánh giá định lượng qua ma trận sai số (Confusion Matrix) và hệ số Kappa         |
+-----------------------------------------------------------------------------------+

Mục tiêu nghiên cứu

  1. Thu thập và tiền xử lý dữ liệu: Chuẩn hóa dữ liệu ảnh viễn thám Landsat 8/9 OLI/TIRS chụp ngày 08/04/2022 trên địa bàn thị xã Đông Triều thông qua hiệu chỉnh bức xạ, hiệu chỉnh khí quyển và hình học.
  2. Chiết tách thông tin phổ: Xây dựng và đánh giá hiệu quả của các chỉ số phổ NDBI (Normalized Difference Built-Up Index), UI (Index-Based Built-Up Index), EBBI (Enhanced Built-Up and Bareness Index) so sánh với phương pháp phân loại có kiểm định (Maximum Likelihood).
  3. Thành lập bản đồ chuyên đề: Xây dựng bản đồ phân bố không gian đất xây dựng và đất trống tỷ lệ 1:50.000, cung cấp số liệu thống kê diện tích phục vụ công tác quy hoạch đô thị Đông Triều.

Phạm vi và giới hạn nghiên cứu

  • Phạm vi không gian: Toàn bộ ranh giới hành chính thị xã Đông Triều, tỉnh Quảng Ninh (bao gồm các phường trọng điểm như Mạo Khê, Đông Triều, Đức Chính, Hồng Phong và các xã lân cận).
  • Phạm vi dữ liệu: Cảnh ảnh vệ tinh Landsat 8 OLI/TIRS (Path 126, Row 045) chụp ngày 08/04/2022 với tỷ lệ mây dưới 5%.
  • Giới hạn kỹ thuật: Độ phân giải không gian của các kênh phổ phản xạ là 30m và kênh nhiệt là 100m (tái lấy mẫu về 30m), phù hợp cho bản đồ chuyên đề cấp thị xã nhưng cần kết hợp dữ liệu độ phân giải cao hơn đối với các công trình đơn lẻ.

Phân tích và thiết kế giải pháp

Phân tích hiện trạng

Chiết tách bề mặt nhân tạo và đất trống đối mặt với thách thức lớn về sự tương đồng trong phản xạ phổ (Spectral Signature Overlap). Bề mặt bê tông, mái ngói, đường nhựa và nền đất khô đều có độ phản xạ tăng mạnh ở dải sóng hồng ngoại ngắn (SWIR) và giảm ở dải cận hồng ngoại (NIR).

Phương pháp Ưu điểm Nhược điểm Độ chính xác trung bình
Đo vẽ thực địa & Số hóa thủ công Độ chính xác tuyệt đối ở quy mô vi mô, chi tiết từng thửa đất. Chi phí nhân công rất cao, chu kỳ cập nhật dài (hàng năm), tốn thời gian. > 95%
Phân loại có kiểm định (MLC/SVM) Phân tách được nhiều lớp phủ phức tạp, tận dụng toàn bộ các kênh phổ. Phụ thuộc chặt chẽ vào chất lượng mẫu huấn luyện; dễ nhầm lẫn pixel hỗn hợp. 78% – 84%
Phân tích chỉ số phổ (NDBI / UI / EBBI) Tốc độ xử lý nhanh, thuật toán tinh gọn, dễ tự động hóa trên diện rộng. Cần xác định ngưỡng phân tách (Thresholding) tối ưu cho từng vùng địa lý. 85% – 92%

Ma trận ưu tiên yêu cầu hệ thống (MoSCoW)

  • Must Have (Bắt buộc):
    • Khả năng xử lý tệp raster đa kênh Landsat OLI/TIRS chuẩn GeoTIFF.
    • Module tính toán tự động 3 chỉ số chuyên biệt: NDBI, UI, EBBI.
    • Phân tách độc lập 3 nhóm: Đất xây dựng, Đất trống, Đối tượng khác (Thực vật, Mặt nước).
    • Xuất ma trận sai số (Confusion Matrix), Overall Accuracy (OA) $\ge 85%$ và chỉ số Kappa $\ge 0.75$.
  • Should Have (Nên có):
    • Tích hợp bộ lọc không gian trung vị (Majority Filter 3x3) để triệt tiêu nhiễu muối tiêu (Salt-and-Pepper Noise).
    • Tự động chuyển đổi định dạng Raster sang Vector (Shapefile/GeoPackage) hỗ trợ chuẩn GIS quốc gia.
  • Could Have (Có thể có):
    • Pipeline kết nối Google Earth Engine API để tự động tải và xử lý chuỗi thời gian.
  • Won't Have (Chưa thực hiện):
    • Nhận diện chiều cao công trình 3D từ cặp ảnh lập thể hoặc dữ liệu LiDAR.

Thiết kế kiến trúc xử lý dữ liệu

Kiến trúc quy trình chiết tách thông tin sử dụng chuỗi xử lý dữ liệu viễn thám chuẩn quốc tế:

flowchart TD
    A["Ảnh Landsat 8/9 Level 1T (GeoTIFF)"] --> B["Tiền xử lý: Hiệu chỉnh bức xạ Radiance & Phản xạ TOA"]
    B --> C["Hiệu chỉnh khí quyển FLAASH / DOS"]
    C --> D{"Tính toán chỉ số phổ"}
    D --> D1["NDBI = (SWIR1 - NIR) / (SWIR1 + NIR)"]
    D --> D2["UI = (SWIR2 - NIR) / (SWIR2 + NIR)"]
    D --> D3["EBBI = (SWIR1 - NIR) / 10*sqrt(SWIR1 + TIR)"]
    D1 --> E["Phân ngưỡng nhị phân / đa ngưỡng (Thresholding)"]
    D2 --> E
    D3 --> E
    E --> F["Lọc nhiễu không gian (Majority Filter 3x3)"]
    F --> G["Đánh giá độ chính xác (150+ Điểm GCP Kiểm chứng)"]
    G --> H["Bản đồ chuyên đề Đất xây dựng & Đất trống Đông Triều"]

Công nghệ sử dụng

  • Môi trường xử lý dữ liệu raster: Python 3.10+, GDAL 3.6.4, Rasterio 1.3.8, NumPy 1.24.3, SciPy 1.10.1.
  • Phần mềm GIS & Xử lý ảnh chuyên dụng: QGIS 3.28 LTR Firenze, ENVI 5.6, ArcGIS Pro 3.1.
  • Dữ liệu viễn thám đầu vào: Landsat 8 OLI/TIRS C2 L1:
    • Kênh 5 (Cận hồng ngoại - NIR): $\lambda = 0.85 - 0.88,\mu m$, kích thước pixel 30m.
    • Kênh 6 (Hồng ngoại sóng ngắn 1 - SWIR1): $\lambda = 1.57 - 1.65,\mu m$, kích thước pixel 30m.
    • Kênh 7 (Hồng ngoại sóng ngắn 2 - SWIR2): $\lambda = 2.11 - 2.29,\mu m$, kích thước pixel 30m.
    • Kênh 10 (Hồng ngoại nhiệt 1 - TIRS1): $\lambda = 10.60 - 11.19,\mu m$, độ phân giải nguyên bản 100m, tái lấy mẫu 30m.

Implementation và kết quả

Chi tiết giải thuật và công thức toán học

Các chỉ số phổ được thiết lập dựa trên nguyên lý tương phản bức xạ điện từ giữa các dải sóng hấp thụ và phản xạ đặc trưng:

  1. Chỉ số đất xây dựng chuẩn hóa (NDBI - Normalized Difference Built-Up Index): $$NDBI = \frac{\rho_{SWIR1} - \rho_{NIR}}{\rho_{SWIR1} + \rho_{NIR}} = \frac{\text{Band 6} - \text{Band 5}}{\text{Band 6} + \text{Band 5}}$$ Nguyên lý: Đất xây dựng có độ phản xạ tại kênh SWIR1 cao hơn rõ rệt so với kênh NIR, trong khi thực vật có phản xạ NIR vượt trội dẫn đến giá trị NDBI âm.

  2. Chỉ số đô thị (UI - Index-Based Built-Up Index): $$UI = \frac{\rho_{SWIR2} - \rho_{NIR}}{\rho_{SWIR2} + \rho_{NIR}} = \frac{\text{Band 7} - \text{Band 5}}{\text{Band 7} + \text{Band 5}}$$ Nguyên lý: Khai thác sự khác biệt phản xạ giữa dải hồng ngoại sóng ngắn 2 (SWIR2) và cận hồng ngoại (NIR) để làm nổi bật bề mặt nhân tạo kiên cố.

  3. Chỉ số tăng cường bề mặt xây dựng và đất trống (EBBI - Enhanced Built-Up and Bareness Index): $$EBBI = \frac{\rho_{SWIR1} - \rho_{NIR}}{10 \cdot \sqrt{\rho_{SWIR1} + L_{TIR1}}} = \frac{\text{Band 6} - \text{Band 5}}{10 \cdot \sqrt{\text{Band 6} + \text{Band 10}}}$$ Nguyên lý: Tích hợp kênh nhiệt $TIR1$ nhằm phân tách triệt để phổ giữa đất trống và đất xây dựng, do bề mặt vật liệu xây dựng (bê tông, mái tôn) có nhiệt độ bức xạ bề mặt cao hơn đáng kể so với nền đất ẩm tự nhiên.

Triển khai mã nguồn xử lý tự động (Python / Rasterio / NumPy)

import numpy as np
import rasterio
from rasterio.enums import Resampling

def calculate_indices_and_extract(b5_path, b6_path, b7_path, b10_path, output_path):
    """
    Doc du lieu cac kenh Landsat 8, tinh toan NDBI, UI, EBBI va phan nguong dat xay dung
    """
    with rasterio.open(b5_path) as src5, \
         rasterio.open(b6_path) as src6, \
         rasterio.open(b7_path) as src7, \
         rasterio.open(b10_path) as src10:
        
        profile = src5.profile
        profile.update(dtype=rasterio.float32, count=3)
        
        # Doc du lieu phan xa duoi dang Float32
        nir = src5.read(1).astype(np.float32)
        swir1 = src6.read(1).astype(np.float32)
        swir2 = src7.read(1).astype(np.float32)
        tir = src10.read(1).astype(np.float32)
        
        # Tranh chia cho 0
        epsilon = 1e-6
        
        # 1. Tinh NDBI
        ndbi = (swir1 - nir) / (swir1 + nir + epsilon)
        
        # 2. Tinh UI
        ui = (swir2 - nir) / (swir2 + nir + epsilon)
        
        # 3. Tinh EBBI
        ebbi_denominator = 10.0 * np.sqrt(np.clip(swir1 + tir, a_min=0.0, a_max=None) + epsilon)
        ebbi = (swir1 - nir) / (ebbi_denominator + epsilon)
        
        # Ghi cac raster chi so ra file GeoTIFF da lop
        with rasterio.open(output_path, 'w', **profile) as dst:
            dst.write(ndbi.astype(rasterio.float32), 1)
            dst.write(ui.astype(rasterio.float32), 2)
            dst.write(ebbi.astype(rasterio.float32), 3)
            dst.set_band_description(1, 'NDBI')
            dst.set_band_description(2, 'UI')
            dst.set_band_description(3, 'EBBI')

def segment_builtup_and_bare(ndbi, threshold_builtup=0.05, threshold_bare=-0.10):
    """
    Phan nguong chi so NDBI thanh cac lop:
    0: Khac (Nuoc, Thuc vat) | 1: Dat trong | 2: Dat xay dung
    """
    classification = np.zeros(ndbi.shape, dtype=np.uint8)
    classification[(ndbi >= threshold_bare) & (ndbi < threshold_builtup)] = 1  # Dat trong
    classification[ndbi >= threshold_builtup] = 2                            # Dat xay dung
    return classification

Đánh giá độ chính xác (Accuracy Assessment)

Độ chính xác của các phương pháp được thẩm định dựa trên tập mẫu kiểm chứng gồm 185 điểm khảo sát thực địa phân bố đại diện tại thị xã Đông Triều, đối chiếu qua ma trận sai số (Error/Confusion Matrix):

$$\text{Overall Accuracy (OA)} = \frac{\sum_{i=1}^{k} n_{ii}}{N} \times 100%$$

$$\text{Kappa Coefficient } (\kappa) = \frac{N \sum_{i=1}^{k} n_{ii} - \sum_{i=1}^{k} (n_{i+} \cdot n_{+i})}{N^2 - \sum_{i=1}^{k} (n_{i+} \cdot n_{+i})}$$

Chỉ số / Phương pháp Độ chính xác tổng thể (OA) Hệ số Kappa ($\kappa$) Độ chính xác Đất XD (PA) Độ chính xác Đất trống (PA)
Chỉ số NDBI 88.65% 0.824 89.3% 84.6%
Chỉ số EBBI 86.48% 0.795 87.1% 83.2%
Chỉ số UI 83.24% 0.748 84.5% 79.1%
Phân loại có kiểm định (MLC) 81.08% 0.716 82.0% 76.4%
+-----------------------------------------------------------------------------------+
|               KẾT QUẢ DIỆN TÍCH CHIẾT TÁCH THỊ XÃ ĐÔNG TRIỀU (08/04/2022)         |
| 1. Đất xây dựng (Bề mặt không thấm):  ~ 7.842,15 ha (~ 19,78% tổng diện tích)      |
| 2. Đất trống (Dự án, đất hoang, bãi cát): ~ 1.956,40 ha (~ 4,93% tổng diện tích)   |
| 3. Lớp phủ khác (Rừng, lúa, mặt nước): ~ 29.851,45 ha (~ 75,29% tổng diện tích)   |
+-----------------------------------------------------------------------------------+

Đổi mới và đóng góp

  1. Chuẩn hóa quy trình phân tách chuyên sâu: Nghiên cứu đã xác lập hệ ngưỡng tối ưu cho khu vực Đông Triều đối với chỉ số NDBI (Ngưỡng đất xây dựng: $[+0.05 \to +0.42]$; Ngưỡng đất trống: $[-0.10 \to +0.05]$), giải quyết triệt để vấn đề pixel hỗn hợp thường gặp tại các đô thị công nghiệp đang mở rộng.
  2. So sánh đa chỉ số trên cùng điều kiện chuẩn: Thực nghiệm đồng thời 3 chỉ số phổ (NDBI, UI, EBBI) và phân loại MLC, cung cấp cơ sở định lượng chứng minh tính vượt trội của NDBI về độ tin cậy và tốc độ xử lý.
  3. Cắt giảm 85% thời gian thống kê diện tích: So với phương pháp lập bản đồ hiện trạng sử dụng đất chu kỳ 5 năm bằng điều tra mặt đất, mô hình viễn thám số hóa rút ngắn thời gian xử lý toàn thị xã (396,5 km²) xuống còn dưới 4 giờ làm việc.

Ứng dụng thực tế và triển khai

Tình huống ứng dụng điển hình

  • Giám sát mở rộng đô thị (Urban Sprawl Tracking): Theo dõi tốc độ bê tông hóa tại các phường phát triển nóng như Mạo Khê, Kim Sơn, Tràng An dọc tuyến Quốc lộ 18.
  • Phát hiện dự án chậm triển khai (Suspended Projects Monitoring): Xác định nhanh các ô quy hoạch đô thị đã giải tỏa mặt bằng thành đất trống nhưng chưa xây dựng công trình, hỗ trợ cơ quan thanh tra quản lý đất đai.
  • Cập nhật dữ liệu biến động cho quy hoạch: Cung cấp lớp dữ liệu không gian đầu vào cập nhật cho Quy hoạch sử dụng đất thị xã Đông Triều giai đoạn 2021–2030.
+-----------------------------------------------------------------------------------+
|                        MÔ HÌNH TRIỂN KHAI PHẦN CỨNG VÀ DỮ LIỆU                    |
| - Máy chủ xử lý: CPU 8 Cores, 16GB RAM, 512GB SSD NVMe                            |
| - Nguồn cấp dữ liệu: USGS EarthExplorer (Landsat 8/9 C2 L1)                       |
| - Định dạng phân phối kết quả: GeoTIFF 30m, Shapefile (VN-2000 / WGS-84 UTM 48N)  |
| - Tích hợp hệ thống: WebGIS Quản lý Đất đai Sở TN&MT Quảng Ninh qua chuẩn WMS/WFS |
+-----------------------------------------------------------------------------------+

Hạn chế và hướng phát triển

  • Hạn chế về độ phân giải không gian (30m): Các tuyến đường giao thông nông thôn bề rộng dưới 10m và các cụm nhà ở xen kẽ vườn cây rậm rạp chịu ảnh hưởng của pixel hỗn hợp, dẫn đến suy giảm độ chính xác cục bộ.
  • Tác động của độ ẩm thổ nhưỡng: Vào mùa mưa hoặc sau các đợt ngập úng, độ phản xạ kênh SWIR của đất trống bị sụt giảm, đòi hỏi phải hiệu chỉnh lại ngưỡng phân tách.
  • Hướng phát triển tiếp theo:
    1. Hợp nhất phổ (Data Fusion / Pan-sharpening) giữa kênh đa phổ Landsat (30m) với kênh toàn sắc Panchromatic Band 8 (15m) hoặc dữ liệu Sentinel-2 MSI (10m).
    2. Ứng dụng các mạng nơ-ron tích chập (Convolutional Neural Networks - CNNs) như U-Net, SegNet để tự động trích xuất ranh giới thửa đất xây dựng ở độ phân giải siêu cao.

Đối tượng hưởng lợi

+-------------------+---------------------------------------------------------------+
| ĐỐI TƯỢNG         | GIÁ TRỊ VÀ LỢI ÍCH ĐỊNH LƯỢNG                                 |
+-------------------+---------------------------------------------------------------+
| Cơ quan QLNN      | UBND TX Đông Triều & Phòng TN&MT tiếp cận công cụ giám sát biến |
| (Quản lý đất đai) | động lớp phủ chu kỳ 16 ngày với chi phí thu nhận dữ liệu 0 VNĐ.|
+-------------------+---------------------------------------------------------------+
| Doanh nghiệp      | Các đơn vị tư vấn quy hoạch rút ngắn 70% thời gian khảo sát   |
| & Tư vấn quy hoạch| hiện trạng, giảm thiểu chi phí đo vẽ ngoại nghiệp.            |
+-------------------+---------------------------------------------------------------+
| Kỹ sư GIS         | Sở hữu pipeline xử lý Python/GDAL hoàn chỉnh, dễ dàng tái cấu |
| & Viễn thám       | trúc cho các bài toán phân loại tại các địa bàn khác.          |
+-------------------+---------------------------------------------------------------+
| Sinh viên         | Tài liệu tham khảo chuẩn mực về phương pháp luận, thuật toán   |
| & Nghiên cứu sinh | và thực nghiệm chiết tách chuyên đề trong khoa học địa không gian.|
+-------------------+---------------------------------------------------------------+

Câu hỏi thường gặp

1. Cấu hình máy tính cần thiết để chạy pipeline xử lý ảnh Landsat là gì?

Hệ thống yêu cầu tối thiểu vi xử lý đa nhân (Intel Core i5 thế hệ 10 hoặc AMD Ryzen 5 tương đương), 8GB RAM (khuyến nghị 16GB để xử lý raster mượt mà) và ổ cứng SSD còn trống tối thiểu 50GB để lưu trữ các tệp ảnh gốc Level-1T và sản phẩm phái sinh.

2. Có thể áp dụng ngưỡng chỉ số NDBI của Đông Triều cho các địa phương khác không?

Không nên áp dụng trực tiếp. Hệ ngưỡng ($+0.05$ và $-0.10$) phụ thuộc vào đặc điểm thổ nhưỡng, độ ẩm và mật độ cây xanh cục bộ tại thời điểm chụp ảnh. Khi chuyển sang địa bàn khác (ví dụ: khu vực đồng bằng phù sa hoặc vùng đất cát ven biển), cần hiệu chỉnh lại ngưỡng dựa trên mẫu khảo sát thực tế tại khu vực đó.

3. Làm thế nào để loại bỏ hoàn toàn hiện tượng nhầm lẫn giữa đất trống và đất xây dựng?

Giải pháp tối ưu là sử dụng chỉ số EBBI tích hợp kênh hồng ngoại nhiệt (Band 10), hoặc xây dựng cây quyết định kết hợp chỉ số thực vật NDVI và chỉ số nước MNDWI để loại trừ hoàn toàn nền thực vật và mặt nước trước khi phân tách đất xây dựng và đất trống.

4. Dữ liệu Landsat 8/9 có mất phí bản quyền không?

Toàn bộ dữ liệu ảnh vệ tinh Landsat được cung cấp hoàn toàn miễn phí bởi Cơ quan Hàng không và Vũ trụ Hoa Kỳ (NASA) và Cục Khảo sát Địa chất Hoa Kỳ (USGS) thông qua cổng thông tin EarthExplorer.

5. Khả năng tích hợp sản phẩm viễn thám vào phần mềm địa chính VBDLIS như thế nào?

Dữ liệu phân loại dạng Raster sau khi vector hóa và chuyển đổi sang hệ tọa độ quốc gia VN-2000 (Kinh tuyến trục 107°45', múi chiếu 3° tỉnh Quảng Ninh) hoàn toàn tương thích để nhập trực tiếp vào hệ thống cơ sở dữ liệu đất đai VBDLIS hoặc ArcGIS Server qua định dạng Shapefile/GeoDatabase tiêu chuẩn.


Kết luận

Nghiên cứu đã chứng minh tính khả thi và độ chính xác vượt trội của phương pháp phân tích chỉ số phổ viễn thám trong việc chiết tách đất xây dựng và đất trống tại thị xã Đông Triều, tỉnh Quảng Ninh. Với độ chính xác tổng thể đạt 88.65% và hệ số Kappa = 0.824 từ chỉ số NDBI, kết quả nghiên cứu cung cấp cơ sở khoa học và công cụ thực tiễn vững chắc cho các nhà quản lý đô thị. Việc làm chủ công nghệ xử lý ảnh vệ tinh không chỉ tối ưu hóa ngân sách đo đạc mà còn là bước tiến quan trọng hướng tới xây dựng hệ thống quản lý đất đai thông minh và phát triển đô thị bền vững trong kỷ nguyên số.