Описание и реконструкция строения порового пространства почвы с помощью корреляционных функций
Апробируется метод описания и реконструкции порового пространства почв с помощью корреляционных функций. Процесс реконструкции является лучшим способом проверки потенциального дескриптора почвенного строения. В качестве объектов исследования выбраны шлифы, отражающие восемь типов строения порового пространства зональных суглинистых почв и почвообразующих пород на территории Русской равнины и содержащих поры различной формы и ориентации. С помощью морфологического анализа и сравнения оригинальных изображений с их реконструкциями по корреляционным функциям, выполненных с помощью алгоритма “отжига”, показано, что на данном этапе метод реконструкции подходит для изометричных типов порового пространства с изрезанными, слабоизрезанными и округлыми порами. Двухточечные корреляционные функции, рассчитанные ортогональным методом, оказались различными для всех типов порового пространства и отразили пористость, удельную поверхность и корреляции на различных расстояниях. Полученные результаты позволяют сделать вывод, что описание порового пространства почв и грунтов с помощью корреляционных функций является перспективным. Показаны дальнейшие направления развития данного метода для описания строения порового пространства и определения физических свойств почв.
Введение
Изучение строения почвы, ее динамики и связи с другими различными факторами и параметрами является основной задачей почвоведения. Упаковка почвенных частиц и микроагрегатов, а также строение порового пространства определяют такие важные свойства почвы, как водо- и воздухопроницаемость, прочность, ОГХ, которые в свою очередь влияют на условия жизнедеятельности почвенной мезофауны и микроорганизмов, протекание различных химических процессов, перенос почвенных растворов, питательных веществ и частиц различного происхождения. С помощью оптической и электронной микроскопии, а в последнее время компьютерной микро-томографии можно получать двух- и трехмерные изображения структуры почвы и отдельных почвенных агрегатов с высоким разрешением (до нескольких микрометров). Для количественного описания структуры и порового пространства в шлифах и сколах почвы используются различные характеристики размеров, формы и ориентации плоских срезов пор и агрегатов.
Для количественной характеристики трехмерной структуры геофизической среды и ее связности можно применять теорию локальной пористости [12], распределение пор по размерам, извилистость и другие параметры [21]. Однако все эти полезные характеристики лишь косвенно описывают структуру пористой среды, универсальный дескриптор (под этим термином мы понимаем количественный показатель, отражающий все необходимые признаки) порового пространства и структуры почвы на данный момент не существует.
Чтобы описать структуру какого-либо объекта в идеале ее необходимо математически представить в виде некоторой функции или набора параметров такой функции. Помимо существования такой постановки задачи, необходимо решать и обратную — восстанавливать структуру по ее функции, что обычно называется реконструкцией [15, 28]. Наличие дескриптора подразумевает возможность восстановления структуры по его параметрам. Качество данной процедуры (сравнение с оригиналом) есть основной критерий отбора идеального дескриптора. Очевидно, что для математической реконструкции структуры почвы
в целом недостаточно стандартного набора почвенных характеристик.
В данной статье на основе численного моделирования порового пространства восьми образцов почвы (двухмерных шлифов) проведена проверка возможности описания и восстановления порового пространства и структуры почвы по двухточечным корреляционным функциям методом “отжига” (simulated annealing, SA), а также показаны сильные стороны и недостатки данного подхода, перспективы его использования в будущем.
Корреляционные функции и процедура реконструкции
Существует много различных корреляционных функций, применимых к описанию многофазных сред (кластерная, линейная, хорды и др.) [13, 28], однако если не оговорено иначе, под этим понятием будем подразумевать вероятность нахождения п точек в одной и той же фазе. Рассмотрим некоторую гипотетическую двухфазную среду в «-мерном Евклидовом пространстве что есть не что иное, как вероятность того, что случайно выбранная точка принадлежит фазе.
Корреляционная функция с п = 2 определяется как (рис. 1):
Рис. 1. Пример расчета двухточечных корреляционных функций в бинаризированной среде: если обе точки„ находятся в одной фазе (фазы помечены индексами b — черная, w — белая), то 52 = 1, если в разных, то S2 = 0.
Для двухфазной статистически однородной среды соотношение корреляционных функций обоих фаз можно записать в следующем виде.
Автокорреляционная функция записывается как (откуда видно, что данная корреляционная функция не различает фазы, то есть улучшение описания среды за счет вычисления независимых корреляционных функций для различных фаз невозможно): и является наиболее важной функцией для описания случайных сред [15]. Данное выражение имеет смысл вероятности нахождения точек.
Автоковариационная функция (нормированная через объемные доли каждой из фаз) может быть записана через автокорреляционную (уравнение (9)) функцию в виде.
Такие функции обладают некоторыми свойствами, наиболее важными из которых можно считать необходимые свойства реализуемости (то есть условия, при которых данная функция может быть реализована в виде двухфазной среды). Наиболее полный анализ таких свойств можно найти в работах. Перечислим здесь лишь некоторые основные свойства двухточечной корреляционной функции.
Остальные свойства не являются достаточно явными и в целом выполняются при восстановлении методом SA. Таким образом, двухточечная корреляционная функция является наиболее изученной, для нее существует множество аналитических решений для упрощенных сред (например, пересекающиеся и не пересекающиеся круги/сферы и пр. [27, 28]).
Считается установленным фактом, что любая структура и ее свойства могут быть полностью описаны бесконечным набором «-точечных корреляционных функций (уравнение (1)). Однако на практике вычисление корреляционных функций сп>2 сопряжено с множеством трудностей, что связано в основном с вычислительными ресурсами. Яонг и Торквато [29] приводят доводы в пользу того, что даже уже при п = 3 сложность расчетов не окупается точностью получаемой информацией, что приводило к тому, что на практике пользовались лишь функцией S2. В последнее время было показано, что для описания и восстановления гетерогенных сред недостаточно двухточечных корреляционных функций [17, 18], что связано с некоторым (в целом небольшим) количеством вырожденных состояний. Учитывая, что под желаемой реконструкцией мы будем понимать создание статистически равнозначной пористой среды (а не точной ее копии), вырожденные состояния, скорее всего, не являются проблемой двухточечных корреляционных функций. Качество реконструкции так же можно улучшать за счет использования дополнительных функций, перечисленных выше [13, 17].
Для начала рассмотрим случай восстановления по заданному набору корреляционных функций, где а — тип функции. Разницу между статистическим описанием для двух сред можно рассчитать в виде суммы квадратов разностей между значениями функций, где /„а и /„а — значения функций для двух сред (при реконструкции это восстановленная и контрольная среды, соответственно). Определенную таким образом “энергию” Е можно рассматривать как некое состояние системы, для минимизации которой отлично подходит метод “отжига”, описанный ниже.
В случае исследования материалов и геофизических пористых сред данные о системе представляют собой цифровое представление в виде пикселей/вокселей. Для реконструкции достаточно взять случайную среду и начать переставлять в ней пиксели, на каждом шаге проверяя, как меняется энергия согласно уравнению (13). Однако ввиду незнания структуры, в начале процедуры имеет смысл делать больше перестановок, даже если они не приводят к уменьшению Е. Для данной цели отлично подходит следующий алгоритм, когда вероятность того, что случайная перестановка принимается, записывается как: где АЕ = Enew — ЕоШ. Начальная температура (трактовка согласно распределению Больцмана в нижней части уравнения (14)) устанавливается так, чтобы вероятность принятия случайной перестановки при АЕ > 0 равнялась 0.5. Считается, что охлаждение системы обратно пропорциональное логарифму Т(к) ос l/ln(fc) приводит к минимуму энергии, однако это требует значительного времени, а потому на практике используется более быстрый план охлаждения, а именно
Т{к)/Т{0) = Хк, где X меньше, но близко к единице [29]. Так, например, для статистически гомогенной и изотропной среды для двухточечной корреляционной функции уравнение (13) примет следующий вид.
Результат конструкдии/реконструкции будет также зависеть от того, каким образом рассчитываются значения функций. Наиболее простым и не требующим серьезных затрат машинного времени является ортогональный алгоритм, предложенный в работе Яонга и Торквато [29]. В таком случае вычисление функций проводится только по ортогональным направлениям перемещением отрезка различной длины.
Еще одним важным моментом конструкции/реконструкции является остановка процесса, то есть условие, когда мы считаем, что мы достигли необходимой точности. Существует два основных подхода: 1) выбор порогового значения АЕ, при котором процесс считается завершенным, или 2) некоторое значительное количество не принятых итераций (например, 105) подряд. Первый метод кажется более предпочтительным, так как позволяет оценить количественно степень соответствия двух сред. Данный показатель энергии может рассматриваться как среднеквадратическая ошибка метода при реконструкции. Джиао с соавт. [16] исследовали пороговые значения разности энергии. Для начала авторы определили, что квант энергии при перестановке одного пикселя составляет: где N — линейный размер системы. Тогда для некоторого порогового значения энергии Eth количество “неправильно” расположенных пикселей составит.
И тогда отношение “неправильно” расположенных пикселей к общему количеству пикселей данной фазы составит.
Вместо порогового значения энергии Eth можно задать значение отношения пикселей у = у/А. Тогда пороговое значение энергии для остановки “отжига” при достижении желаемой точности реконструкции можно получить из уравнение (18) в виде.
Общий ход работы можно разделить на несколько основных этапов: 1) проводим стандартную обработку типового изображения порового пространства почвы, получаем его микроморфо-логическую картину (оригинал); 2) осуществляем расчет корреляционной функции для данного типа порового пространства; 3) случайно заполняем реконструируемое пространство пикселями, представляющими собой твердую фазу или поро-вое пространство; 4) случайно переставляем два разнофазных пикселя и, рассчитав корреляционную функцию для полученного изображения (реконструкция), сравниваем ее с функцией оригинала (уравнение (15)), принимаем либо отклоняем
перестановку (см. уравнение (14)); 5) повторяем пункт 4 до тех пор, пока не достигнем необходимой точности реконструкции; 6) проводим сравнение оригинального и реконструированного изображений на основе морфометрического анализа. Подчеркнем, что процедура реконструкции выполняется нами лишь для оценки качества статистического дескриптора строения порового пространства (в настоящей работе это двухточечная корреляционная функция); примером ее прикладного использования может являться восстановление трехмерной структуры почвы по двухмерным срезам. Для описания общей структуры почвы и оценки различных ее свойств достаточно произвести расчет корреляционных функций реального изображения (шлифа).
Объекты и методы
В качестве объектов исследования выбраны шлифы, отражающие восемь типов строения порового пространства (ПП) зональных суглинистых почв и почвообразующих пород на территории Русской равнины. Данные типы ПП специфичны для основных педогенных и литогенных агрегатных структур в шлифах: массивной (не разделенной на агрегаты в пределах шлифа), трещиновато-массивной, массивно-трещиноватой, комковатой, зернистой, ореховатой, плитчатой и пластинчатой (рис. 2) [7]. Таким образом, отобранные шлифы охватывают широкий диапазон конфигураций порового пространства и отлично подходят для проверки потенциального дескриптора.
Каждое исследованное изображение порового пространства было получено цифровым сканированием шлифа с разрешением 0.01 мм на пиксель. В дальнейшем -из общего изображения вырезали наиболее репрезентативный участок размером 2000 на 2000 пикселей. Для каждого шлифа производили разделение порового пространства на две фазы — поры и твердое вещество. Разделение производили подборкой одиночного порогового значения по гистограмме интенсивностей пикселей в виде среднего значения между двумя пиками, соответствующими этим фазам. Чтобы исключить влияние прозрачных минеральных зерен и сохранить для подсчета тонкие макропоры D = = 0.02—2.0 мм, отличающиеся максимальным разнообразием формы и ориентации, из полученных бинарных изображений удаляли все детали шириной менее 20 пикселей. Разрешение полученного в результате изображения уменьшали таким образом, чтобы его размер составлял 200 на 200 пикселей, и именно по нему рассчитывали корреляционные функции и производили процедуру реконструкции изображения. Уменьшение изображения связано с ресурсоемкостью вычислений, размер в 200 на 200 пикселей является неплохим компромиссом между качеством картинки и скоростью реализации вычислительного алгоритма.
Рис. 2. Типы порового пространства пород и почв в шлифах. Структура: I — массивная; II — трещиновато-массивная; III — массивно-трещиноватая; IV — комковатая; V — зернистая; VI — ореховатая; VII — плитчатая; VIII — пластинчатая. Поры белые.
Расчет корреляционных функций и алгоритм “отжига” реализовывался с помощью специально разработанного программного кода, однопроцессорные вычисления проводили на узлах высокопроизводительных компьютерных систем (суперкомпьютер Чебышев в НИ ВЦ МГУ). Время, необходимое для реконструкции одного изображения, в среднем составляло порядка нескольких суток. Для каждого типа ПП производили реконструкцию в двойной повторности (рис. 3). Пороговое значение энергии устанавливалось равным 10-7, что согласно уравнение (19) соответствует доле “неправильно” расставленных пикселей от 1 до 2% (в зависимости от массовой доли порового пространства для каждого типа ПП) от их общего количества.
Во всех оригинальных и восстановленных бинарных изображениях проводили стандартный морфометрический анализ. Для каждой детали изображения (поры) измеряли площадь (S), периметр (Р), продольный (L) и поперечный (D) габариты, угол отклонения длинной оси детали от вертикали (показатель ориентации). Фиксировали также количество измеренных деталей (общий объем выборки). На основании измерений для каждой детали рассчитывали показатель R = 4nS/P2, характеризующий отличие контура детали от окружности, и показатель /= D/L, характеризующий изометричность контура. Кроме того, для каждой детали рассчитывали обобщенный фактор формы F= (4kS/P2 + D/L)/2. Этот показатель имеет определенные преимущества при характеристике формы пор в шлифах и позволяет различать поры с различными срезами от трещиновидных до округлых [8]. В связи с особенностями программного комплекса для расчета параметров порового пространства полученные реконструкции и оригинальные изображения размером 200 на 200 увеличивали до размеров 360 на 360 пикселей.
Для комплексной характеристики строения порового пространства в оригинальных и восстановленных изображениях были определены следующие показатели:
1. Распределение пор по форме. Его описывали пятью частотами: содержанием трещиновидных (0 < F< 0.2), вытянутых изрезанных (0.2 < F< 0.4), изометричных изрезанных (0.4 < F< 0.6), изомет-ричных слабоизрезанных (0.6 < F < 0.8) и округлых (0.8 < F< 1.0) макропор (в процентах от общей численности измеренных пор).
2. Распределение пор по ориентации. Его описывали тремя частотами: содержанием пор с вертикальной, наклонной и горизонтальной ориентировкой (в процентах от общей численности измеренных пор). Первую группу составляли поры с показателями ориентации в интервалах 0—30° и 150—180°, вторую группу — с показателями ориентации в интервалах 30—60° и 120—150°, третью группу — с показателями ориентации в интервале 60-120°.
3. Средние арифметические значения и величины стандартного отклонения площади, периметра и габаритов пор.
Результаты и обсуждение
Визуальное сравнение оригинальных изображений и их реконструкций показало, что разработанный метод описания и восстановления (реконструкции) строения порового пространства почвы во многих случаях проявил себя достаточно успешно (рис. 3). К аналогичному выводу приводит сравнение формы и ориентации пор в изображениях. Так, из восьми исследованных типов ПП в двух типах (I и ГУ) оригинальные и реконструированные изображения имеют однотипные распределения пор по форме и ориентации, в трех типах (И, III и V) различия касаются либо формы, либо ориентации пор. В то же время для трех типов (VI, VII и VIII) различия между оригинальными изображениями и их реконструкциями оказались весьма существенными (табл. 1).
Для количественной оценки сходства и различия изображений в пределах каждого типа ПП был проведен кластерный анализ по трем объектам: оригиналу и двум реконструкциям. Анализ проводили по программам Рожкова [6]. Согласно полученным данным во всех типах ПП внутреннее сходство между реконструкциями равно или выше, чем между реконструкциями и оригиналом (табл. 2). При этом сходство между реконструкциями и оригиналами возрастает в следующем порядке: горизонтально вытянутые поры < поры, вытянутые в различных направлениях изометричные поры.
Помимо показателей формы и ориентации имеет смысл сопоставить прямые морфометрические параметры пор в разных изображениях (табл. 3). Измерения показали, что, как правило, в оригинале средние арифметические значения параметров меньше или равны реконструкциям. А варьирование морфометрических параметров зависит от формы пор. Если в изображении присутствуют трещиновидные или сильно вытянутые поры (типы III, VI, VII), то в оригинале стандартные отклонения превышают или равны реконструкциям. Если поры имеют изометричную форму (типы IV и V), то в оригинале стандартные отклонения меньше или равны реконструкциям. В целом по морфометрическим параметрам пор наиболее удачная реконструкция получена для массивной почвенной структуры с разрозненными порами округлой формы (тип I).
Рис. 3. Результаты реконструкции для восьми разных морфологических типов порового пространства почв:
А — оригинальное изображение пор в шлифе; Б и В — две повторности реализации моделирования корреляционными функциями. Поры черные.
Детальное изучение полученных реконструкций показывает возможности и ограничения предложенного метода расчета корреляционных функций. Так, для типов I и IV ПП реконструкции хорошо соответствуют оригиналам визуально и по набору поровых показателей. В то же время для типов V—VIII качество реконструкции существенно снижается благодаря большому количеству разно-ориентированных трещин на оригинальном изображении.
Таблица 1. Распределение пор по форме и ориентации в оригинальных изображениях и реконструкциях строения порового пространства в шлифах
Указанные недостатки реконструкций являются результатом осреднения корреляционных функций по пространству, что приводит к потере информации о форме и ориентации сильно вытянутых пор. Помимо осреднения функций, отсутствие трещино-видных пор на реконструкциях может быть вызвано иерархичностью порового пространства (то есть его статистической не стационарностью), однако данная тема требует более детальной проработки в будущем.
К числу недостатков метода можно также отнести занижение общего количества и завышение
Таблица 2. Уровни сходства между оригинальными изображениями и реконструкциями строения порового пространства в шлифах (%). Результаты кластерного анализа по восьми показателям (5 классов формы и 3 класса ориентации) площади пор на реконструкциях по сравнению с оригиналами.
Таблица 3. Морфометрические показатели пор в оригинальных изображениях и реконструкциях строения порового пространства в шлифах
Первое объясняется присутствием на реконструкциях 1—2% (от общей площади изображения) однопиксельных образований, которые были в дальнейшем удалены при морфометрическом анализе. Причиной второго недостатка может быть изменение габаритов вследствие округления формы исходно изрезанных и вытянутых изрезанных пор, а также изменение площади пор в результате принятой смены масштабирования.
В то же время на рис. 4 показано хорошее совпадение корреляционных функций оригинального и реконструированного изображений для всех типов ПП, что указывает на перспективные возможности алгоритма “отжига”. Все выявленные недостатки реконструированных изображений в первую очередь являются следствием изотропного расчета корреляционных функций (переход от уравнения (6) к уравнению (7)) и ортогональным методом дискретизации аргумента функции. Такой подход обладает некоторыми значительными недостатками [20], однако он позволяет успешно рассчитывать корреляционные функции для изометричных многофазных сред. Эффективного метода расчета с учетом анизотропии на данный момент не существует. Следует отметить, что хотя процесс реконструкции требует значительного машинного времени, расчет самих функций может быть с легкостью реализован на обычном персональном компьютере.
Сравнение корреляционных и автоковариационных функций (рис. 5) показало, что даже без учета анизотропии данные показатели значительно различались для всех типов порового пространства.
Рис. 4. Рассчитанные корреляционные функции для оригинальных (7) и реконструированных (реконструкции 1 и 2 — 2, 3) изображений для восьми типов ПП. Очевидно отличное совпадение значений во всех случаях.
На рисунке значение ^(О) соответствует пористости, наклон функции (производная в точке £2(0)) является показателем удельной поверхности [11], а флуктуации на всем протяжении функций отражают изменение корреляций в зависимости от расстояния. Наклон графиков функций S2 в точках г = 0 (рис. 4, 5) хорошо соотносится с визуальной оценкой поверхности раздела порового пространства и твердой фазы в
шлифах (рис. 2, 3). Использование полученных таким образом данных для определения различных свойств почв и грунтов — электрической проводимости, проницаемости, теплопроводности, модуля упругости, коэффициента диффузии, и т.п. с использованием так называемых точных пределов представляет несомненный интерес для физики почв.
С помощью процедур реконструкции можно решать еще одну важную задачу — восстановление трехмерной структуры по двухмерным данным [30].
Рис. 5. Сравнение корреляционных (А) и автоковариационных (Б) функций для восьми (1—8) типов порового пространства (рассчитаны по оригинальным шлифам).
Такая методика может использоваться в ситуациях, где разрешения томографов недостаточно, или размер образцов слишком велик и разнороден (при неоднородности разных масштабов можно сочетать данные томографии с 2D SEM изображениями). Не менее интересной и важной задачей является создание материалов и сред (например, идеальных для сельского хозяйства структур почвы в разных горизонтах), обладающих некоторыми необходимыми свойствами, при помощи аналитически выведенных функций [15, 28]. Уменьшение итоговой энергии ниже 10-7 было затруднительно ввиду не идеальности алгоритма случайной перестановки пикселей. В будущем возможно значительное улучшение точности реконструкции за счет эксплуатации модифицированного DPN алгоритма [26].
Таким образом, хотя на нынешнем этапе метод математической реконструкции порового пространства почвы далек от совершенства, его можно с успехом использовать для почв и пород с изометричными (в двумерных срезах) порами. Такие поры характерны для массивной структуры суглинистых почвообразующих пород с биогенными порами округлой формы и для комковатой структуры гумусовых горизонтов почвы с изометричными изрезанными порами упаковки комковатых агрегатов. Улучшение описания структуры и порового пространства почв и грунтов с помощью корреляционных функций можно развивать по трем направлениям: 1) совместное использование других типов функций (кластерная, линейная, хорды и др.), 2) разработка метода учета анизотропии и 3) использование многоточечных корреляционных функций (с п > 2). В будущем необходимо проведение реконструкции и анализа изображений без перемасштабирования на больших размерах расчетной сетки, в том числе в трех измерениях по данным рентгеновской микро-томографии.
Заключение
Описание порового пространства почв с помощью корреляционных функций является перспективным и позволяет реконструировать строение почвенной массы по данным значениям функции. Процесс реконструкции является лучшим способом проверки потенциального дескриптора почвенной структуры. Морфометрический анализ оригинальных (в шлифах) и реконструированных изображений показал, что на данном этапе разработки методов вычисления корреляционной функции наиболее полное описание возможно для изотропных почвенных структур с изометричными изрезанными, изометричными слабоизрезанными и округлыми порами. Двухточечные корреляционные функции, рассчитанные ортогональным методом для восьми изображений шлифов почвы, оказались различными для всех типов ПП и отразили пористость, удельную поверхность и корреляции на различных расстояниях. Стандартный морфометрический анализ позволяет проводить оценку качества проведенной реконструкции. Использованный подход к поиску дескриптора порового пространства и структуры почвы может найти практическое применение при моделировании процессов влагопереноса и газообмена, при мониторинге почвенного структурного состояния, для количественной диагностики и оценки степени деградации почвы.
Авторы сердечно благодарят к. т. н. А.Н. Ляхова, к. ф.-м. н. П.А. Брызгалова и проф. В.В. Воеводина за помощь в организации высокопроизводительных вычислений.