Рус Eng Cn Перевести страницу на:  
Please select your language to translate the article


You can just close the window to don't translate
Библиотека
ваш профиль

Вернуться к содержанию

Кибернетика и программирование
Правильная ссылка на статью:

Экспериментальное сравнение алгоритмов кластеризации в задаче группировки данных о грозовых разрядах

Беликова Марина Юрьевна

старший преподаватель, ФГБОУ ВО «Горно-Алтайский государственный университет»

649000, Россия, Республика Алтай, г. Горно-Алтайск, ул. Ленкина, 1

Belikova Marina Yur'evna

Senior Lecturer of the Department of Physics and Informatics, Gorno-Altai State University

649000, Russia, respublika Altai, g. Gorno-Altaisk, ul. Lenkina, 1

BelikovaMY@yandex.ru
Другие публикации этого автора
 

 
Каранина Светлана Юрьевна

кандидат физико-математических наук

доцент, ФГБОУ ВО «Горно-Алтайский государственный университет»

649000, Россия, Республика Алтай, г. Горно-Алтайск, ул. Ленкина, 1

Karanina Svetlana Yur'evna

PhD in Physics and Mathematics

Associate Professor, Department of Physics and Informatics, Gorno-Altai State University

649000, Russia, respublika Altai, g. Gorno-Altaisk, ul. Lenkina, 1

krechetovas@yandex.ru
Другие публикации этого автора
 

 
Глебова Алена Викторовна

старший преподаватель, ФГБОУ ВО «Горно-Алтайский государственный университет»

649000, Россия, Республика Алтай, г. Горно-Алтайск, ул. Ленкина, 1

Glebova Alena Viktorovna

Senior Lecturer, Department of Physics and Informatics, Gorno-Altai State University

649000, Russia, respublika Altai, g. Gorno-Altaisk, ul. Lenkina, 1

glebova-alena-1991@yandex.ru
Другие публикации этого автора
 

 

DOI:

10.25136/2644-5522.2018.1.25261

Дата направления статьи в редакцию:

19-01-2018


Дата публикации:

26-01-2018


Аннотация: В данной статье представлены результаты экспериментального сравнения кластерного анализа данных о грозовых разрядах алгоритмами k-means, dbscan и иерархическими агломеративными алгоритмами, где для вычисления межкластерного расстояния используются методы ближайшего соседа, полной и средней связи и метод Уорда. Произведена оценка влияния нормировочных параметров на количество кластеров, определяемых рассматриваемыми алгоритмами на тестовой выборке. В качестве тестовой выборки были использованы данные о времени регистрации и координаты грозовых разрядов, зафиксированные Всемирной сетью локализации молний (World Wide Lightning Location Network, WWLLN). Построение группировочных решений выбранными алгоритмами кластеризации производились с помощью па-кетов кластерного анализа Nbclust, dbscan и fpc, разработанных на языке R. Показано, что выбор значений нормировочных параметров имеет существенное влияние на количество выделяемых кластеров из рассматриваемой выборки при использовании иерархических алгоритмов кластеризации (особенно для метода ближайшего соседа). Выбор нормировочных параметров практически не влияет или оказывает незначительное влияние на результаты кластеризации грозовых разрядов с помощью алгоритмов k-means и dbscan. Наилучшее согласование с экспертной оценкой получено для алгоритма dbscan при нормировочных параметрах соответствующих линейным размерам грозовой конвективной ячейки 100 км и времени существования от 30 минут до часа.


Ключевые слова:

интеллектуальный анализ данных, алгоритмы кластеризации, k-means, dbscan, иерархические алгоритмы, качество кластеризации, индекс силуэта, нормализация данных, грозовой разряд, WWLLN

Исследование выполнено при финансовой поддержке фонда РФФИ в рамках научного проекта № 16-47-040081 р_а.

Abstract: The authors present the results of an experimental comparison of the cluster analysis of thunderstorm data using the algorithms of k-means, dbscan and hierarchical agglomerative algorithms, where closest neighbor, full and medium coupling methods and the Ward method are used to calculate the intercluster distance. The influence of the normalization parameters on the number of clusters determined by the algorithms under consideration on the test sample is estimated. Data on the time of registration and the coordinates of lightning discharges recorded by the World Wide Lightning Location Network (WWLLN) were used for test purposes. The construction of grouping solutions by the chosen clustering algorithms was carried out with the help of the Nbclust, dbscan, and fpc cluster analysis packages developed in the R language. The article showns that the choice of the values of the normalization parameters has a significant effect on the number of clusters allocated from the sample under consideration using hierarchical clustering algorithms (especially for method of the nearest neighbor). The choice of the normalizing parameters has practically no effect or has a negligible effect on the results of lightning cluster clustering using the k-means and dbscan algorithms. The best agreement with expert judgment was obtained for the dbscan algorithm with normalizing parameters corresponding to linear dimensions of a thunderstorm convective cell of 100 km and a period of time of 30 minutes to an hour.


Keywords:

data mining, clustering algorithms, k-means, dbscan, hierarchical algorithms, cluster validity, average silhouette width (asw), data normalization, lightning, WWLLN

Введение

В настоящее время отмечается значительный интерес к применению методов кластерного анализа для обработки данных грозопеленгационных систем (ГПС) [1-6]. ГПС регистрируют электромагнитный сигнал отдельных грозовых разрядов и определяют их координаты, время и, возможно, некоторые характеристики, например, полярность разряда, амплитуда тока, количество компонент молнии. Для решения задач в области атмосферного электричества необходима информация не только о разрядах молний и их плотности, но и информация о грозах, которые рассматриваются как кучево-дождевые облака, продуцирующие молнии. В частности, к таким задачам относится задача оценки вклада грозовых генераторов в глобальную электрическую цепь [5, 7]. Особенностью гроз является наличие одного или нескольких разрядов за время их развития и существования. Для получения информации о грозах необходима группировка данных о грозовых разрядах, регистрируемых ГПС. Методом решения задачи группировки данных о молниевых разрядах является кластерный анализ.

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

Результаты кластеризации и, следовательно, выводы, основанные на этих результатах, могут существенным образом зависеть от заданных значений параметров нормировки (в том случае если они используются), выбранного алгоритма и набора его входных параметров. Для подбора оптимального параметра (параметров) алгоритма могут быть использованы так называемые внутренние критерии качества группировки, значения которых зависят от разброса объектов внутри полученных кластеров и расстояний между кластерами. Наилучшему значению критерия качества соответствует оптимальный параметр (параметры) рассматриваемого алгоритма кластеризации и группировка объектов, полученная при этом, признается наилучшей [8]. Одним из наиболее часто используемых внутренних критериев качества для алгоритмов кластеризации является индекс «силуэта» (average silhoutte index, asw) [9]. С помощью индекса asw можно сравнивать и результаты группировки объектов, полученные разными алгоритмами.

Задача кластерного анализа грозовых разрядов изначально не выглядит сложной и выбор алгоритма кластеризации для ее решения зависит в основном от предпочтений исследователя. Для кластеризации данных о грозовых разрядах использовались следующие алгоритмы: алгоритм на основе метода разбиений k-means [1, 2], алгоритм на основе сеточного метода [3], алгоритм на основе графового метода [4], алгоритм на основе оценки плотности dbscan [5] и комбинированный алгоритм на основе метода ближайшего соседа и модального анализа [6]. Отметим, что кластеризация данных о грозовых разрядах с разными входными параметрами проводилась только для алгоритма dbscan [5].

В настоящей работе рассматриваются результаты кластеризации молниевых разрядов с помощью алгоритмов k-means [10], dbscan [11] и иерархических агломеративных алгоритмов (single - метод ближайшего соседа, complete – метод полной связи, average – метод средней связи, ward.D2 – метод Уорда), которые различаются методами вычисления межкластерного расстояния [12]. Готовые программные реализации этих алгоритмов являются одними из самых доступных. Например, в широко применяемой в настоящее время для анализа данных среде R перечисленные алгоритмы кластеризации, а также индексы качества кластеризации, в том числе и индекс asw, реализованы в виде функций, включенных в пакеты кластерного анализа.

Предыдущие исследования показали необходимость дополнительных вычислительных экспериментов для определения оптимальных входных параметров алгоритма dbscan, с целью повышения согласованности результатов кластеризации молниевых разрядов с региональными особенностями грозовой актиминости (среднегодовое количество и среднегодовая продолжительность гроз) [13]. Использование иерархических алгоритмов для кластеризации грозовых разрядов интересно с позиции определения иерархической структуры грозы (многоячейковые грозы). Ранее эти алгоритмы для группировки грозовых разрядов не применялись. Алгоритм k-means является наиболее часто используемым для решения различных прикладных задач кластеризации, несмотря на известные недостатки этого алгоритма: необходимость задавать количество кластеров, чувствительность к «шуму» и выбросам, неспособность выделять кластеры разной структуры. Для данных о грозовых разрядах есть большая вероятность появления одиночных разрядов на значительном расстоянии от некоторой группы разрядов. Так как алгоритм k-means выделяет кластеры центра, являющегося средним геометрическим характеристик объектов, то может быть получена не лучшая в смысле интерпретируемости результатов группировка грозовых разрядов. Поэтому вопрос о применимости алгоритма k-means для группировки данных о грозовых разрядах требует дополнительного уточнения.

Проведение кластерного анализа указанными алгоритмами в среде R предполагает нормировку исходных данных. В качестве нормировочных параметров при кластеризации грозовых разрядов, могут быть использованы линейные размеры (от 25 до 100 км) и время существования (от 20 до 60 минут) грозовой конвективной ячейки. В работе [6] предложено в качестве нормировочных параметров использовать линейные размеры равные 50 км и время существования 30 минут. Сравнение результатов кластеризации грозовых разрядов с другими значениями нормировочных параметров не производилось.

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

Исходные данные и методы

Для вычислительного эксперимента были использованы данные о грозовых разрядах, зарегистрированные Всемирной сетью локализации молний (World Wide Lightning Location Network, WWLLN) [5, 14] на территории Республики Алтай 7 августа 2017 года.

Сетью WWLLN на территории республики 7 августа 2017 года было зарегистрировано 2030 грозовых разрядов. Для упрощения визуального анализа результатов кластеризации на предмет их соответствия грозовым облакам был выбран часовой промежуток времени, соответствующий спутниковым данным по облачности грануле продукта MYD06 (MODIS/Aqua) [15] в 7:45 по всемирному скоординированному времени (ВСВ). Тематические продукты MODIS/Aqua являются результатом сканирования подспутниковой области в течение пяти минут. Однако, для тестирования результатов кластерного анализа молниевых разрядов сложно выбрать область так, чтобы было достаточное количество разрядов именно в пятиминутном интервале. Выбор временного промежутка в один час также обусловлен тем, что характер облачного покрова за указанный временной интервал значительно не изменился (по данным изображений облачности геостационарного спутника Meteosat [16]).

Для промежутка времени с 07:15 по 08:15 ВСВ (14:15-15:15 по местному времени) в выбранном регионе сетью WWLLN было зарегистрировано 34 молниевых разряда. Все молниевые разряды соответствуют облакам с низкой температурой верхней границы облаков и с содержанием ледяных частиц (рис. 1). Это характерные параметры грозовых облаков.

Визуальный анализ расположения грозовых разрядов позволяет выделить 5-7 групп кластеров. В северной части республики (участок между 52° и 53° с.ш. и 87° и 88° в.д.) наблюдается «сложный» кластер («северный» кластер), который состоит из трех групп разрядов, скорее всего, соответствующих многоячейковому грозовому облаку. В западной (участок между 50° и 52° с.ш. и 84° и 86° в.д.) и восточной (участок между 49° и 51° с.ш. и 89° и 90° в.д.) частях рассматриваемой территории, имеется по два кластера молниевых разрядов («западные» и «восточные» кластеры, соответственно). Эти кластеры, скорее всего, принадлежат разным конвективным ячейкам, которые предположительно определяются по параметрам облачности (рис.1). Вопрос объединения каждой пары «западных» и «восточных» кластеров в один «сложный» кластер является спорным, так как расстояние между парами групп разрядов достигает 100 км.

Рисунок 1. Параметры облачности по данным МYD06 (07:45) и молниевые разряды за временной интервал с 07:15 до 08:15 ВСВ

Сетью WWLLN регистрируются следующие показатели грозовых разрядов: дата, время, широта, долгота, погрешность и количество станций, в которых был зарегистрирован электромагнитный импульс. Для кластерного анализа были использованы координаты долготы, широты и времени.

В среде R описания выбранных 34 грозовых разрядов представлены в виде матрицы размерностью 34×3. Первые два столбца матрицы содержат вещественные числа, соответствующие координатам разрядов в десятичных градусах. Третий столбец матрицы содержит дату-время регистрации молниевых разрядов формате POSIXct. Формат POSIXct определяет дату-время в виде секунд, истекших с 1 января 1970 г. Данные о грозовых разрядах были преобразованы в секунды с 1 января 2017 г.

Построение группировочных решений рассматриваемыми алгоритмами кластеризации с оценкой их качества производилось с помощью пакетов кластерного анализа Nbclust [17], dbscan [18] и fpc [19], разработанных на языке R. Указанные пакеты содержат «комплексные» функции, реализующие алгоритмы кластеризации и позволяющие определить лучшие входные параметры алгоритмов с помощью внутренних индексов качества. Были использованы следующие функции (функция/пакет), возвращающие оптимальное число кластеров в соответствии с выбранным критерием качества asw: Nbclust/Nbclust для иерархических алгоритмов и алгоритма k-means, dbscan/dbscan для алгоритма dbscan и cluster.stats/fpc для вычисления индекса asw.

Вычислительный эксперимент и результаты

Целью вычислительного эксперимента является определение степени влияния на результат кластеризации априорных нормировочных значений и способа нормировки выборки молниевых разрядов, соответствующих временному интервалу один час (c 07:15 по 08:15 UTC).

При проведении кластерного анализа объектов исследования необходимо задание априорных нормировочных параметров dA - линейные размеры грозовой конвективной ячейки и tA - среднее время существования грозовой конвективной ячейки, которые показывают, что грозовые разряды расстояние между которыми меньше dA и разница во времени регистрации меньше tA должны образовывать одну группу. Выше было указано, что линейные размеры грозовой конвективной ячейки могут быть от 25 до 100 км и время существования - от 20 до 60 минут. С учетом точности определения места возникновения молниевого разряда сетью WWLLN, которая в среднем составляет 5.4 км, при кластеризации грозовых разрядов целесообразно использовать следующие значения априорных нормировочных параметров dA= 50 км, 75 км, 100 км и tA= 20 мин, 30 мин, 40 мин, 50 мин, 60 мин.

Нормирование априорными параметрами производится двумя способами. Первый способ подразумевает, что значения нормировочных параметров (dA,dA,tA) задаются в явном виде [6]. Во втором способе значения нормировочных параметров задаются как (1,1,C), где C определяется по формуле C = dA / tA и является коэффициентом связывающим расстояние в пространстве между событиями (грозовыми разрядами) и временем их регистрации [20]. Отметим, что вычисление расстояния между нормированными характеристиками объектов производится с помощью евклидовой метрики.

Результаты группировки молниевых разрядов для первого и второго способа задания априорных нормировочных параметров совпали. Таким образом, способ нормировки ([6] или [20]) не влияет на результаты кластеризации.

В таблице 1 для каждого набора априорных нормировочных параметров (dA,dA,tA), (1,1,C) и для каждого алгоритма приведены полученные в результате кластеризации количество кластеров (k) и величина индекса качества asw. Для алгоритма dbscan дополнительно указано количество разрядов, определенных как «шум» (noise).

Значения индекса asw для всех нормировочных параметров и для всех алгоритмов попадают в интервал от 0.71 до 1, что свидетельствует о наличии сильной (strong) структуры кластеров [10].

Таблица 1 – Результаты кластеризации данных о грозовых разрядах, зарегистрированных WWLLN 7 августа 2017 года на территории Республики Алтай

Иерархический алгоритм single

tA = 20 min

tA = 30 min

tA = 40 min

tA = 50 min

tA = 60 min

k

asw

k

asw

k

asw

k

asw

k

asw

50 км

15

0.82

3

0.85

4

0.87

4

0.88

4

0.89

75 км

15

0.84

15

0.82

3

0.83

3

0.85

4

0.87

100 км

15

0.86

15

0.84

15

0.82

3

0.83

3

0.85

Иерархический алгоритм complete

50 км

3

0.80

3

0.85

4

0.87

4

0.88

4

0.89

75 км

19

0.82

3

0.81

3

0.83

3

0.85

4

0.87

100 км

19

0.83

19

0.82

3

0.80

3

0.83

3

0.85

Иерархический алгоритм average

50 км

3

0.80

3

0.85

4

0.87

4

0.88

4

0.89

75 км

15

0.84

3

0.81

3

0.83

3

0.85

4

0.87

100 км

19

0.83

19

0.82

3

0.80

3

0.83

3

0.85

Иерархический алгоритм ward.D2

50 км

3

0.80

3

0.85

4

0.87

4

0.88

4

0.89

75 км

19

0.83

3

0.81

3

0.83

3

0.85

4

0.87

100 км

20

0.83

19

0.82

3

0.80

3

0.83

3

0.85

kmeans

50 км

3

0.80

3

0.85

3

0.86

3

0.87

3

0.88

75 км

3

0.75

3

0.81

3

0.83

3

0.85

3

0.86

100 км

4

0.72

3

0.78

3

0.80

3

0.83

3

0.85

dbscan

50 км

7

noise=5

0.78

6

noise=2

0.76

6

noise=2

0.78

6

noise=2

0.78

6

noise=2

0.78

75 км

7

noise=5

0.78

6

noise=2

0.72

6

noise=2

0.74

6

noise=2

0.76

6

noise=2

0.78

100 км

8

noise=3

0.76

7

noise=0

0.66

7

noise=0

0.70

6

noise=0

0.69

6

noise=0

0.70

Для используемых в работе [6] нормировочных параметров dA = 50 км и tA = 30 мин все рассматриваемые алгоритмы, кроме dbscan, разбивают выборку на три «крупных» кластера, разнесенных в пространстве друг от друга на большое расстояние. При этом пространственно полученные кластеры, объединили в себе все разряды в северной, западной и восточной частях рассматриваемой территории (например, рис.2). Таким образом, в один кластер были отнесены достаточно удаленные друг от друга разряды.

Алгоритм k-means только для параметров dA = 100 км и tA = 20 мин разбивает выборку на четыре кластера, при всех других параметрах нормировки – на три кластера. Это свидетельствует о малой зависимости результатов кластеризации от нормировочных параметров. Пространственно три выделенных кластера соответствуют группам разрядов в северной, западной и восточной частях территории республики (рис. 2).

В целом для всех рассматриваемых иерархических алгоритмов выявлена зависимость результатов кластеризации от нормировочных параметров, так как для разного набора нормировочных параметров получено различное количество кластеров. При этом набольшую зависимость от нормировочных параметров показал метод single (метод ближайшего соседа). Иерархические алгоритмы при tA = 20 минут и dA =75 км и 100 км выделяют более 15 кластеров, что не является приемлемым результатом для рассматриваемой выборки. При tA = 30 минут и dA = 75 км и 100 км результаты группировки аналогичны результатам алгоритма k-means (3 кластера, аналогичных представленным на рисунке 2). Близким к результатам экспертной оценки группировки грозовых разрядов можно считать выделение 4 кластеров при нормировочных параметрах dA = 50 км и tA = 40, 50 и 60 минут. При этом в один кластер объединяются два кластера в западной части и три кластера в северной части территории республики (эти кластеры можно отнести к одному многоячейковому грозовому облаку), восточные кластеры полностью согласуются с экспертной оценкой (рис.3).

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

Иерархические агломеративные алгоритмы рассматривались в работе на предмет их использования для группировки грозовых разрядов в виде иерархической структуры, с целью определения одноячейковых и многоячейковых грозовых облаков. Экспериментально выявлено, что иерархические агломеративные алгоритмы не подходят для таких целей. Кроме этого необходимо отметить и основные недостатки этих алгоритмов: высокая вычислительная сложность (порядка O(N3)) и, следовательно, их неприменимость для обработки больших массивов данных, а также неспособность выделять кластеры разной структуры.

Рисунок 3. Пример разбиения на кластеры молниевых разрядов с помощью алгоритма average (метод средней связи) и их соответствие грозовым облакам

Алгоритм dbscan выделяет разное количество кластеров (от 6 до 8 кластеров) при tA = 20 минут для всех рассматриваемых значений параметра dA и при dA = 100 км для всех параметров tA. Количество выделяемых алгоритмом dbscan кластеров близко к экспертной оценке в 5-7 кластеров. Следовательно, наблюдается слабая зависимость результатов кластеризации грозовых разрядов от параметров нормировки.

Алгоритм dbscan выделил соответственно экспертной оценке по два отдельных кластера в восточной части территории республики Алтай (рис. 4 а,б). Северный кластер разбивается на 3 кластера, что, в принципе соответствует возможности определения многоячейковой грозовой области. В западной части разряды одного из кластеров определены как «шум» при dA = 50 и 75 км и tA не менее 30 минут (рис. 4а). «Шумовые» разряды не выделяются при нормировочных параметрах dA = 100 км и tA не менее 30 минут (рис. 4 б), что наилучшим образом согласуется с экспертной оценкой.

А) при = 50 км и = 30 минут

Б) при = 100 км и = 30

Рисунок 4. Пример разбиения на кластеры молниевых разрядов с помощью алгоритма dbscan и их соответствие грозовым облакам

В ходе проведения кластерного анализа грозовых разрядов с помощью представленных в работе алгоритмов, авторы пришли к выводу о том, что процесс группировки грозовых разрядов должен предполагать два этапа: на первом этапе выделяются небольшие достаточно плотные кластеры и на втором этапе близкие кластеры объединяются в сложный кластер. При этом результаты вычислительного эксперименты показали, что для вычисления межкластерного расстояния лучше использовать метод средней связи (average method) или метод Уорда (Ward method).

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

Заключение

Исследовано влияние нормировочных параметров на результаты кластеризации данных о грозовых разрядах иерархическими агломеративными алгоритмами и алгоритмами k-means, dbscan.

На примере рассматриваемой выборки грозовых разрядов и алгоритмов кластеризации k-means, dbscan и иерархических алгоритмов продемонстрировано, что, несмотря на положительный опыт использования нормировочных параметров, соответствующих линейным размерам грозовой конвективной ячейки 50 км и времени ее существования 30 минут [6], эти параметры не позволили получить количество кластеров, согласующееся с экспертной оценкой.

Экспериментально показано, что выбор значений нормировочных параметров имеет существенное влияние на количество выделяемых кластеров с помощью иерархических алгоритмов кластеризации (особенно для метода ближайшего соседа). При этом иерархические алгоритмы имеют ряд недостатков и не могут быть использованы дл выделения иерархической структуры молниевых разрядов (одноячейковые или многоячейковые грозы). Выбор нормировочных параметров практически не влияет на результаты кластеризации грозовых разрядов с помощью алгоритмов k-means и dbscan. При этом количество кластеров молниевых разрядов, полученных в результате кластеризации алгоритмом k-means, не согласуется с их экспертной оценкой. Наилучшее согласование результатов кластеризации с экспертной оценкой получено алгоритмом dbscan при нормировочных параметрах соответствующих линейным размерам грозовой конвективной ячейки 100 км и времени существования от 30 минут до часа.

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

Библиография
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.