Полное руководство: Расчет коэффициента R² в машинном обучении дистанционного зондирования

Расчет значений R² для анализа сельскохозяйственных угодий с использованием спутниковых данных и машинного обучения

Для расчета значений R² при анализе сельскохозяйственных угодий с использованием спутниковых данных и машинного обучения необходимо реализовать коэффициент детерминации между прогнозами модели и фактическими значениями с помощью функции r2_score из библиотеки scikit-learn или ручных вычислений с использованием numpy. Для обработки спутниковых изображений с использованием rasterio и xgboost это включает извлечение спектральных индексов и признаков из растровых каналов, обучение модели на размеченных данных ферм, а затем оценку точности прогнозирования с использованием R² в качестве ключевого метрического показателя валидации, который указывает на долю дисперсии целевой переменной, объясненную вашей моделью.

Содержание

• Понимание R² в контексте дистанционного зондирования
• Подготовка спутниковых данных с помощью Rasterio
• Инженерия признаков для анализа сельскохозяйственных угодий
• Обучение модели машинного обучения с XGBoost
• Расчет значений R²
• Валидация и интерпретация
• Практический пример реализации

Понимание R² в контексте дистанционного зондирования

R-квадрат (R²), также известный как коэффициент детерминации, является фундаментальным статистическим показателем в проектах машинного обучения с использованием дистанционного зондирования. В контексте анализа сельскохозяйственных угодий R² представляет собой долю дисперсии целевой переменной (такой как урожайность, влажность почвы или здоровье растительности), которая может быть объяснена признаками, полученными со спутниковых данных. Согласно исследованиям в области дистанционного зондирования, R² рассчитывается как квадрат коэффициента корреляции (CC) и служит ключевым метрическим показателем для выбора наилучших комбинаций моделей машинного обучения [источник].

В приложениях для анализа сельскохозяйственных угодий значения R² обычно варьируются от 0 до 1, где:

• R² = 1 указывает на идеальное прогнозирование (вся дисперсия объяснена)
• R² = 0 указывает на то, что модель работает не лучше, чем среднее значение
• Отрицательный R² предполагает, что модель работает хуже, чем использование среднего значения

Важность R² в сельскохозяйственном дистанционном зондировании обусловлена его способностью количественно оценивать, насколько хорошо ваша модель машинного обучения улавливает связь между признаками спутниковых изображений и реальными результатами в сельском хозяйстве. Это делает его особенно ценным для прогнозирования урожайности, оценки здоровья почвы и мониторинга растительности с использованием мультиспектральных спутниковых данных.

Подготовка спутниковых данных с помощью Rasterio

При работе со спутниковыми изображениями для анализа сельскохозяйственных угодий rasterio является основной библиотекой для чтения, обработки и извлечения данных из геопространственных растровых файлов. Каждый растровый набор данных может содержать несколько каналов, и каждый канал по сути обрабатывается как массив NumPy, что позволяет выполнять эффективные численные вычисления [источник].

import rasterio
import numpy as np

def load_satellite_data(file_path):
    """Загрузка спутниковых изображений с помощью rasterio"""
    with rasterio.open(file_path) as src:
        # Чтение всех каналов
        bands = src.read()
        # Получение метаданных, включая преобразование и CRS
        transform = src.transform
        crs = src.crs
        return bands, transform, crs

Основные этапы предварительной обработки для анализа сельскохозяйственных угодий включают:

1. Очистка данных: обработка пропущенных значений (маскирование облаков, пиксели без данных)
2. Геометрическая коррекция: обеспечение правильного выравнивания между разными растровыми слоями
3. Атмосферная коррекция: преобразование цифровых чисел в значения отражательной способности
4. Извлечение подмножеств: фокусировка на регионах сельскохозяйственного интереса

Библиотека rasterio предоставляет мощные возможности для этих операций, позволяя эффективно работать с большими наборами спутниковых данных, сохраняя при этом геопространственную справочную информацию. Это критически важно для анализа сельскохозяйственных угодий, где пространственные отношения и точная геопривязка являются существенными.

Инженерия признаков для анализа сельскохозяйственных угодий

Эффективная инженерия признаков является критически важной при работе со спутниковыми изображениями для сельскохозяйственных приложений. Данные дистанционного зондирования содержат богатую информацию, которая может быть преобразована в значимые признаки для моделей машинного обучения. Распространенные подходы включают:

Спектральные индексы

Наиболее широко используемым вегетационным индексом в дистанционном зондировании является нормализованный разностный вегетационный индекс (NDVI), рассчитываемый с использованием ближнего инфракрасного и красного диапазонов длин волн для измерения зеленой массы или здоровья растительности [источник]. Другие важные индексы для анализа сельскохозяйственных угодий включают:

def calculate_indices(red, nir, swir1):
    """Расчет распространенных вегетационных и сельскохозяйственных индексов"""
    ndvi = (nir - red) / (nir + red)
    evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * swir1 + 1)
    savi = (nir - red) / (nir + red + 0.5) * (1 + 0.5)
    return ndvi, evi, savi

Текстурные и статистические признаки

• Среднее значение, стандартное отклонение: спектральная изменчивость внутри полей
• Процентили: характеристика распределения
• Метрики текстуры: признаки матрицы совместного встречения уровней серого (GLCM)

Временные признаки

Для многовременного анализа следует учитывать:

• Фенологические метрики: сезонные паттерны роста растительности
• Обнаружение изменений: различия между временными периодами
• Кривые роста: временная эволюция вегетационных индексов

Эти признаки служат основой для надежных моделей машинного обучения в анализе сельскохозяйственных угодий, улавливая как пространственные, так и временные паттерны в сельскохозяйственных ландшафтах.

Обучение модели машинного обучения с XGBoost

XGBoost (Extreme Gradient Boosting) особенно хорошо подходит для анализа сельскохозяйственных угодий с использованием спутниковых данных благодаря своей способности обрабатывать сложные отношения в многомерном пространстве признаков. При работе с растровыми данными XGBoost может эффективно обучаться на разнообразных спектральных, текстурных и временных признаках, извлеченных со спутниковых изображений [источник].

Подготовка данных для XGBoost

import pandas as pd
from sklearn.model_selection import train_test_split

# Подготовка данных для XGBoost
def prepare_xgboost_data(features_df, target_variable):
    """
    Подготовка данных для обучения XGBoost
    features_df: DataFrame, содержащий все созданные признаки
    target_variable: Название целевого столбца (например, 'yield', 'health_score')
    """
    X = features_df.drop(target_variable, axis=1)
    y = features_df[target_variable]
    
    # Разделение данных
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.25, random_state=42
    )
    
    return X_train, X_test, y_train, y_test

Конфигурация обучения XGBoost

import xgboost as xgb

def train_xgboost_model(X_train, y_train, params=None):
    """Обучение модели XGBoost для анализа сельскохозяйственных угодий"""
    if params is None:
        params = {
            'objective': 'reg:squarederror',
            'max_depth': 6,
            'learning_rate': 0.1,
            'n_estimators': 100,
            'subsample': 0.8,
            'colsample_bytree': 0.8,
            'random_state': 42
        }
    
    model = xgb.XGBRegressor(**params)
    model.fit(X_train, y_train)
    return model

Возможности регуляризации XGBoost помогают предотвратить переобучение при работе с потенциально зашумленными спутниковыми данными, что делает его идеальным для сельскохозяйственных приложений, где данные полевых измерений могут иметь ограничения.

Расчет значений R²

Расчет значений R² является основным этапом валидации для ваших моделей машинного обучения в анализе сельскохозяйственных угодий. Существует несколько подходов к вычислению R² при работе со спутниковыми изображениями и моделями xgboost.

Метод 1: Использование Scikit-learn (Рекомендуется)

from sklearn.metrics import r2_score

def calculate_r2(model, X_test, y_test):
    """Расчет оценки R² с помощью scikit-learn"""
    predictions = model.predict(X_test)
    r2 = r2_score(y_test, predictions)
    return r2, predictions

# Использование
r2_score, predictions = calculate_r2(xgboost_model, X_test, y_test)
print(f"Оценка R²: {r2_score:.4f}")

Метод 2: Ручной расчет с использованием NumPy

def calculate_r2_manual(y_true, y_pred):
    """
    Ручной расчет R² с использованием numpy
    Формула: R² = 1 - (SS_res / SS_tot)
    где SS_res = Σ(y_true - y_pred)²
    и SS_tot = Σ(y_true - mean(y_true))²
    """
    y_true = np.array(y_true)
    y_pred = np.array(y_pred)
    
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    
    r2 = 1 - (ss_res / ss_tot)
    return r2

Метод 3: Расчет R² на пиксель для пространственного анализа

При работе с растровыми данными может потребоваться вычисление значений R² пространственно:

def calculate_spatial_r2(observed_raster, predicted_raster):
    """Расчет R² для каждого пикселя по нескольким временным периодам"""
    # Выравнивание массивов для расчета
    obs_flat = observed_raster.flatten()
    pred_flat = predicted_raster.flatten()
    
    # Расчет R²
    ss_res = np.sum((obs_flat - pred_flat) ** 2)
    ss_tot = np.sum((obs_flat - np.mean(obs_flat)) ** 2)
    r2 = 1 - (ss_res / ss_tot)
    
    return r2

Для анализа сельскохозяйственных угодий часто бывает ценно рассчитывать значения R² для разных типов культур, условий почвы или методов управления, чтобы понять производительность модели в различных сельскохозяйственных контекстах.

Валидация и интерпретация

Стратегии кросс-валидации

При работе со спутниковыми данными для анализа сельскохозяйственных угодий надежная валидация является обязательной из-за пространственной автокорреляции и потенциальной смещенности выборки. Рассмотрите следующие подходы:

from sklearn.model_selection import cross_val_score, GroupKFold

def cross_validate_model(model, X, y, groups=None):
    """Выполнение кросс-валидации для моделей анализа сельскохозяйственных угодий"""
    if groups is not None:
        # Использование GroupKFold для обеспечения того, чтобы поля не разделялись между обучением и тестом
        cv = GroupKFold(n_splits=5)
    else:
        cv = 5
    
    scores = cross_val_score(
        model, X, y, cv=cv, 
        scoring='r2', n_jobs=-1
    )
    
    return scores

Интерпретация значений R² в сельскохозяйственном контексте

Интерпретация R² для анализа сельскохозяйственных угодий должна учитывать специфические факторы сельского хозяйства:

• Высокий R² (>0.8): Отличная производительность модели, подходящая для применения точного земледелия
• Умеренный R² (0.5-0.8): Полезен для общего сельскохозяйственного мониторинга и планирования
• Низкий R² (<0.5): Может потребоваться дополнительные признаки или другие подходы к моделированию

Сельскохозяйственные системы по своей природе сложны и зависят от погодных условий, почвы, методов управления и других факторов. Поэтому даже умеренные значения R² могут быть ценны для практического управления фермой в сочетании с другими метриками и отраслевыми знаниями.

Дополнительные метрики оценки

Дополните R² другими метриками для комплексного анализа сельскохозяйственных угодий:

from sklearn.metrics import mean_squared_error, mean_absolute_error

def evaluate_model_performance(y_true, y_pred):
    """Комплексная оценка производительности модели для анализа сельскохозяйственных угодий"""
    r2 = r2_score(y_true, y_pred)
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    mae = mean_absolute_error(y_true, y_pred)
    
    return {
        'R²': r2,
        'RMSE': rmse,
        'MAE': mae,
        'RMSE_relative': rmse / np.mean(y_true) * 100
    }

Практический пример реализации

Вот полная рабочая流程 для расчета значений R² в проекте анализа сельскохозяйственных угодий с использованием указанного технологического стека:

import rasterio
import numpy as np
import pandas as pd
import xgboost as xgb
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split

class FarmlandAnalyzer:
    def __init__(self):
        self.model = None
        self.feature_names = None
        
    def load_satellite_data(self, file_paths):
        """Загрузка и предварительная обработка спутниковых изображений"""
        all_bands = []
        for file_path in file_paths:
            with rasterio.open(file_path) as src:
                bands = src.read()
                all_bands.append(bands)
        
        # Стекинг каналов со всех дат
        stacked_bands = np.concatenate(all_bands, axis=0)
        return stacked_bands
    
    def extract_features(self, stacked_bands, field_boundaries):
        """Извлечение признаков для каждого сельскохозяйственного поля"""
        features = []
        
        for field_id, bounds in field_boundaries.items():
            # Извлечение растровых данных для этого поля
            field_data = self._extract_field_data(stacked_bands, bounds)
            
            # Расчет признаков
            field_features = self._calculate_field_features(field_data)
            field_features['field_id'] = field_id
            features.append(field_features)
        
        return pd.DataFrame(features)
    
    def _extract_field_data(self, stacked_bands, bounds):
        """Извлечение растровых данных в пределах границ поля"""
        # Реализация зависит от вашей конкретной системы координат и формата границ
        pass
    
    def _calculate_field_features(self, field_data):
        """Расчет спектральных и статистических признаков для поля"""
        features = {}
        
        # Расчет индексов для каждого временного периода
        num_time_periods = field_data.shape[0] // 4  # Предполагаем 4 канала на период
        
        for t in range(num_time_periods):
            start_idx = t * 4
            red = field_data[start_idx + 2]  # Обычно красный канал
            nir = field_data[start_idx + 3]  # Обычно ближний ИК канал
            swir1 = field_data[start_idx + 4] if field_data.shape[0] > start_idx + 4 else None
            
            if swir1 is not None:
                ndvi, evi, savi = self._calculate_indices(red, nir, swir1)
            else:
                ndvi = self._calculate_ndvi(red, nir)
                evi, savi = None, None
            
            # Хранение временно-специфичных признаков
            features[f'ndvi_t{t}'] = np.mean(ndvi)
            features[f'ndvi_std_t{t}'] = np.std(ndvi)
            if evi is not None:
                features[f'evi_t{t}'] = np.mean(evi)
                features[f'savi_t{t}'] = np.mean(savi)
        
        return features
    
    def _calculate_indices(self, red, nir, swir1):
        """Расчет вегетационных индексов"""
        ndvi = (nir - red) / (nir + red + 1e-8)
        evi = 2.5 * (nir - red) / (nir + 6 * red - 7.5 * swir1 + 1 + 1e-8)
        savi = (nir - red) / (nir + red + 0.5) * (1 + 0.5)
        return ndvi, evi, savi
    
    def _calculate_ndvi(self, red, nir):
        """Расчет NDVI"""
        return (nir - red) / (nir + red + 1e-8)
    
    def train_model(self, features_df, target_col):
        """Обучение модели XGBoost для анализа сельскохозяйственных угодий"""
        X = features_df.drop(target_col, axis=1)
        y = features_df[target_col]
        
        # Разделение данных
        X_train, X_test, y_train, y_test = train_test_split(
            X, y, test_size=0.25, random_state=42
        )
        
        # Обучение модели
        params = {
            'objective': 'reg:squarederror',
            'max_depth': 6,
            'learning_rate': 0.1,
            'n_estimators': 100,
            'subsample': 0.8,
            'colsample_bytree': 0.8,
            'random_state': 42
        }
        
        self.model = xgb.XGBRegressor(**params)
        self.model.fit(X_train, y_train)
        self.feature_names = X.columns.tolist()
        
        # Расчет R²
        predictions = self.model.predict(X_test)
        r2 = r2_score(y_test, predictions)
        
        return {
            'model': self.model,
            'r2_score': r2,
            'feature_importance': dict(zip(
                self.feature_names, 
                self.model.feature_importances_
            ))
        }
    
    def predict_r2(self, new_features_df):
        """Прогнозирование и расчет R² для новых данных"""
        if self.model is None:
            raise ValueError("Модель еще не обучена")
        
        predictions = self.model.predict(new_features_df)
        
        # Если доступны истинные значения, расчет R²
        if 'target' in new_features_df.columns:
            true_values = new_features_df['target']
            r2 = r2_score(true_values, predictions)
            return predictions, r2
        else:
            return predictions, None

Пример использования

# Инициализация анализатора
analyzer = FarmlandAnalyzer()

# Загрузка спутниковых данных (пример путей к файлам)
satellite_files = [
    'landsat_2021_05.tif',
    'landsat_2021_06.tif',
    'landsat_2021_07.tif'
]

stacked_bands = analyzer.load_satellite_data(satellite_files)

# Определение границ полей (пример)
field_boundaries = {
    'field_1': {'minx': 1000, 'maxx': 1200, 'miny': 2000, 'maxy': 2200},
    'field_2': {'minx': 1500, 'maxx': 1700, 'miny': 2500, 'maxy': 2700}
}

# Извлечение признаков
features_df = analyzer.extract_features(stacked_bands, field_boundaries)

# Добавление целевой переменной (пример: урожайность)
features_df['yield'] = np.random.normal(5.0, 1.0, len(features_df))  # Заменить на реальные данные

# Обучение модели и получение R²
result = analyzer.train_model(features_df, 'yield')
print(f"Оценка R² модели: {result['r2_score']:.4f}")

# Анализ важности признаков
print("Топ-5 наиболее важных признаков:")
sorted_features = sorted(
    result['feature_importance'].items(), 
    key=lambda x: x[1], 
    reverse=True
)
for feature, importance in sorted_features[:5]:
    print(f"{feature}: {importance:.4f}")

Эта комплексная реализация предоставляет полную рабочую流程 для расчета значений R² в анализе сельскохозяйственных угодий с использованием спутниковых данных и указанного технологического стека.

Источники

1. Количественная оценка осадков по данным дистанционного зондирования с использованием регрессионных моделей машинного обучения - Исследование использования R-квадрата для выбора модели дистанционного зондирования
2. R - использование случайных лесов, опорных векторов и нейронных сетей для пиксельной классификации мультиспектральных изображений Sentinel-2 - Подходы машинного обучения для спутниковых изображений
3. Python-ориентированная рабочая流程 для классификации земного покрова с использованием Geemap, Rasterio, Geopandas, Numpy, Matplotlib и Scikit-learn - Практическая реализация на Python для обработки спутниковых данных
4. Работа с растровыми данными — Документация по геопространственному анализу с Python и R 2020 - Технические детали обработки растровых данных
5. Анализ изменений земного покрова с помощью Python и Rasterio - Учебное пособие - Комплексное учебное пособие по обработке растровых данных
6. XGBoost для классификации земного покрова в R - GIS-Blog.com - Реализация XGBoost для геопространственных данных
7. Rasterio — ThirdEye Data - Документация и возможности библиотеки Rasterio
8. Обработка спутниковых изображений с помощью Python и R с использованием данных Landsat 9 OLI/TIRS и SRTM DEM в Кот-д’Ивуаре, Западная Африка - Продвинутые техники обработки спутниковых изображений

Заключение

Расчет значений R² для анализа сельскохозяйственных угодий с использованием спутниковых данных и машинного обучения требует системного подхода, который объединяет обработку растровых данных, инженерию признаков и правильную оценку модели. Основные выводы включают:

1. Используйте функцию r2_score из scikit-learn для надежного расчета R² или реализуйте ручные расчеты с использованием numpy при работе с пользовательскими метриками оценки.

2. Правильная предварительная обработка данных с помощью rasterio является обязательной - убедитесь, что ваши спутниковые изображения правильно геопривязаны, атмосферно скорректированы и правильно выровнены перед извлечением признаков.

3. Инженерия признаков должна включать спектральные индексы (NDVI, EVI, SAVI), статистические показатели и временные признаки для улавливания сложных отношений в сельскохозяйственных системах.

4. Модели XGBoost хорошо работают для анализа сельскохозяйственных угодий благодаря своей способности обрабатывать многомерные пространства признаков и обеспечивать надежные прогнозы даже с зашумленными спутниковыми данными.

5. Стратегии валидации должны учитывать пространственную автокорреляцию в сельскохозяйственных данных, используя такие техники, как GroupKFold для надежной оценки модели.

6. Интерпретируйте значения R² в контексте сложности сельскохозяйственных систем - даже умеренные значения R² могут быть ценны для практического управления фермой в сочетании с другими метриками и отраслевыми знаниями.

Для успешной реализации начните с простой базовой модели и постепенно увеличивайте сложность по мере валидации каждого компонента. Помните, что сельскохозяйственные системы зависят от множества факторов, которые не могут быть полностью уловлены спутниковыми изображениями, поэтому значения R² должны интерпретироваться вместе с другими показателями производительности и полевыми знаниями.