Метод расчета профилей распределения пористости и объемных долей материалов в пористой среде с помощью анализа данных рентгеновской микротомографии

УДК 532.546

Н. И. Рыжиков1'2, Д. Н. Михайлов2, В. В. Шако2

1 Московский физико-технический институт (государственный университет) 2Московский научно-исследовательский центр «Шлюмберже»

Метод расчета профилей распределения пористости и объемных долей материалов в пористой среде с помощью анализа данных рентгеновской микротомографии

В статье предложен новый метод расчета профилей распределения пористости и объемных долей материалов в пористой среде, основанный на анализе формы гистограмм отдельных сечений трехмерного изображения рентгеновской компьютерной микротомографии. В отличие от широко используемого метода установления порога, доля различных веществ в образце вычисляется без задания какого-либо априорного отсечения. Представлены расчеты профилей пористости и концентрации слагающих минералов для двух различных искусственных пористых образцов. Полученные результаты хорошо согласуются с макроскопическими параметрами исследуемых образцов и расчетами, выполненными с помощью коммерческого программного пакета СТАп. На примере исследования образца после закачки слабоконцентрированного глинистого раствора, показано, что разработанный метод применим для анализа распределения веществ малой концентрации.

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

1. Введение

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

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

В данной работе предложен новый метод расчета профилей распределения пористости и объемных долей материалов, в том числе привнесенных, в пористой среде с помощью анализа данных рентгеновской компьютерной микротомографии (X-rav micro Computed Tomography — XmCT).

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

При обработке микрокомпьютерной томографии чаще всего используется пороговая сегментация [6, 7], при которой, однако, выбор порога является физически необоснованным.

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

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

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

Расчет концентрации привнесенной в образец примеси являлся частью экспериментальной работы но исследованию динамики проникновения глинистого раствора и суспензий твердых частиц в пористую среду. Рентгеновская компьютерная микротомография образцов выполнялась с помощью прибора БкуБсап 1172 (производство БкуБсап).

2. Методика анализа данных

Трехмерное изображение образца, полученное с помощью рентгеновской компьютерной микротомографии, состоит из последовательных двумерных изображений поперек оси образца сечений. Расстояние между сечениями равно размеру пикселя. В данной работе размер пикселя составлял 2.5 5 мкм. Каждое сечение характеризуется 16-битными градациями серого цвета, которые соответствуют коэффициентам затухания рентгеновского луча в области расположения вокселя.

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

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

Распределение (или гистограмма) градаций серого на выделенном сечении представляется в виде зависимости числа пикселей от яркости серого цвета, рис. 2.

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

етвующих в данном сечении (1).

р(г) = ^ WiPi{z),

г=1

здесь г — значение яркости серого цвета, р(х) — функция плотности вероятности от всего сечения изображения, Рг(г) — функция плотности вероятности для г-й подобласти, п — общее число подобластей, на которое разбивается изображение, — весовой коэффициент для г-й подобласти.

300000 2Б0000 200000 1-50000 100000 50000 О

Песок /

/ ная гяия ¿3

Л / /

/

/

7000 6000 11000 13000 15000

Яркость серого

17000

10000

Рис. 2. Пример гистограммы сечения рентгеновской компьютерной микротомографии

Рис. 3. Аппроксимация нормализованной гистограммы отдельного материала с помощью нормального распределения

При этом весовые коэффициенты для функции плотности вероятности от определенной подобласти соответствуют относительной площади, занимаемой данной подобластью.

Например, на рис. 2 выделяются три подобласти, соответствующие участкам сечения занимаемым водой, песком и модифицированной глиной соответственно.

Суперпозиция (1), при известных функциях плотности вероятности для каждой подобласти, дает возможность рассчитать относительные площади каждой подобласти на сечении, решая обратную задачу (2):

р

г=1

0.

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

Форма гистограммы отдельного материала на сечении рентгеновской компьютерной микротомографии часто [9] имеет форму нормального распределения (гауесиана):

I (г) = А ехр

№

В нашем случае адекватность аппроксимации гистограммы гауееианом была проверена экспериментально на различных веществах и для различных образцов (пример представлен на рис. 3).

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

3. Модель переходной зоны

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

Рис. 4. Переходная зона

Рассмотрим поведение границы между материалами на сечении рентегновской микротомографии (рис. 4). Можно отметить, что профиль значений серого цвета в переходной зоне имеет линейный или близко к линейному вид. Тогда, зная значения серого цвета для пикселей граничащих материалов 22), значение серого цвета г^ ъ г-м пикселе переходной зоны может быть выражено как (4):

(т + 1 - г)г\ + гг2 , ,

гг = -Л-' №

т + 1

где т — число пикселей в переходной зоне.

Значение серого цвета для пикселей граничащих материалов — случайные величины, а следовательно, и их сумма тоже случайная величина. Для выделенного пикселя из переходной зоны функция плотности вероятности выражается через параметры функции гауссианов граничащих материалов:

р(*г) =-ехр ( -(г М

о%\/ 2тт

= (т + 1 - г)м1 + гм2 т + 1

2 (т + 1 - г)а2 + га2

я? =-н-,

г т + 1

где а1, Мь а2, М2 и а^ м — средние значения и дисперсии гауссианов, соответствующих первому граничащему материалу, второму и переходной зоне.

Функция плотности вероятности для всей переходной зоны будет представлять собой сумму функций плотности вероятности каждого пикселя по ширине переходной зоны:

т

Р^апзи

г

Зная формы функций плотности вероятности всех частей изображения — отдельных материалов и переходных зон, — можно найти их параметры, минимизируя следующую невязку:

к

£

з=1

п

П 2

- ^^ (У, А В, П) -

Н

г=1

^ 0,

С(х, А, В, П) = А ехр ^- (^ ,

где Н(г) — гистограмма анализируемого сечения, п — количество материалов (подобла-

на гистограмме. Суммарное число параметров минимизации (7) составляет п + 0.5п(п - 1).

Как отмечалось ранее, зная параметры каждой функции плотности вероятности, можно рассчитать соответствующие им площади.

Что касается переходных зон, то в литературе обсуждаются различные способы их интерпретации. Например, авторы [10] предложили распределять переходную зону между материалами пропорционально площадям, занятым «чистыми» материалами. Мы предполагаем, что поскольку количество пикселей по обе стороны от границы раздела материалов в среднем одинаково, то и площадь переходной зоны должна делиться пополам между граничащими материалами. Таким образом, площадь, соответствующая отдельным материалам на сечении, равна

1 п

где Я — площадь, занимаемая г-м материалом па сечении, Ба — площадь, соответствующая функции плотности вероятности г-го материала, втг]1 — площадь переходной зоны

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

чальное приближение при минимизации для ( I + 1)-го сечения. Кроме того, па результат

минимизации для (I + 1)-го сечения накладывается ограничение в виде 10% отклонения от результатов 1-го, т.е. для любого параметра минимизации А должно выполняться условие |1 — А1+1/А11 ^ 0.1, где I — порядковый помер сечения, на котором А\ является решением обратной задачи (7).

Такое ограничение физически оправдано, поскольку распределение материалов на двух соседних сечениях рентгеновской микротомографии должно различаться незначительно: расстояние между соседними сечениями (1-5 ц,т) значительно меньше характерных размеров неоднородностей пористой среды — пор и зерен скелета (20-100 ^т).

4. Результаты расчетов

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

Были использованы образцы цилиндрической формы длиной 3-5 см и диаметром 8 мм. Такой размер был специально подобран для достижения необходимого разрешения рентгеновской микротомографии.

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

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

Достоверность рассчитанных профилей концентрации проверялась путем сравнения (рис. 5) с данными, полученными обычным методом установки порога с помощью коммерческого программного пакета CTAn — Computed Tomography Analyser (фирма SkvScan, skvscan.be).

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

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

Образец № 2: искусственный образец с пористостью 33.2% и проницаемостью 8 Дарси, в который в ходе фильтрационных экспериментов был закачан однопроцентный раствор бентонитовой глины. Затем с противоположного конца осуществлена промывка образца водой, в результате чего его проницаемость частично восстановилась до 4 Дарси.

Не полное восстановление проницаемости свидетельствует о наличии в образце остаточной привнесенной глины.

Рентгеновская компьютерная микротомография образца показала, что количество глины, удерживаемой в поровом пространстве, небольшое, и ее гистограмма перекрывается гистограммами переходных зон. В этом случае для получения концентрации методом анализа гистограмм влияние переходной зоны предварительно уменьшалось путем использования специального фильтра Кувахара (Kuwahara, fill). Рассчитанный профиль концентрации глины представлен на рис. 6.

Отметим, что левый торец образца, через который осуществлялась закачка, был частично механически поврежден в ходе эксперимента. Поврежденный участок образца расположен от 0 до 1 мм (рис. 6). Далее наблюдается плавное снижение концентрации внесенной глины по мере удаления от точки ее закачки.

Полученные аналогичным способом профили концентрации примесных агентов (глина, твердые частицы и т.д.), закачанных в образец пористой среды, в дальнейшем могут быть

1

0,9 0.8 ¡2 0 7

га о.е

go5

<и

g04 О

0.3 0.2 0 1 о

1

0,9 0.8 S 0.7

« о.е

Ё 0,5

<U

§04 О йй 0,3

0.2

0.1

о

. «

tj в/ it \ 1 Vv

/ if /1 у

// -СТАп целестит

...............А— Л :тит

\ > ................iV у4

Та j % 1

т V*7 V -—»

4 1 4 2.4 3.4 4 Длина 4 5.4 6 образца, мм 4 7.4 8 4

A ft» ¡Г^ у |гуJ

варц РЧ

1 ЧХ л

и и

V Ч____ -

34 4.4 5.4

Длина образца, мм

0.9

0.8

S 0.7

0.6

£ 0.5

1>

0.4

о

и 0.3

0.2

0

..................j i * 1/

i i

т /У

ЭЛИТ ИТ

АГгал

3.4 4.4 5.4 64

Длина образца, мм

Рис. 5. Профили концентрации для трех минералов в модельном составном образце, рассчитанные методом установки порога

0.045

0.04

0.035

К э 0.03

р. 0.025

к И 0.02

0.015

0.01

0.005

0

f

\ t\ %

г ш W - ,

/ / ** v

/ 1 t v i ^ л 1 л ft*

'J

Г \ 1 4 У

- - - -

2

Длина,

Рис. 6. Профиль концентрации глины в пористом образце, рассчитанный методом анализа гистограмм. Левый торец образца (до 1 мм) частично поврежден в ходе эксперимента

использованы для проверки различных моделей кольматации (например, [12 15|).

На рис. 7 представлено сравнение профиля пористости, рассчитанного методом установки порога (коммерческий программный пакет СТАп), с профилем пористости, рассчитанным методом анализа гистограмм. Ключевым фактом в данном случае является не просто соответствие уровней пористости, которое легко может быть обеспечено выбором схютветствующего порога, а схожий тренд изменения пористости.

Рис. 7. Сравнение профилей пористости, полученных с помощью метода анализа гистограмм и метода установки порога

Высокая амплитуда вариации пористости, полученной методом установки порога, является следствием того, что анализировалась лишь небольшая область сечения, поскольку размер области, с которой может работать СТАп, охраничен оперативной памятью. Напротив, при анализе гистограмм были использованы все имеющиеся сечения (при анализе гистограмм использовались более 3 млн пикселей, а для метода установки порога на том же компьютере могло быть использовано только 360 тысяч пикселей).

5. Заключение

Разработан новый метод анализа рентгеновской микротомографии пористых сред для расчета профилей пористости, а также объемных долей слагающих материалов. Метод основан на анализе формы гистограмм отдельных сечений трехмерного изображения рентгеновской компьютерной микротомографии и, в отличие от стандартного метода установления порога, позволяет рассчитать долю различных веществ в образце не прибегая к заданию порога apriori.

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

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

Литература

1. Al-Abduwani F.A.H. fet. al]. Filtration of micron-sized particles in granular media revealed by x-ray computed tomography /7 Review of scientific instruments. 2005. V. 76. art. no.103704.

2. De. Zwart. A.H. Experiment and theoretical investigation of clogging processes near production wells using X-ray Tomography /7 SPE. 2008. 116411.

3. Jamankulov T. Experimental investigation of the conformance of deep bed filtration using computed tomography technique MSe thesis, Faculty of Geoteehnology. Delft University of Technology. Delft. 2005.

4. Saraf A. [et. al.J. Analysis of the effect of residual oil and particle trapping during produced water reinjection using X-ray tomography // SPE. — 2009. — 122137.

5. Alt M. A. J. [et. al.J. Measurement of the particle deposition profile in deep-bed filtration during produced water reinjection // SPE. — 2005. — 93056.

6. Cole M.E. Jet. al.J. Characterization of reservoir core using computed microtomographv // SPE. - 1996. - 30542.

7. Arns C.H. Jet. al.J. Pore scale characterization of carbonates using X-ray Microtomographv 11 SPE. - 2004. - 90368.

8. Gonzalez R.C., Woods R.E. Digital Image Processing //Addison-Wesley Pub. — 1992. — 730 p.

9. Gupta L., Sortrakul T. A Gaussian-mixture-based image segmentation algorithm // Patter Recognition. - 1998. - V. 21. - 315-325.

10. Kobayashi Y. [et. al.J. Evaluation of porosity in geomaterials treated with biogrout considering partial volume effect // Advances in computed tomography for geomaterials (GeoX). - 2010. - P. 289-294.

11. Kuwahara M. Processing of ri-angiocardiographic images // Digital Processing of Biomedical Images. — 1976. — P. 187-203.

12. Ives K.D. Rapid filtration review paper // Water. Res. — 1967. — N 4. — P. 201.

13. Iwasaki T. Some notes on sand filtration //J. Am. Water Works Ass. — 1937. — V. 29, N 10. - P. 1591-1602.

14. Herzig J.P., Leclerc D.M., Le Goff P. Flow of suspensions through porous media — application to deep filtration // Industrial and Engineering Chemistry. — 1970. — V. 62, N 5. - P. 8-35.

15. Шехтман Ю.М. Фильтрация малоконцентрированных суспензий. - М.: Недра, 1961. - 237 с.

16. Santos A., Bedrikovetsky P.G. Stochastic model for particulate suspension flow in porous media // Transport in Porous Media. — 2006. — V. 62. - P. 23-53.

Поступила в редакцию 25.06.2012