Линейная регрессия

Earth Engine имеет несколько методов выполнения линейной регрессии с использованием редукторов:

  • ee.Reducer.linearFit()
  • ee.Reducer.linearRegression()
  • ee.Reducer.robustLinearRegression()
  • ee.Reducer.ridgeRegression()

Самый простой преобразователь линейной регрессии — это linearFit() , которая вычисляет оценку методом наименьших квадратов линейной функции одной переменной с постоянным членом. Для более гибкого подхода к линейному моделированию используйте один из редукторов линейной регрессии, которые допускают переменное количество независимых и зависимых переменных. linearRegression() реализует обычную регрессию наименьших квадратов (OLS). robustLinearRegression() использует функцию стоимости, основанную на остатках регрессии, для итеративного снижения веса выбросов в данных ( О'Лири, 1990 ). ridgeRegression() выполняет линейную регрессию с регуляризацией L2.

Регрессионный анализ с помощью этих методов подходит для сокращения объектов ee.ImageCollection , ee.Image , ee.FeatureCollection и ee.List . Следующие примеры демонстрируют применение каждого из них. Обратите внимание, что linearRegression() , robustLinearRegression() и ridgeRegression() имеют одинаковые входные и выходные структуры, но linearFit() ожидает двухполосный входной сигнал (X, за которым следует Y), а ridgeRegression() имеет дополнительный параметр ( lambda , необязательный ) и выходной сигнал ( pValue ).

ee.ImageCollection

linearFit()

Данные должны быть настроены как двухканальное входное изображение, где первый канал является независимой переменной, а второй канал — зависимой переменной. В следующем примере показана оценка линейного тренда будущих осадков (после 2006 года по данным NEX-DCP30 ), прогнозируемого с помощью климатических моделей. Зависимая переменная — это прогнозируемые осадки, а независимая переменная — это время, добавляемое перед вызовом linearFit() :

Редактор кода (JavaScript) 

// This function adds a time band to the image.
var createTimeBand = function(image) {
  // Scale milliseconds by a large constant to avoid very small slopes
  // in the linear regression output.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};

// Load the input image collection: projected climate data.
var collection = ee.ImageCollection('NASA/NEX-DCP30_ENSEMBLE_STATS')
  .filter(ee.Filter.eq('scenario', 'rcp85'))
  .filterDate(ee.Date('2006-01-01'), ee.Date('2050-01-01'))
  // Map the time band function over the collection.
  .map(createTimeBand);

// Reduce the collection with the linear fit reducer.
// Independent variable are followed by dependent variables.
var linearFit = collection.select(['system:time_start', 'pr_mean'])
  .reduce(ee.Reducer.linearFit());

// Display the results.
Map.setCenter(-100.11, 40.38, 5);
Map.addLayer(linearFit,
  {min: 0, max: [-0.9, 8e-5, 1], bands: ['scale', 'offset', 'scale']}, 'fit');

Обратите внимание, что выходные данные содержат две полосы: «смещение» (пересечение) и «масштаб» («масштаб» в этом контексте относится к наклону линии, и его не следует путать с входным параметром масштаба для многих редукторов, который является пространственным масштабом). Результат с областями возрастающего тренда, выделенными синим цветом, нисходящего тренда красным и отсутствием тренда зеленым цветом, должен выглядеть примерно так, как показано на рисунке 1.


Рисунок 1. Выходные данные linearFit() примененные к прогнозируемым осадкам. Районы, где прогнозируется увеличение количества осадков, показаны синим цветом, а уменьшение количества осадков - красным.

linearRegression()

Например, предположим, что имеются две зависимые переменные: осадки и максимальная температура, а также две независимые переменные: константа и время. Коллекция идентична предыдущему примеру, но перед сокращением необходимо вручную добавить константную полосу. Первые две полосы входных данных — это переменные «X» (независимые), а следующие две полосы — это переменные «Y» (зависимые). В этом примере сначала получите коэффициенты регрессии, затем сгладьте изображение массива, чтобы извлечь интересующие полосы:

Редактор кода (JavaScript) 

// This function adds a time band to the image.
var createTimeBand = function(image) {
  // Scale milliseconds by a large constant.
  return image.addBands(image.metadata('system:time_start').divide(1e18));
};

// This function adds a constant band to the image.
var createConstantBand = function(image) {
  return ee.Image(1).addBands(image);
};

// Load the input image collection: projected climate data.
var collection = ee.ImageCollection('NASA/NEX-DCP30_ENSEMBLE_STATS')
  .filterDate(ee.Date('2006-01-01'), ee.Date('2099-01-01'))
  .filter(ee.Filter.eq('scenario', 'rcp85'))
  // Map the functions over the collection, to get constant and time bands.
  .map(createTimeBand)
  .map(createConstantBand)
  // Select the predictors and the responses.
  .select(['constant', 'system:time_start', 'pr_mean', 'tasmax_mean']);

// Compute ordinary least squares regression coefficients.
var linearRegression = collection.reduce(
  ee.Reducer.linearRegression({
    numX: 2,
    numY: 2
}));

// Compute robust linear regression coefficients.
var robustLinearRegression = collection.reduce(
  ee.Reducer.robustLinearRegression({
    numX: 2,
    numY: 2
}));

// The results are array images that must be flattened for display.
// These lists label the information along each axis of the arrays.
var bandNames = [['constant', 'time'], // 0-axis variation.
                 ['precip', 'temp']]; // 1-axis variation.

// Flatten the array images to get multi-band images according to the labels.
var lrImage = linearRegression.select(['coefficients']).arrayFlatten(bandNames);
var rlrImage = robustLinearRegression.select(['coefficients']).arrayFlatten(bandNames);

// Display the OLS results.
Map.setCenter(-100.11, 40.38, 5);
Map.addLayer(lrImage,
  {min: 0, max: [-0.9, 8e-5, 1], bands: ['time_precip', 'constant_precip', 'time_precip']}, 'OLS');

// Compare the results at a specific point:
print('OLS estimates:', lrImage.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: ee.Geometry.Point([-96.0, 41.0]),
  scale: 1000
}));

print('Robust estimates:', rlrImage.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: ee.Geometry.Point([-96.0, 41.0]),
  scale: 1000
}));

Просмотрите результаты и обнаружите, что выходные данные linearRegression() эквивалентны коэффициентам, оцененным редуктором linearFit() , хотя выходные данные linearRegression() также содержат коэффициенты для другой зависимой переменной, tasmax_mean . Коэффициенты устойчивой линейной регрессии отличаются от оценок OLS. В примере сравниваются коэффициенты разных методов регрессии в определенной точке.

ee.Image

В контексте объекта ee.Image редукторы регрессии можно использовать вместе с reduceRegion или reduceRegions для выполнения линейной регрессии на пикселях в регионе(ах). Следующие примеры демонстрируют, как рассчитать коэффициенты регрессии между полосами Landsat в произвольном многоугольнике.

linearFit()

В разделе руководства, описывающем диаграммы данных массива, показана точечная диаграмма корреляции между диапазонами Landsat 8 SWIR1 и SWIR2. Здесь рассчитываются коэффициенты линейной регрессии для этой зависимости. Возвращается словарь, содержащий свойства 'offset' (пересечение оси Y) и 'scale' (наклон).

Редактор кода (JavaScript) 

// Define a rectangle geometry around San Francisco.
var sanFrancisco = ee.Geometry.Rectangle([-122.45, 37.74, -122.4, 37.8]);

// Import a Landsat 8 TOA image for this region.
var img = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318');

// Subset the SWIR1 and SWIR2 bands. In the regression reducer, independent
// variables come first followed by the dependent variables. In this case,
// B5 (SWIR1) is the independent variable and B6 (SWIR2) is the dependent
// variable.
var imgRegress = img.select(['B5', 'B6']);

// Calculate regression coefficients for the set of pixels intersecting the
// above defined region using reduceRegion with ee.Reducer.linearFit().
var linearFit = imgRegress.reduceRegion({
  reducer: ee.Reducer.linearFit(),
  geometry: sanFrancisco,
  scale: 30,
});

// Inspect the results.
print('OLS estimates:', linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

linearRegression()

Здесь применяется тот же анализ, что и в предыдущем разделе linearFit , за исключением того, что на этот раз используется функция ee.Reducer.linearRegression . Обратите внимание, что регрессионное изображение строится из трех отдельных изображений: постоянного изображения и изображений, представляющих полосы SWIR1 и SWIR2 из одного и того же изображения Landsat 8. Имейте в виду, что вы можете объединить любой набор полос, чтобы создать входное изображение для уменьшения региона с помощью ee.Reducer.linearRegression , они не обязательно должны принадлежать одному и тому же исходному изображению.

Редактор кода (JavaScript) 

// Define a rectangle geometry around San Francisco.
var sanFrancisco = ee.Geometry.Rectangle([-122.45, 37.74, -122.4, 37.8]);

// Import a Landsat 8 TOA image for this region.
var img = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318');

// Create a new image that is the concatenation of three images: a constant,
// the SWIR1 band, and the SWIR2 band.
var constant = ee.Image(1);
var xVar = img.select('B5');
var yVar = img.select('B6');
var imgRegress = ee.Image.cat(constant, xVar, yVar);

// Calculate regression coefficients for the set of pixels intersecting the
// above defined region using reduceRegion. The numX parameter is set as 2
// because the constant and the SWIR1 bands are independent variables and they
// are the first two bands in the stack; numY is set as 1 because there is only
// one dependent variable (SWIR2) and it follows as band three in the stack.
var linearRegression = imgRegress.reduceRegion({
  reducer: ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
  }),
  geometry: sanFrancisco,
  scale: 30,
});

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var b0 = ee.List(coefList.get(0)).get(0); // y-intercept
var b1 = ee.List(coefList.get(1)).get(0); // slope

// Extract the residuals.
var residuals = ee.Array(linearRegression.get('residuals')).toList().get(0);

// Inspect the results.
print('OLS estimates', linearRegression);
print('y-intercept:', b0);
print('Slope:', b1);
print('Residuals:', residuals);

Возвращается словарь, содержащий свойства 'coefficients' и 'residuals' . Свойство 'coefficients' представляет собой массив размеров (numX, numY); каждый столбец содержит коэффициенты для соответствующей зависимой переменной. В этом случае массив состоит из двух строк и одного столбца; первая строка, первый столбец — это точка пересечения оси Y, вторая строка, первый столбец — это наклон. Свойство 'residuals' представляет собой вектор среднеквадратических остатков каждой зависимой переменной. Извлеките коэффициенты, преобразуя результат в виде массива, а затем вырезая нужные элементы или преобразуя массив в список и выбирая коэффициенты по позиции индекса.

ee.FeatureCollection

Предположим, вы хотите узнать линейную зависимость между отражательной способностью Sentinel-2 и Landsat 8 SWIR1. В этом примере для расчета взаимосвязи используется случайная выборка пикселей, отформатированная как набор точек. Создается диаграмма рассеяния пар пикселей вместе с линией наилучшего соответствия по методу наименьших квадратов (рис. 2).

Редактор кода (JavaScript) 

// Import a Sentinel-2 TOA image.
var s2ImgSwir1 = ee.Image('COPERNICUS/S2/20191022T185429_20191022T185427_T10SEH');

// Import a Landsat 8 TOA image from 12 days earlier than the S2 image.
var l8ImgSwir1 = ee.Image('LANDSAT/LC08/C02/T1_TOA/LC08_044033_20191010');

// Get the intersection between the two images - the area of interest (aoi).
var aoi = s2ImgSwir1.geometry().intersection(l8ImgSwir1.geometry());

// Get a set of 1000 random points from within the aoi. A feature collection
// is returned.
var sample = ee.FeatureCollection.randomPoints({
  region: aoi,
  points: 1000
});

// Combine the SWIR1 bands from each image into a single image.
var swir1Bands = s2ImgSwir1.select('B11')
  .addBands(l8ImgSwir1.select('B6'))
  .rename(['s2_swir1', 'l8_swir1']);

// Sample the SWIR1 bands using the sample point feature collection.
var imgSamp = swir1Bands.sampleRegions({
  collection: sample,
  scale: 30
})
// Add a constant property to each feature to be used as an independent variable.
.map(function(feature) {
  return feature.set('constant', 1);
});

// Compute linear regression coefficients. numX is 2 because
// there are two independent variables: 'constant' and 's2_swir1'. numY is 1
// because there is a single dependent variable: 'l8_swir1'. Cast the resulting
// object to an ee.Dictionary for easy access to the properties.
var linearRegression = ee.Dictionary(imgSamp.reduceColumns({
  reducer: ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
  }),
  selectors: ['constant', 's2_swir1', 'l8_swir1']
}));

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var yInt = ee.List(coefList.get(0)).get(0); // y-intercept
var slope = ee.List(coefList.get(1)).get(0); // slope

// Gather the SWIR1 values from the point sample into a list of lists.
var props = ee.List(['s2_swir1', 'l8_swir1']);
var regressionVarsList = ee.List(imgSamp.reduceColumns({
  reducer: ee.Reducer.toList().repeat(props.size()),
  selectors: props
}).get('list'));

// Convert regression x and y variable lists to an array - used later as input
// to ui.Chart.array.values for generating a scatter plot.
var x = ee.Array(ee.List(regressionVarsList.get(0)));
var y1 = ee.Array(ee.List(regressionVarsList.get(1)));

// Apply the line function defined by the slope and y-intercept of the
// regression to the x variable list to create an array that will represent
// the regression line in the scatter plot.
var y2 = ee.Array(ee.List(regressionVarsList.get(0)).map(function(x) {
  var y = ee.Number(x).multiply(slope).add(yInt);
  return y;
}));

// Concatenate the y variables (Landsat 8 SWIR1 and predicted y) into an array
// for input to ui.Chart.array.values for plotting a scatter plot.
var yArr = ee.Array.cat([y1, y2], 1);

// Make a scatter plot of the two SWIR1 bands for the point sample and include
// the least squares line of best fit through the data.
print(ui.Chart.array.values({
  array: yArr,
  axis: 0,
  xLabels: x})
  .setChartType('ScatterChart')
  .setOptions({
    legend: {position: 'none'},
    hAxis: {'title': 'Sentinel-2 SWIR1'},
    vAxis: {'title': 'Landsat 8 SWIR1'},
    series: {
      0: {
        pointSize: 0.2,
        dataOpacity: 0.5,
      },
      1: {
        pointSize: 0,
        lineWidth: 2,
      }
    }
  })
);


Рисунок 2. Диаграмма рассеяния и линия линейной регрессии по методу наименьших квадратов для выборки пикселей, представляющих коэффициент отражения TOA Sentinel-2 и Landsat 8 SWIR1.

ee.List

Столбцы двумерных объектов ee.List могут быть входными данными для редукторов регрессии. Следующие примеры предоставляют простые доказательства; независимая переменная является копией зависимой переменной, производящей точку пересечения по оси Y, равную 0, и наклон, равный 1.

linearFit()

Редактор кода (JavaScript) 

// Define a list of lists, where columns represent variables. The first column
// is the independent variable and the second is the dependent variable.
var listsVarColumns = ee.List([
  [1, 1],
  [2, 2],
  [3, 3],
  [4, 4],
  [5, 5]
]);

// Compute the least squares estimate of a linear function. Note that an
// object is returned; cast it as an ee.Dictionary to make accessing the
// coefficients easier.
var linearFit = ee.Dictionary(listsVarColumns.reduce(ee.Reducer.linearFit()));

// Inspect the result.
print(linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

Транспонируйте список, если переменные представлены строками, путем преобразования в ee.Array , транспонирования его, а затем обратного преобразования в ee.List .

Редактор кода (JavaScript) 

// If variables in the list are arranged as rows, you'll need to transpose it.
// Define a list of lists where rows represent variables. The first row is the
// independent variable and the second is the dependent variable.
var listsVarRows = ee.List([
  [1, 2, 3, 4, 5],
  [1, 2, 3, 4, 5]
]);

// Cast the ee.List as an ee.Array, transpose it, and cast back to ee.List.
var listsVarColumns = ee.Array(listsVarRows).transpose().toList();

// Compute the least squares estimate of a linear function. Note that an
// object is returned; cast it as an ee.Dictionary to make accessing the
// coefficients easier.
var linearFit = ee.Dictionary(listsVarColumns.reduce(ee.Reducer.linearFit()));

// Inspect the result.
print(linearFit);
print('y-intercept:', linearFit.get('offset'));
print('Slope:', linearFit.get('scale'));

linearRegression()

Применение ee.Reducer.linearRegression() аналогично приведенному выше примеру LinearFit() , за исключением того, что включена постоянная независимая переменная.

Редактор кода (JavaScript) 

// Define a list of lists where columns represent variables. The first column
// represents a constant term, the second an independent variable, and the third
// a dependent variable.
var listsVarColumns = ee.List([
  [1, 1, 1],
  [1, 2, 2],
  [1, 3, 3],
  [1, 4, 4],
  [1, 5, 5]
]);

// Compute ordinary least squares regression coefficients. numX is 2 because
// there is one constant term and an additional independent variable. numY is 1
// because there is only a single dependent variable. Cast the resulting
// object to an ee.Dictionary for easy access to the properties.
var linearRegression = ee.Dictionary(
  listsVarColumns.reduce(ee.Reducer.linearRegression({
    numX: 2,
    numY: 1
})));

// Convert the coefficients array to a list.
var coefList = ee.Array(linearRegression.get('coefficients')).toList();

// Extract the y-intercept and slope.
var b0 = ee.List(coefList.get(0)).get(0); // y-intercept
var b1 = ee.List(coefList.get(1)).get(0); // slope

// Extract the residuals.
var residuals = ee.Array(linearRegression.get('residuals')).toList().get(0);

// Inspect the results.
print('OLS estimates', linearRegression);
print('y-intercept:', b0);
print('Slope:', b1);
print('Residuals:', residuals);