CC BY
45
УДК 528.854(630*2)
DOI: 10.37482/0536-1036-2022-2-84-104
повышение информативности мультиспектральных спутниковых изображений с использованием данных текстурного анализа
С.А. Зотов1, аспирант; ResearcherID: AAN-7199-2020, ORCID: https://orcid.org/0000-0002-3483-1449
Е.В. Дмитриев1,2, канд. физ.-мат. наук, ст. науч. сотр.; ResearcherID: E-4794-2014,
ORCID: https://orcid.org/0000-0001-5363-3934
П.Г. Мельник3,4, канд. с.-х. наук, доц.; ResearcherID: E-7644-2014,
ORCID: https://orcid.org/0000-0002-2802-7614
Т.В. Кондранин1, д-р физ.-мат. наук; ResearcherID: K-9407-2013,
ORCID: https://orcid.org/0000-0003-3565-3194
'Московский физико-технический институт (национальный исследовательский университет); Институтский пер., д. 9, г. Долгопрудный, Московская обл., Россия, 141701; e-mail: [email protected], [email protected]
2Институт вычислительной математики РАН, ул. Губкина, д. 8, Москва, Россия, 119333; e-mail: [email protected]
3Московский государственный технический университет им. Н.Э. Баумана (Мытищинский филиал), 1-я Институтская ул., д. 1, г. Мытищи, Московская обл., Россия, 141005; e-mail: [email protected]
4Институт лесоведения РАН, ул. Советская, д. 21, с. Успенское, Московская обл., Россия, 143030; e-mail: [email protected]
оригинальная статья / Поступила в редакцию 26.03.20 / Принята к печати 03.08.20
Аннотация. Рассматривается задача повышения информативности тематической обработки мультиспектральных изображений среднего (10-30 м) и высокого (1-4 м) пространственного разрешения, получаемых зарубежными и отечественными спутниковыми системами дистанционного зондирования Земли, за счет привлечения дополнительной текстурной информации с панхроматических спутниковых изображений сверхвысокого пространственного разрешения (^(1-0,4) м). Объектом исследования стали изображения тестовых полигонов на территории Савватьевского лесничества (Тверская обл.) со спутников Landsat 8, Sentinel 2 и WorldView 2, оснащенных мультиспектраль-ной аппаратурой. Для валидации результатов расчетов использованы геопривязанные данные наземных обследований. В качестве спектральных признаков взяты значения спектральной отражательной способности в каналах видимого и ближнего инфракрасного диапазонов, нормализованные на соответствующую интегральную характеристику. в целях извлечения текстурных признаков на основе распределения совместной встречаемости уровней серого (текстурные признаки Харалика) в пределах скользящего окна, пробегающего с заданным пространственным шагом изображение, вычисляли статистические характеристики. Проведен корреляционный анализ текстурных признаков с учетом изменений расстояния и угла смежности. Для выбранных определяющих признаков: автокорреляции, асимметрии, контраста и корреляции - показано, что первые три могут быть использованы с произвольным выбранным направлением смежности, а последний необходимо рассматривать в двух различных направлениях. Установ-
Данная статья опубликована в режиме открытого доступа и распространяется на условиях лицензии Creative Commons «Attribution» («Атрибуция») 4.0 Всемирная (CC BY 4.0) • Авторы заявляют об отсутствии конфликта интересов
лено, что при решении задачи распознавания природных и антропогенных объектов все рассмотренные алгоритмы классификации обеспечивают существенное повышение точности при совместном использовании спектральных и текстурных признаков по сравнению с традиционной спектральной классификацией. Этот результат продемонстрирован для всех полученных различными спутниками изображений тестовых полигонов. Можно сделать предварительный вывод, что предлагаемый комплексный подход тематической обработки позволяет повысить качество распознавания объектов в случае использования изображений как среднего, так и высокого пространственного разрешения. Оценки, полученные при построении тематических карт доминантных и субдоминантных отделов лесной растительности, показали близкие точности классификации по различным исходным мультиспектральным изображениям (с разбросом не более 5 % около среднего значения 85 %). По большей части это связано с наличием специфических погрешностей результатов наземной лесной таксации и говорит о необходимости их актуализации с применением данных дистанционного спутникового зондирования.
Для цитирования: Зотов С.А., Дмитриев Е.В., Мельник П.Г., Кондранин Т.В. Повышение информативности мультиспектральных спутниковых изображений с использованием данных текстурного анализа // Изв. вузов. Лесн. журн. 2022. № 2. С. 84-104. DOI: 10.37482/0536-1036-2022-2-84-104
Финансирование: Работа выполнена при финансовой поддержке РФФИ, проект № 20-07-00370 «Фундаментальные проблемы повышения информативности обработки данных оптоэлектронных авиакосмических устройств высокого пространственного и спектрального разрешения».
Ключевые слова: распознавание образов, информативность, дистанционное зондирование, мультиспектральные изображения, почвенно-растительный покров, тематическая обработка, лесная таксация.
INCREASING THE INFORMATIVITY OF MULTISPECTRAL SATELLITE IMAGES USING TEXTURE ANALYSIS DATA
Sergey A. Zotov1, Postgraduate Student; ResearcherID: AAN-7199-2020, ORCID: https://orcid.org/0000-0002-3483-1449
Egor V. Dmitriev12, Candidate of Physics and Mathematics, Senior Research Scientist; ResearcherID: E-4794-2014, ORCID: https://orcid.org/0000-0001-5363-3934 Petr G. Melnik3,4, Candidate of Agriculture, Assoc. Prof.; ResearcherID: E-7644-2014, ORCID: https://orcid.org/0000-0002-2802-7614
Timofey V. Kondranin1, Doctor of Physics and Mathematics; ResearcherID: K-9407-2013, ORCID: https://orcid.org/0000-0003-3565-3194
'Moscow Institute of Physics and Technology, Institutskiy per., 9, Dolgoprudny, Moscow Region, 141701, Russian Federation; e-mail: [email protected], [email protected] 2Marchuk Institute of Numerical Mathematics RAS, ul. Gubkina, 8, Moscow, 119333, Russian Federation; e-mail: [email protected]
3Mytischi Branch of Bauman Moscow State Technical University, ul. 1-ya Institutskaya, 1, Mytischi, Moscow Region, 141005, Russian Federation; e-mail: [email protected] 4Institute of Fоrеst Science, Russian Academy of Sciences, ul. Sovetskaya, 21, s. Uspenskoe, Moscow Region, 143030, Russian Federation; e-mail: [email protected]
Original article /Received on March 26, 2020 /Accepted on August 3, 2020
This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC BY 4.0) license • The authors declare that there is no conflict of interest
Abstract. The article considers the problem of incresing the informativity of multispectral images of medium (10-30 m) and high (1-4 m) spatial resolution obtained from foreign and national satellite remote sensing systems by involving additional textural information from panchromatic satellite images of very high spatial resolution (<(1-0.4) m). The images of test sites on the territory of Savvat'yevo forestry (Tver region) from Landsat 8, Sentinel 2 and WorldView 2 satellites equipped with multispectral instruments were an object of this research. Geo-referenced ground survey data were used to validate the calculation results. We used the values of the spectral reflectance in the visible and near-infrared channels normalized to the appropriate integral characteristic as spectral features. Statistical characteristics were calculated in order to extract texture features based on the distribution of the co-occurrence of gray levels (Haralick texture features) within a moving window running the image with a given spatial step. A correlation analysis of textural features was carried out considering changes in distance and angle of adjacency. It was shown that for the selected leading features (autocorrelation, asymmetry, contrast and correlation) the first three can be used with an arbitrary direction of adjacency, while the latter needs to be considered in two different directions. Also we have found that all the considered classification algorithms provide a significant increase of accuracy when both spectral and textural features are used, in comparison with the traditional spectral classification. This result was shown for all images of test sites obtained by different satellites. It is possible to make a preliminary conclusion that the proposed integrated approach of thematic processing can improve the quality of object recognition in the case of using images of both medium and high spatial resolution. Estimates obtained during the thematic mapping of dominant and subdominant forest species showed close classification accuracies for different initial multispectral images (with a scatter of no more than 5 % around the average value of 85 %). Mostly this is due to the presence of specific errors in the ground-based forest inventory data and indicates the necessity of their updating with the use of satellite remote sensing images.
For citation: Zotov S.A., Dmitriev E.V., Melnik P.G., Kondranin T.V. Increasing the Informativity of Multispectral Satellite Images Using Texture Analysis Data. Lesnoy Zhurnal [Russian Forestry Journal], 2022, no. 2, pp. 84-104. DOI: 10.17238/issn0536-1036-2022-2-84-104
Funding: The work was carried out with the financial support of the Russian Foundation for Basic Research; project No. 20-07-00370 "Fundamental Problems of Improving the Informativity of Data Processing of Optoelectronic Aerospace Devices with High Spatial and Spectral Resolution".
Keywords: pattern recognition, informativity, remote sensing, multispectral images, soil and vegetation cover, thematic processing, forest inventory.
Введение
В настоящее время мультиспектральные спутниковые изображения среднего (10-30 м) и высокого (1-4 м) пространственного разрешения (ПР) широко используются для дистанционного мониторинга почвенно-раститель-ного покрова [14]. Преимуществами изображений среднего разрешения являются возможность регулярно получать информацию о состоянии исследуемых объектов на выбранной территории и большая скорость тематической обработки [19]. Спутниковые изображения высокого разрешения обеспечивают высокую детализацию объектов, однако при практическом использовании таких данных возникают связанные с их значительной коммерческой стоимостью трудности.
При оценке состояния растительности широкое распространение получили технологии обработки мультиспектральных изображений в каналах видимого и ближнего инфракрасного диапазонов с применением концепции вегетационных индексов [1, 5, 9]. Многочисленные исследования показали, что одной из наиболее сложных проблем при использовании данной концепции становится выявление структурных особенностей и динамики изменений лесного полога. В частности, в работе [9] предложен статистический метод восстановления сомкнутости полога и индекса листовой поверхности хвойных лесов центральной части Канады с использованием относительного и нормализованного разностного вегетационных индексов. Точность дистанционного определения сомкнутости полога и индекса листовой поверхности, как было доказано, существенно варьирует в зависимости от периода вегетации и условий наблюдения (углов визирования), вызывая трудности интерпретации результатов тематической обработки.
Альтернативный подход основан на решении обратной задачи переноса излучения в системе «лесной полог - атмосфера». Одной из первых, в которой продемонстрирован данный подход, была работа [12], где индекс листовой поверхности восстанавливался по мультиспектральным изображениям со спутника Landsat ETM+ (ПР - 30 м) на основе прямого моделирования отражательной способности лесного полога в спектральных каналах прибора ЕТМ.
В последующие годы с использованием мультиспектральных данных Landsat TM и Landsat ETM+ был разработан стандартный продукт Landsat TCC [21], представляющий собой глобальную карту распределения площади проективного покрытия лесной растительности. Валидация Landsat TCC с применением данных наземных обследований полигонов на территории Финляндии [16] показала достаточно приемлемые с потребительской точки зрения точности определения проективного покрытия леса (« 85 %). Важно подчеркнуть, что при оценке качества определения структурных параметров древостоев, в том числе и проективного покрытия, необходимо учитывать возможность ошибки в данных наземных обследований [20].
В последние годы в эксплуатацию введены современные спутниковые системы, такие как WorldView, Quickbird, Pléiades. На спутниках этих серий установлены мультиспектральная (8 (WorldView) и 4 (Quickbird и Pléiades) канала диапазоном 400-1000 нм, ПР - 2-3 м в надире) и панхроматическая аппаратура (ПР ~ 0,5-0,7 м), система Sentinel (12 каналов диапазоном 400-2500 нм, ПР - 10, 20 и 60 м), а также система Cartosat, имеющая только панхроматическую стереокамеру (ПР - 2,7 м). открываются новые возможности совместного использования мультиспектральных и панхроматических данных, получаемых синхронно в режиме реального времени, причем подчеркнем, что ПР муль-тиспектральных изображений систем WorldView, Quickbird и Pléiades в разы выше ПР классических систем типа Landsat и современных Sentinel.
К настоящему времени выполнен ряд исследований, посвященных использованию высокодетальной текстурной информации при тематической обработке мультиспектральных изображений; рассмотрены различные схемы совместной спектрально-текстурной обработки дистанционных данных применительно к оценке состояния растительного покрова [6, 11, 13]. Показана эффективность использования дополнительной текстурной информации
для повышения точности классификации характеристик лесной растительности в различных климатических условиях - во влажных тропических лесах [13] и смешанных лесах центральных регионов России [11]. Важно отметить, что получаемые при таком подходе более надежные оперативные данные о ле-сотаксационных характеристиках (размерах крон, высоте и диаметре деревьев) используются в дальнейшем для решения важных прикладных задач лесопользования: определения биопродукционных параметров лесных территорий на большой по размеру площади.
Цель работы - исследование возможностей повышения точности тематической обработки мультиспектральных спутниковых изображений среднего и высокого пространственного разрешения совместно с высокодетальной текстурной информацией, получаемой современными оперативными спутниковыми системами в виде координатно-привязанных панхроматических изображений сверхвысокого пространственного разрешения. Также ставится задача сравнения эффективности применения мультиспектральных данных среднего и высокого пространственного разрешения для мониторинга почвенно-растительного покрова при использовании различных методов автоматизированной обработки.
Объекты и методы исследования
Мультиспектральные спутниковые системы дистанционного зондирования в заданном для конкретной аппаратуры наборе узких спектральных каналов обеспечивают получение изображений наблюдаемых объектов, представляемых в виде так называемого спектрального куба данных с двумя пространственными и одной спектральной координатой.
В качестве признаков спектральной координаты в данной работе используется набор спектральных яркостей пикселей в каналах конкретной спутниковой аппаратуры. Для каждого пикселя производится нормализация спектральных яркостей на интегральную яркость в канале, что позволяет существенно сгладить разброс данных, связанный, например, с изменением условий освещения.
При вычислении текстурных признаков применяется статистический подход [4], согласно которому описание текстуры изображения основано на вычислении двумерной функции вероятности совместной встречаемости уровней серого в заданном направлении [17]:
ОЬСМ(г,»
P(i. j) =-
X GLCM(i, j)
i.j=
Элементами матрицы GLCM (i, j) (англ. Gray-Level Co-occurrence Matrix) являются частоты появления пары пикселей с яркостями i и j (N - число уровней серого). Пиксели расположены на расстоянии d (в пикселях) друг от друга и образуют угол смежности а с заданным референсным направлением.
Матрица GLCM рассчитывается для частей обрабатываемого панхроматического изображения, попадающих внутрь скользящего окна. Центр окна пробегает пиксели основного панхроматического изображения. Ширина окна подбирается в зависимости от характерных размеров обследуемых объектов и
пространственного разрешения дистанционных данных. При обработке изображений лесных территорий ряд авторов рекомендуют использование скользящего окна шириной от 20 до 80 пикс. [2, 3].
На основе p (i, j) рассчитываются статистики, известные как текстурные признаки Харалика [18]. Нами рассматривается 19 текстурных признаков Харалика, наиболее часто используемых при текстурном анализе изображений [8, 17]: автокорреляция, островершинность, асимметрия, контраст, корреляция, энтропия разности, дисперсия разности, неоднородность, энергия, энтропия, однородность, однородность 2, информационная мера корреляции 1, информационная мера корреляции 2, максимальная вероятность, среднее суммы, энтропия суммы, сумма квадратов, дисперсия суммы.
Для эффективного применения текстурных признаков Харалика при классификации объектов подстилающей поверхности требуется достаточно высокое ПР обрабатываемых изображений, т. е. размер пикселя должен быть существенно меньше характерного размера текстурных элементов обследуемых объектов.
Мультиспектральные изображения, получаемые с современных спутниковых систем, как правило, имеют более низкое ПР, чем панхроматические. Однако применение более крупного шага движения панхроматического окна позволяет повысить эффективность совместной обработки обоих изображений. В алгоритмах, реализованных в настоящей работе, центр скользящего окна пробегает каждый пиксель многоспектрального изображения с одновременным «вырезанием» части соответствующего панхроматического изображения. В результате изображения текстурных признаков добавляются к спектральным признакам соответствующего пикселя изображения из центра скользящего окна. Далее классификация проводится по спектральным и текстурным признакам совместно [6].
В основе алгоритмов классификации лежат методы машинного обучения с учителем, суть которых в том, что определенные числовые параметры алгоритма оцениваются по так называемой обучающей выборке - множеству из m пар признаков х е X и соответствующих им наименований y е Y.
В качестве классификаторов рассмотрены следующие алгоритмы: k-ближайших соседей (k-nearest neighbors (kNN)) [10, 15], «случайный лес» (random forest (RF)) [7], квадратичного дискриминантного анализа (quadratic discriminant analysis (QDA) [23].
В алгоритме kNN объект входит в класс, к которому принадлежит большинство из k-ближайших соседей обучающей выборки. В формализованном виде это условие записывается так:
m
a(x) = arg max y(i) = y] w(i, x),
где a(x) е Y - номер класса, определяемого по вектору признаков; [y^ = y] - индикаторная функция; w (i, x) - весовая функция i-го элемента обучающей выборки по отношению к классифицируемому вектору признаков x. В нашем случае после перегруппировки обучающей выборки в порядке увеличения расстояния от х весовая функция представима в виде w (i, x) = [i < k].
Необходимым элементом kNN-алгоритма является определение «расстояния» между объектами, для чего в общем случае в пространстве признаков объектов вводится метрика из набора существующих, например Евклидова, Манхэттенская или Чебышевская. С учетом специфики координатной привязки
изображений в работе используется Евклидова метрика, весовые коэффициенты принимаются одинаковыми для всех соседей.
Не менее важной проблемой kNN-алгоритма является правильный выбор количества ближайших соседей k, которое подбиралось нами таким образом, чтобы при обработке эталонных объектов точность классификации была наибольшей. kNN-алгоритм прост для реализации, но в то же время для обеспечения приемлемых точностей классификации необходимы повышенная репрезентативность обучающей выборки и, как следствие, высокие вычислительные затраты.
Алгоритм «случайный лес» основан на построении «леса» из N «решающих деревьев» заданной глубины, каждое из которых обучается на случайной подгруппе обучающей выборки [7]. Композиция из N-алгоритмов bx(x), ..., bN(x) является, по существу, регрессионным алгоритмом, усредняющим результаты обработки каждого из этих алгоритмов. Композицию алгоритмов ЬДх), ..., bN(x) можно представить в виде
1 N
a( х) = sign — X bn (х)
N n=1
Количество деревьев в лесу и их глубина подбираются таким образом, чтобы точность классификации, отрабатываемой на эталонных объектах, была наибольшей. RF-алгоритм имеет относительно малую склонность к переобучению, в минимальной степени подвержен влиянию шумов и выбросов. К его недостаткам, помимо указанных для kNN-алгоритма, следует отнести высокую сложность реализации.
Алгоритм QDA - это частный случай байесовского классификатора. Формализм QDA в общем случае описывается в виде
a( x) = argmax Pypy (x).
yeY
Здесь Py - это априорная вероятность класса y, т. е. вероятность события появления (существования) этого класса вне зависимости от причин, формирующих реализацию x. Распределение признаков внутри класса y характеризуется функцией условной плотности распределения вероятности py(x).
В случае, если признаки каждого класса описываются нормальными распределениями, что имеет место при решении многих прикладных задач,
a(х) = arg max[2In Ру - (х - цу f I-1 (х - цу)- In (det(Sу))], байесовский классификатор сводится к формуле
где ц - математическое ожидание вектора признаков; Ey - ковариационная матрица признаков, определяемых для класса y [22].
По сравнению с алгоритмами kNN и RF алгоритм QDA требует меньших вычислительных затрат.
Для отработки методики использовали мультиспектральные изображения среднего и высокого ПР и панхроматические сверхвысокого ПР территорий Сав-ватьевского лесничества (Тверская обл.), для которых есть данные наземной лесной таксации, необходимые для валидации результатов тематической обработки спутниковых изображений.
На рис. 1 представлено изображение полигона с таксационными выделами и контурами, используемыми в численных экспериментах по совместной обработке мультиспектральных и панхроматических изображений.
Рис. 1. Тестовый участок на территории Савватьевского лесничества (Тверская обл.): желтые линии - границы таксационных выделов, синяя и фиолетовая рамки - контуры используемых изображений WorldView 2
Fig. 1. Test site on the territory of Savvat'yevo forestry (Tver region): yellow lines -boundaries of inventory plots; blue and purple frames - bounding boxes of the WorldView 2
images used
Согласно данным лесной таксации на территории Савватьевского лесничества преобладают две породы: береза и сосна. Примесь остальных пород незначительна, что затрудняет их представление в обучающей выборке. Поэтому деревья классифицировали на лиственные и хвойные породы, представленные в подавляющем большинстве березой и сосной. Рассматриваемый тестовый участок содержал 604 выдела, в 460 из которых доминировали хвойные породы, а в оставшихся 144 - лиственные.
В качестве источника спектральных данных о поверхности выбранной территории использованы мультиспектральные снимки трех разных спутников: Landat 8, Sentinel 2, WorldView 2.
В табл. 1 приведены каналы, их номера в списке каналов соответствующего спутника (N), спектральные диапазоны (ДА,) и пространственные разрешения (Ах). Все эти значения использовались далее в качестве информативных признаков для классификации.
Для вычисления текстурных признаков в рамках численного эксперимента взяты два снимка от 25.06.2016 с разрешением 0,46 м, сделанных панхроматической камерой WorldView 2.
Таблица 1
каналы мультиспектральных спутниковых систем, выбранные в качестве информативных признаков
Канал Landsat 8 Sentinel 2 WorldView 2
N ДХ, мкм Дх, м N ДХ, мкм Дх, м N ДХ, мкм Дх, м
Blue 2 0,4520,512 30 2 0,4580,523 10 2 0,4500,510 1,84
Green 3 0,5330,590 30 3 0,5430,578 10 3 0,5100,580 1,84
Red 4 0,6360,673 30 4 0,6500,680 10 5 0,6300,690 1,84
NIR-1 5 0,8510,879 30 8 0,7880,903 10 7 0,7700,895 1,84
NIR-2 - - - - - - 8 0,8600,940 1,84
SWIR-1 6 1,5661,651 30 11 1,5651,655 20 - - -
SWIR-2 7 2,1072,294 30 12 2,1002,280 20 - - -
Результаты исследования и их обсуждение
Комбинация мультиспектральных и панхроматических снимков полностью отражает полигон. Перед вычислением текстурных признаков яркости обрабатываемого панхроматического изображения градуируются таким образом, чтобы минимальное значение яркости было 0, а максимальное - 100. Тестирование показало, что дальнейшее увеличение числа градаций яркости существенно не влияет на текстурные признаки.
Каждый из 19 текстурных признаков возможно вычислить при разных расстояниях d и углах а. Такое количество признаков при классификации бывает избыточным, поэтому проводится отбор информативных признаков. Текстурные признаки проверяются на информативность: вычисляются их значения при а = 0°, d = 1 по объектам, входящим в сбалансированную обучающую выборку, которая включает в себя объекты, характерные для выбранной местности: лиственный лес, хвойный лес, луговая растительность, вода, песок, массивы жилых домов. Термин «сбалансированная выборка» означает равенство количества объектов каждого класса. Также отметим, что в ходе тестирования различных моделей был установлен наилучший результат при ширине скользящего окна 80 пикс., поэтому данная ширина использовалась для вычисления текстурных признаков.
Для определения информативности каждого признака по всем статистикам вычисляется корреляционная матрица Спирмена, представленная на рис. 2.
Автокорреляция Среднее суммы Островершинность Корреляция Энергия Энтропия
Информационная мера корреляции 1 Информационная мера корреляции 2 Максимальная вероятность Энтропия суммы Сумма квадратов Дисперсия суммы Асимметрия Контраст Энтропия разности! Дисперсия разности Неоднородность Однородность Однородность 2
к
Ml
Рис. 2. Корреляционная матрица Спирмена в абсолютных значениях, полученная по статистикам текстурных признаков с углом а = 0° и расстоянием d = 1. Текстурные признаки вычислены по объектам сбалансированной обучающей выборки
Fig. 2. Absolute Spearman correlation matrix derived from the texture feature statistics with angle а = 0° and distance d = 1. The texture features were calculated for the objects of the
balanced training set
Если коэффициент корреляции Спирмена двух признаков близок к +1 или —1, то эти признаки имеют линейную зависимость и только один из них является информативным. Корреляция Спирмена была выбрана, потому что она более устойчива к выбросам. Также корреляция Спирмена в отличие от корреляции Пирсона не меняется значимо при различных сбалансированных обучающих выборках.
Считаем, что признаки линейно зависимы при корреляции Спирмена между наборами их значений более 0,9 в абсолютном значении. Тогда при вычислении значений признаков по сбалансированной обучающей выборке можно выделить следующие группы линейно зависимых признаков (в том числе группы, состоящие из одного признака): автокорреляцию, среднее суммы; островершинность, корреляцию, энергию, энтропию, информационную меру корреляции 1, информационную меру корреляции 2, максимальную вероятность, энтропию суммы, сумму квадратов, дисперсию суммы; асимметрию;
контраст, энтропию разности, дисперсии разности, неоднородность, однородность, однородность 2.
для проведения спектрально-текстурной тематической обработки из каждой группы выбирается по одному признаку. По нашим предварительным оценкам, изменение признаков, выбираемых из выделенных групп, не влечет за собой значимых изменений результатов расчетов. Таким образом, выбор из группы произволен. При отборе наиболее информативных признаков будут участвовать следующие текстурные признаки Харалика: автокорреляция, контраст, асимметрия, корреляция.
Эти текстурные признаки будут называться базовыми. Для оценки влияния параметра а на значение каждого из них вычисляется корреляционная матрица Спирмена по статистикам с разным углом а (0, 45, 90 и 135°), расстоянием d, равным 1. Получившиеся матрицы представлены на рис. 3, а.
Угол а, по нашему мнению, существенно влияет на текстурный признак, вычисляемый при коэффициентах корреляции Спирмена больше 0,9 в абсолютном значении.
Анализ корреляционных матриц показывает: автокорреляция не зависит от значения а; контраст и асимметрия не зависят существенно от значения а; корреляция при а = 135° существенно отличается от значений при остальных а. При этом в среднем корреляция этого признака при а = 135° одинакова с корреляцией при других а (значения корреляционных матриц могут колебаться, т. к. балансировка тестовой выборки осуществляется случайным образом).
Поэтому в отборе признаков для классификации будут участвовать отдельно значения текстурного признака «корреляция», вычисленные при а = 0° и а = 135° (далее помечено как а = 0° и а = 135° соответственно).
Аналогичный анализ значимости параметра d при вычислении базовых признаков позволяет обнаружить, что базовые признаки существенно не зависят от d, если оно равно 1, 2, 3. Корреляционные матрицы Спирмена, вычисленные по статистикам базовых текстурных признаков с расстояниями d от 1 до 3, представлены на рис. 3, б. По этим матрицам был подтвержден вывод об отсутствии существенной зависимости базовых текстурных признаков от d.
Итак, в отборе признаков для классификации будут участвовать следующие текстурные признаки харалика: автокорреляция, контраст, асимметрия, корреляция (а = 0°), корреляция (а = 135°).
Далее для каждого набора данных и метода классификации производится поэтапный отбор признаков. На первом этапе методом полного перебора подбираются спектральные признаки по максимальной точности на некоторой тестовой выборке. на втором этапе, если в рамках тестируемой модели спектральные признаки используются с текстурными, то аналогично подбирается набор базовых текстурных признаков, дополняющих отобранные спектральные признаки.
Все тестируемые в рамках численного эксперимента модели различаются по методам классификации / алгоритмам (3 варианта), используемым наборам спектральных данных (3 варианта), использованию (или неиспользованию) текстурных данных совместно со спектральными (2 варианта).
Автокорреля
Контраст
б
Рис. 3. Корреляционные матрицы Спирмена в абсолютных значениях, полученные по статистикам базовых текстурных признаков: а - с разными углами смежности а и расстоянием смежности d = 1; б - при угле а = 0 и разных расстояниях d. Текстурные признаки
вычислялись по объектам сбалансированной обучающей выборки Fig. 3. Absolute Spearman correlation matrix derived from the statistics of the basic texture features: а - with different angles а and distance d = 1; б - with the angle of adjacency а = 0 and different adjacency distances d. The texture features were calculated for the objects of the
balanced training set
Варианты наборов спектральных данных в соответствии с табл. 1: муль-тиспектральные изображения Landsat 8, мультиспектральные изображения Sentinel 2, мультиспектральные изображения WorldView 2.
Итого в рамках численного эксперимента тестируются 18 моделей.
а
Для каждого метода классификации путем полного перебора производится отбор признаков по точности на отдельной тестовой выборке: спектральных -для моделей, использующих только спектральные признаки, а также спектральных и текстурных - для моделей, использующих оба вида признаков. Метод полного перебора предполагает выбор набора признаков из всех возможных комбинаций по точности предсказания на тестовой выборке. Результаты представлены в табл. 2. Эти наборы данных не являются единственно верными и могут отличаться при существенных изменениях тестовых выборок.
Таблица 2
признаки, отобранные для каждой тестируемой модели
Метод классификации Спутник Признаки Отобранные признаки
RF Landsat 8 Спектральные SWIR-1, SWIR-2
Спектральные и текстурные SWIR-1, SWIR-2, автокорреляция, асимметрия, контраст, корреляция (а = 0°)
Sentinel 2 Спектральные Green, SWIR-1, SWIR-2
Спектральные и текстурные Green, SWIR-1, SWIR-2, автокорреляция, асимметрия, контраст, корреляция (а = 0°), корреляция (а = 135°)
Worldview 2 Спектральные Blue, Green, Red, NIR-1, NIR-2
Спектральные и текстурные Blue, Green, Red, NIR-1, NIR-2, автокорреляция, контраст, асимметрия, корреляция (а = 0°), корреляция (а = 135°)
kNN Landsat 8 Спектральные SWIR-1, SWIR-2
Спектральные и текстурные SWIR-1, SWIR-2, автокорреляция, асимметрия
Sentinel 2 Спектральные Blue, SWIR-1, SWIR-2
Спектральные и текстурные Blue, SWIR-1, SWIR-2, асимметрия, контраст
Worldview 2 Спектральные Blue, Green, Red, NIR-1, NIR-2
Спектральные и текстурные Blue, Green, Red, NIR-1, NIR-2, автокорреляция, контраст, асимметрия, корреляция (а = 135°)
QDA Landsat 8 Спектральные Green, Red, SWIR-1, SWIR-2
Спектральные и текстурные Green, Red, SWIR-1, SWIR-2, автокорреляция, контраст, корреляция (а = 0°), корреляция (а = 135°)
Sentinel 2 Спектральные Red, SWIR-1, SWIR-2
Спектральные и текстурные Red, SWIR-1, SWIR-2, автокорреляция -0,0, контраст, корреляция (а = 0°)
Worldview 2 Спектральные Green, Red, NIR-1, NIR-2
Спектральные и текстурные Green, Red, NIR-1, NIR-2, автокорреляция, контраст
Порядок численного эксперимента, применимый для каждой модели из табл. 2, следующий. На первом этапе происходит обучение выбранного классификатора по группированному набору спектральных или спектрально-текстурных признаков. Далее проводится попиксельная классификация. Для каждого таксационного выдела определяется преобладание в его границах пикселей, классифицированных как лиственный или хвойный лес, т. е. выявляется доминантный класс. Если в контур выдела не попадает ни один пиксель, определенный как лиственный или хвойный лес, считается, что доминантный класс в этом выделе отсутствует, а если в границы выдела не попадает центр ни одного пикселя, то - что доминантный класс определить не удается. Также для каждого выдела устанавливается субдоминантный класс - второй класс по количеству пикселей в контуре выдела из классов «лиственный лес» или «хвойный лес». Если в границах контура нет пикселей недоминантного класса, то субдоминантный класс для данного выдела отсутствует.
Точности моделей определяются долей выделов, для которых доминантный (или субдоминантный) класс вычислен верно, исключая выделы, в которых нет доминантного (или субдоминантного) класса, и выделы, для которых доминантный (или субдоминантный) класс не был определен. Результаты тестирования моделей, в том числе точности построения тематических карт доминантных и субдоминантных отделов, а также усредненная точность распознавания всего рассматриваемого спектра объектов, устанавливаемая методом 10-блочной кросс-валидации, сведены в табл. 3.
Метод 10-блочной кросс-валидации показывает, что увеличение точности при переходе от спектральной к спектрально-текстурной классификации существенно превосходит различие точности рассматриваемых методов.
Дополнительно мы провели визуальный анализ качества результатов попиксельной спектрально-текстурной классификации изображения всей обследуемой территории. В качестве примера на рис. 4 представлены результаты классификации методом kNN, полученные на основе мультиспектральных изображений Landsat 8, Sentinel 2 и WorldView 2, а также текстурное изображение высокого разрешения для сравнения. Видно, что при использовании данных Landsat 8, имеющих наиболее низкое пространственное разрешение, в большей степени нарушается геометрическая структура восстанавливаемых объектов, в особенности это касается зон застройки. Данные тематической обработки изображений WorldView 2 имеют наименьшее количество артефактов.
Сделанные выводы согласуются с результатами кросс-валидации, которые представлены в табл. 3 для метода kNN. Спектрально-текстурная классификация изображений Landsat 8 и Sentinel 2 имеет значимо более низкую точность по сравнению с классификацией изображений WoldView 2. Данный эффект отсутствует при использовании метода QDA. Часть артефактов, возникающих при обработке данных Sentinel 2 и WorldView 2, объясняется наличием дорожной сети, которая в большинстве случаев дешифрования может быть классифицирована как зона застройки.
Таблица 3
результаты классификации наборов данных разными методами
Точность
Метод спутник Признаки по доминантам по субдоминантам по кросс-валидации
Landsat 8 Спектральные 0,834 0,709 0,732
Спектральные и текстурные 0,853 0,762 0,852
RF Sentinel 2 Спектральные 0,861 0,763 0,788
Спектральные и текстурные 0,836 0,719 0,930
WorldView 2 Спектральные 0,852 0,767 0,624
Спектральные и текстурные 0,857 0,770 0,924
Landsat 8 Спектральные 0,840 0,708 0,774
Спектральные и текстурные 0,842 0,728 0,777
Sentinel 2 Спектральные 0,867 0,767 0,793
kNN Спектральные и текстурные 0,851 0,741 0,794
WorldView 2 Спектральные 0,876 0,793 0,666
Спектральные и текстурные 0,882 0,815 0,928
Landsat 8 Спектральные 0,835 0,708 0,778
Спектральные и текстурные 0,847 0,744 0,890
QDA Sentinel 2 Спектральные 0,880 0,801 0,800
Спектральные и текстурные 0,872 0,777 0,859
WorldView 2 Спектральные 0,885 0,800 0,602
Спектральные и текстурные 0,865 0,781 0,875
Рис. 4. Изображения, полученные в результате попик-сельной классификации методом kNN по спектральным и текстурным признакам, сделанные на основе мультиспек-тральных изображений Landsat 8 (а), Sentinel 2 (б) и WorldView 2 (в), и панхроматическое изображение высокого разрешения (г)
Fig. 4. Results of pixel-by-pixel classification by kNN algorithm with the use of both spectral and texture features extracted from multispectral images of Landsat 8 (а), Sentinel 2 (б) and WorldView 2 (в), and a high-resolution panchromatic image (г)
На рис. 5 показаны доминантные классы в контурах выделов. Тематические карты были построены по результатам спектрально-текстурной классификации методом kNN мультиспектральных изображений, полученных со спутников Landsat 8, Sentinel 2, WorldView 2, а также по данным наземной лесной таксации. Контуры выделов на рисунке залиты цветами, обозначающими доминантные классы. Аналогичным образом показаны тематические карты субдоминантных отделов (рис. 6). Следует отметить, что спутниковые тематические карты, в особенности построенные с использованием изображений Sentinel 2 и WorldView 2, имеют большую схожесть.
Рис. 5. Доминантные классы в контурах выделов лесной таксации Савватьевского лесничества, определенные по классификации методом kNN по спектральным и текстурным признакам, сделанные на основе мультиспектральных изображений Landsat 8 (а), Sentinel 2 (б), WorldView 2 (в), а также по наземным данным лесной таксации (г). Синяя и фиолетовая рамки - контуры используемых изображений WorldView 2
Fig. 5. Dominant classes within the boundaries of forest inventory plots of Savvat'yevo forestry calculated from classification by kNN algorithm with the use of spectral and texture features extracted from multispectral images of Landsat 8 (a), Sentinel 2 (б), WorldView 2 (в), as well as obtained from ground-based forest inventory data (г). Blue and purple frames are bounding boxes of the WorldView 2 images used
Согласно табл. 3 точность определения доминантных пород выделов Савватьевского лесничества с помощью всех рассматриваемых моделей составляет приблизительно 83.. .89 %. При этом нет выраженной разницы между результатами моделей, использующих только спектральные признаки при классификации и спектральные совместно с текстурными. Добавление текстурных признаков, как видно из таблицы, увеличивает качество кросс-валидации на обучающей выборке, особенно в случаях с моделями, использующими мультиспектраль-ные фотографии WorldView 2. Это связано с тем, что применение текстурных признаков позволяет яснее отличать луговую растительность и другие типы
подстилающей поверхности от лесной, но при распознании лиственной и хвойной растительности ключевыми являются именно спектральные признаки.
Также по табл. 3 можно проследить зависимость точности определения доминантных классов от использованных мультиспектральных снимков, которые, в свою очередь, отличаются друг от друга ПР. Точность моделей на основе изображений Landsat 8 с ПР 30 м меньше, чем точность моделей изображений Sentinel 2 с ПР 10 м, в среднем на 2 %. При этом модели, работающие с изображением Sentinel 2, показали точности определения доминантных классов хуже на < 1 % по сравнению с моделями, работающими с изображением WorldView 2 (ПР ~ 2 м). Метод RF дает наименьшую точность определения доминантных классов, а kNN и QDA отличаются незначительно.
Точность классификации по доминантам не превышает 89 %. Вероятно, это связано не только с ошибками моделей классификации, но и с устареванием данных наземной лесной таксации, ошибками в них. Некоторые модели, как показано на рис. 5, определяют, что для части выделов доминанты нет, т. е. в этих выделах леса нет вообще (хотя по данным лесной таксации он должен быть), и такие выделы повторяются в результатах классификации разными методами. Выделов без леса 17. Тщательное рассмотрение позволяет установить выделы, которые могут давать ошибку при классификации: в 11 выделах лес занимает менее 60 % площади, в 7 - полнота очень низкая, в 27 - были посажены лесные культуры, в 26 - предположительно смешанный лес, поэтому точное определение класса таких выделов затруднено.
Рис. 6. Субдоминантные классы в контурах выделов лесной таксации Савватьевского лесничества, определенные по классификации методом kNN по спектральным и текстурным признакам, сделанные на основе мультиспектральных изображений Landsat 8 (а), Sentinel 2 (б), WorldView 2 (в), а также по наземным данным лесной таксации (г). Синяя и фиолетовая рамки - контуры используемых изображений WorldView 2
Fig. 6. Subdominant classes within the boundaries of the inventory plots of Savvat'yevo forestry calculated from by classification by kNN algorithm with the use of spectral and textural features extracted from multispectral images of Landsat 8 (a), Sentinel 2 (б), WorldView 2 (в), as well as obtained from ground-based forest inventory data (г). Blue and purple frames are bounding boxes of the WorldView 2 images used
Согласно рис. 6, модели классификации, работающие с данными меньшего пространственного разрешения, определяют большее количество выделов как выделы без субдоминантной породы, что во многих случаях соотносится с данными наземных исследований. Это связано с тем, что попиксельная классификация по данным высокого разрешения с большей вероятностью определит хотя бы один пиксель из выдела как другую породу, которой там быть не должно. Такая проблема решается введением порога доли площади выдела, его необходимо преодолеть, тогда порода будет признана субдоминантной для данного выдела.
Заключение
Предложенный метод тематической обработки с использованием мультиспектральных изображений среднего и высокого пространственного разрешения совместно с данными текстурного анализа панхроматических изображений сверхвысокого разрешения позволяет - по сравнению с традиционной тематической обработкой с использованием только спектральных признаков - повысить точность сегментации лесных территорий на фоне различных природных и антропогенных объектов, а также классифицировать хвойные и лиственные леса. В целом улучшение точности тематической обработки, оцениваемой методом кросс-валидации, характерно для всех рассмотренных алгоритмов классификации. Наиболее значимое оно при использовании ансамблевого алгоритма классификации RF (случайный лес), наименьшее увеличение точности соответствует методу квадратичного дискриминантного анализа. Применение традиционного подхода к обработке мультиспектральных изображений WorldView 2 снижает точность классификации вследствие увеличения пространственных вариаций спектральных сигнатур между освещенными и затененными участками лесного полога. Однако при добавлении текстурной информации точность существенно повышается и превосходит соответствующие точности тематической обработки изображений среднего пространственного разрешения. Сравнение тематических карт доминантных и субдоминантных отделов лесной растительности, построенных по данным спектрально-текстурной обработки, с наземной информацией не позволяет выявить значимого эффекта как от привлечения текстурных признаков, так и от использования ансамблевого подхода: разброс по точности составляет не более 5 % около среднего значения 85 %. Анализ изображений выделов, дающих ошибку, показывает, что основной причиной несоответствий являются специфические погрешности данных наземной лесной таксации, и говорит о необходимости их актуализации.
СПИСОК ЛИТЕРАТУРЫ / REFERENCES
1. Воронина П.В., Мамаш Е.А. Классификация тематических задач мониторинга сельского хозяйства с использованием данных дистанционного зондирования MODIS // Вычислительные технологии. 2014. Т. 19, № 3. C. 76-102. Voronina P.V, Mamash E.A. Classification of Thematic Monitoring for Agriculture Problems Using Remote Sensing MODIS Data. Vychislitelnye tekhnologii [Computational Technologies], 2014, vol. 19, no. 3, pp. 76-102.
2. Колодникова Н.В. Обзор текстурных признаков для задач распознавания образов // Докл. ТУСУР. 2004. № 1(9). С. 113-124. Kolodnikova N.V. Overview of Textural
Features for Pattern Recognition Problems. Doklady TUSUR [Proceedings of TUSUR University], 2004, no. 1(9), pp. 113-124.
3. Потапов Л.А. Новые информационные технологии на основе вероятностных текстурных и фрактальных признаков в радиолокационном обнаружении малоконтрастных целей // Радиотехника и электроника. 2003. Т. 48, № 9. С. 1101-1119. Potapov A.A. New Information Technology in Radar Detection of Low-Contrast Targets Based on Probabilistic Texture and Fractal Features. Radiotehnika i elektronika [Journal of Communications Technology and Electronics], 2003, vol. 48, no. 9, pp. 1101-1119.
4. Фраленко В.П. Методы текстурного анализа изображений, обработка данных дистанционного зондирования Земли // Программные системы: теория и приложения. 2014. Т. 5, № 4(22). С. 19-39. Fralenko V.P. Methods of Image Texture Analysis, Earth Remote Sensing Data Processing. Programmnyye sistemy: teoriya i prilozheniya [Program Systems: Theory and Applications], 2014, no. 4(22), pp. 19-39.
5. Шабанов Н.В., Барталев С.А., Ерошенко Д.А., Плотников Ф.В. Развитие возможностей дистанционной оценки индекса листовой поверхности по данным MODIS // Современные проблемы дистанционного зондирования Земли из космоса. 2018. Т. 15, № 4. С. 166-178. Shabanov N.V., Bartalev S.A., Eroshenko F.V, Plotnikov D.E. Development of Capabilities for Remote Sensing Estimate of Leaf Area Index from MODIS Data. Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa [Current problems in remote sensing of the Earth from space], 2018, vol. 15, no. 4, pp. 166-178. DOI: https://doi.org/10.21046/2070-7401-2018-15-4-166-178
6. Beguet B., Guyon D., Boukir S., Chehata N. Automated Retrieval of Forest Structure Variables Based on Multi-Scale Texture Analysis of VHR Satellite Imagery. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, vol. 96, pp. 164-178. DOI: https://doi. org/10.1016/j.isprsjprs.2014.07.008
7. Breiman L. Random Forests. Machine Learning, 2001, vol. 45, iss. 1, pp. 5-32. DOI: https://doi.org/10.1023/A:1010933404324
8. Biynolfsson P., Nilsson D., Torheim T., Asklund T., Thellenberg Karlsson C., Tiygg J., Nyholm T., Garpebring A. Haralick Texture Features from Apparent Diffusion Coefficient (ADC) MRI Images Depend on Imaging and Pre-Processing Parameters. Scientific Reports, 2017, vol. 7, art. 4041. DOI: https://doi.org/10.1038/s41598-017-04151-4
9. Chen J.M., Cihlar J. Retrieving Leaf Area Index of Boreal Conifer Forests Using Landsat TM images. Remote Sensing of Environment, 1996, vol. 55, iss. 2, pp. 153-162. DOI: https://doi.org/10.1016/0034-4257(95)00195-6
10. Cover T.M., Hart P.E. Nearest Neighbor Pattern Classification. IEEE Transactions on Information Theory, 1967, vol. 13, no. 1, pp. 21-27. DOI: https://doi.org/10.1109/ TIT.1967.1053964
11. Dmitriev E.V, Sokolov A.A., Kozoderov W, Delbarre H., Melnik P.G., Donskoi S.A. Spectral Texture Classification of High-Resolution Satellite Images for the State Forest Inventory in Russia. Remote Sensing for Agriculture, Ecosystems, and Hydrology XXI, 2019, vol. 11149, art. 111491J. DOI: https://doi.org/10.1117/12.2532965
12. Eklundh L., Harrie L., Kuusk A. Investigating Relationships between Landsat ETM+ Sensor Data and Leaf Area Index in a Boreal Conifer Forest. Remote Sensing of Environment, 2001, vol. 78, iss. 3, pp. 239-251. DOI: https://doi.org/10.1016/S0034-4257(01)00222-X
13. Ferreira M.P., Wagner F.H., Aragâo L.E.O.C., Shimabukuro Y.E., de Souza Fil-ho C.R. Tree Species Classification in Tropical Forests Using Visible to Shortwave Infrared WorldView-3 Images and Texture Analysis. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, vol. 149, pp. 119-131. DOI: https://doi.org/10.1016/j.isprsjprs.2019.01.019
14. Franklin S.E. Remote Sensing for Sustainable Forest Management. Boca Raton, FL, CRC press, 2001, p. 424. DOI: https://doi.org/10.1201/9781420032857
15. Goldberger J., Hinton G.E., Roweis S., Salakhutdinov R.R. Neighbourhood Components Analysis. Advances in Neural Information Processing Systems 17: Proceedings of the 2004 Conference. Ed. by L. Saul, Y. Weiss, L. Bottou. Cambridge, MA, MIT Press, 2004, pp. 513-520.
16. Hadi H., Korhonen L., Hovi A., Ronnholm P., Rautiainen M. The Accuracy of Large-Area Forest Canopy Cover Estimation Using Landsat in Boreal Region. International Journal of Applied Earth Observation and Geoinformation, 2016, vol. 53, pp. 118-127. DOI: https://doi.org/10.1016/ijag.2016.08.009
17. Haralick R.M. Statistical and Structural Approaches to Texture. Proceedings of the IEEE, 1979, vol. 67, no. 5, pp. 786-804. DOI: https://doi.org/10.1109/PR0C.1979.11328
18. Haralick R.M., Shanmugam K., Dinstein I. Textural Features for Image Classification. IEEE Transactions on Systems, Man, and Cybernetics, 1973, vol. SMC-3, no. 6, pp. 610-621. DOI: https://doi.org/10.1109/TSMC.1973.4309314
19. Lechner A.M., Foody G.M., Boyd D.S. Applications in Remote Sensing to Forest Ecology and Management. One Earth, 2020, vol. 2, iss. 5, pp. 405-412. DOI: https://doi. org/10.1016/j.oneear.2020.05.001
20. Paletto A., Tosi V. Forest Canopy Cover and Canopy Closure: Comparison of Assessment Techniques. European Journal of Forest Research, 2009, vol. 128, iss. 3, pp. 265-272. DOI: https://doi.org/10.1007/s10342-009-0262-x
21. Sexton J.O., Song X.-P., Feng M., Noojipady P., Anand A., Huang C., Kim D.-H., Collins K.M., Channan S., DiMiceli C., Townshend J.R. Global, 30-m Resolution Continuous Fields of Tree Cover: Landsat-Based Rescaling of MODIS Vegetation Continuous Fields with Lidar-Based Estimates of Error. International Journal of Digital Earth, 2013, vol. 6, iss. 5, pp. 427-448. DOI: https://doi.org/10.1080/17538947.2013.786146
22. Srivastava S., Gupta M.R., Frigyik B.A. Bayesian Quadratic Discriminant Analysis. Journal of Machine Learning Research, 2007, vol. 8, pp. 1277-1305.
23. Tharwat A. Linear vs. Quadratic Discriminant Analysis Classifier: A Tutorial. International Journal of Applied Pattern Recognition, 2016, vol. 3, no. 2, pp. 145-180. DOI: https://doi.org/10.1504/IJAPR.2016.079050
вклад авторов: Все авторы в равной доле участвовали в написании статьи Authors' Contribution: All authors contributed equally to the writing of the article
