О пространственной интерполяции температуры грунта по данным термометрии скважин

Введение

Для неасимптотического численного моделирования нестационарной теплопроводности нужно знать начальные условия. Если теплопроводящей средой является грунт, то обычно используют данные термометрии скважин (таблично заданные зависимости температуры от глубины) для начального момента времени. Затем можно либо применить интерполяцию по рассеянным данным [1], либо решить стационарную задачу теплопроводности с граничными условиями первого рода в точках измерения температуры. Во втором случае, если задача нелинейная, могут возникнуть проблемы со сходимостью метода.

Интерполяция температуры грунта по данным термометрии скважин

Рисунок 1 – Интерполяция температуры грунта в программе Frost 3D Universal

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

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

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

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

  4. Вдоль термометрических скважин интерполяция сохраняет форму (монотонность и выпуклость) данных [3].

Постановка задачи

Пусть \Pi=[x_0,x_1]\times[y_0,y_1]\neq\emptyset и h:\Pi\rightarrow\mathbb{R} – поле высот, определяющее часть \left\{\left.\left(A,h(A)\right)\right|A\in\Pi\right\} поверхности \partial S\subset\mathbb{R}^3 бесконечно протяженного и бесконечно глубокого грунта S\subset\mathbb{R}^3.

Исследуем задачу интерполяции неизвестного температурного поля T:S\rightarrow\mathrm{Conv}\left\{\left.T_{ij}\right|j=1,2,\ldots,m_i;i=1,2,\ldots,n\right\} (\mathrm{Conv}\,P обозначает выпуклую оболочку множества точек P) по данным термометрии

    \[\left[\begin{matrix}T\left(A_i,z_{i1}\right) \\ T\left(A_i,z_{i2}\right) \\ \vdots \\ T\left(A_i,z_{im_i}\right) \end{matrix}\right]=\left[\begin{matrix}T_{i1} \\ T_{i2} \\ \vdots \\ T_{im_i} \end{matrix}\right]\!\!,~i=1,2,\ldots,n,\]

где z_{i1},z_{i2},\ldots,z_{im_i} – высоты точек измерения \left(A_i,z_{i1}\right),\left(A_i,z_{i2}\right),\ldots,\left(A_i,z_{im_i}\right) в i-ой термометрической скважине W_i=\left\{A_i\right\}\times\left(-\infty,\hbar(A_i)\right], A_i\in\mathbb{R}^2, i=1,2,\ldots,n. Здесь \hbar:\mathbb{R}^2\rightarrow\mathrm{Conv}\left(h(\Pi)\right) – «разумное» продолжение (см. стадию 3 ниже) функции h:\Pi\rightarrow\mathbb{R} на всю плоскость \mathbb{R}^2. Точные определения бесконечно протяженного и глубокого грунта S и его поверхности \partial S – следующие:

    \[S=\bigcup\limits_{A\in\mathbb{R}^2}\{A\}\times\left(-\infty,\hbar(A)\right],~\partial S=\left\{\left.\left(A,h(A)\right)\right|A\in\mathbb{R}^2\right\}.\]

Заметим, что \bigcup\limits_{i=1}^nW_i\subset S.

Алгоритм

С поставленной интерполяционной задачей можно справиться в несколько стадий. Перечислим их.

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

  2. Слияние дублирующихся точек измерения с заменой нескольких значений температуры их средним арифметическим. Перенумерация обновленного набора точек измерения из W_i и переопределение m_i, i=1,2,\ldots,n.

  3. Построение функции \hbar:\mathbb{R}^2\rightarrow\mathrm{Conv}\left(h(\Pi)\right) такой, что для всех точек A\in\mathbb{R}^2 имеем

        \[\hbar(A)=h(A_\Pi),\]

    где A_\Pi – ближайшая к A\in\mathbb{R}^2 точка \Pi. В частности, если A\in\Pi, то A_\Pi=A и \hbar(A)=h(A).

  4. Построение поля глубин d:S\rightarrow\left[0,\infty\right),

        \[d(A,z)\equiv\hbar(A)-z.\]

    Очевидно, d\left(\partial S\right)=\{0\}.

  5. Для каждого непустого 2-мерного слоя S_z=\left\{\left.A\in\mathbb{R}^2\right|(A,z)\in S\right\} (\left(\mathbb{R}^2\times\{z\}\right)\cap S\neq\emptyset) грунта получение сглаженной (с пропорциональной глубине интенсивностью) версии \tilde{\hbar}_z:S_z\rightarrow\mathrm{Conv}(h(\Pi)),

        \[\tilde{\hbar}_z=B_{\lambda d(A,z)}\hbar\]

    продолженного поля высот \hbar. Здесь

        \[\left(B_{\sigma}f\right)(A)\equiv\begin{cases} \frac{1}{2\pi\sigma^2}\displaystyle{\int\limits_{\mathbb{R}^2}}e^{-\frac{\left(|X-A|/\sigma\right)^2}{2}}f(X)dX, & \sigma>0; \\ f(A), & \sigma=0;\end{cases}\]

    B_\sigma – размытие по Гауссу [2] со среднеквадратическим отклонением (радиусом) \sigma\geq0. Используется

        \[\sigma=\lambda d\left(A,z\right),~\left(A,z\right)\in S,\]

    где физически безразмерная константа \lambda>0 не зависит от \left(A,z\right) и подобрана из эмпирических соображений.

    Заметим, что для всех \left(A,z\right)\in\partial S справедливо

        \[\tilde{\hbar}_z(A)=(B_0\hbar)(A)=\hbar(A).\]

  6. Построение сглаженного поля глубин \tilde{d}:S\rightarrow\mathbb{R},

        \[\tilde{d}\left(A,z\right)\equiv\hbar_z(A)-z.\]

    Для всех \left(A,z\right)\in\partial S имеем

        \[\tilde{d}\left(A,z\right)=\tilde{\hbar}_z(A)-\hbar(A)=\hbar(A)-\hbar(A)=0,\]

    т. е. \tilde{d}\left(\partial S\right)=\{0\}, как и для d. В нужных подобластях – коррекция поля \tilde{d} так, чтобы имело место строгое убывание \tilde d\left(A,z\right) по z\in\left(-\infty,\hbar(A)\right] для всех A\in\mathbb{R}^2, но, как и ранее, было верно \tilde{d}\left(\partial S\right)=\{0\}. Отсюда мы получаем, что \tilde{d}\geq0.

  7. Вычисление сглаженных глубин

        \[\tilde{d}_{ij}=\tilde{d}\left(A_i,z_{ij}\right)\geq0,~j=1,2,\ldots,m_i,~i=1,2,\ldots,n\]

    точек измерения температуры. В рамках одной термометрической скважины W_i эти глубины попарно различны, благодаря строгой монотонности \tilde{d}\left(A_i,z\right) по z\in\left(-\infty,\hbar(A_i)\right], i=1,2,\ldots,n.

    Построение линейных интерполяционных (для экстраполяции применяется метод ближайшего соседа) сплайнов T_i:\left[0,\infty\right)\rightarrow\mathrm{Conv}\left\{T_{i1},T_{i2,},\ldots,T_{im_i}\right\} таких, что

        \[\left[\begin{matrix}T_i\left(\tilde{d}_{i1}) \\ T_i\left(\tilde{d}_{i2}) \\ \vdots \\ T_i\left(\tilde{d}_{im_i}) \end{matrix}\right]=\left[\begin{matrix} T_{i1} \\ T_{i2} \\ \vdots \\ T_{im_i} \end{matrix}\right]\!\!,~i=1,2,\ldots,n.\]

  8. Если n=1, то полагаем

        \[T\left(A,z\right)\equiv T_1\left(\tilde{d}\left(A,z\right)\right)\]

    и завершаем работу.

  9. Построение выпуклой оболочки C=\mathrm{Conv}\left\{A_1,A_2,\ldots,A_n\right\} и триангуляции Делоне того же множества точек \left\{A_1,A_2,\ldots,A_n\right\}.

  10. Вычисление в каждой фиксированной точке \left(A^*,z^*\right)\in S (обозначим ее сглаженную глубину \tilde{d}\left(A^*,z^*\right) через \tilde{d}^*) искомой температуры T^* по следующему алгоритму. Пусть A_C^* – ближайшая к A^* точка выпуклой оболочки C. В частности, если A^*\in C, то A_C^*=A^*. Если точка A_C^* лежит в некотором невырожденном треугольнике Делоне A_iA_jA_k, i,j,k\in\left\{1,2,\ldots,n\right\} и поэтому существуют однозначно определяемые скаляры (барицентрические координаты) \alpha,\beta,\gamma\geq0 такие, что A_C^*=\alpha A_i+\beta A_j+\gamma A_k и \alpha+\beta+\gamma=1, то полагаем

        \[T^*=\alpha T_i\left(\tilde{d^*}\right)+\beta T_j\left(\tilde{d^*}\right)+\gamma T_k\left(\tilde{d^*}\right).\]

    Возможен также случай, когда триангуляция Делоне пуста (не содержит ни одного невырожденного треугольника). Тогда, учитывая, что n\geq2, заключаем, что точка A_C^* лежит на некотором невырожденном отрезке A_iA_j (i,j\in\{1,2,\ldots,n\}), т. е. A_C^*=\alpha A_i+\beta A_j, \alpha+\beta=1, где выбор скаляров \alpha,\beta\geq0 однозначен. Полагаем

        \[T^*=\alpha T_i\left(\tilde{d^*}\right)+\beta T_j\left(\tilde{d^*}\right).\]

Как видно из построения, изоповерхности результирующего температурного поля T адаптируются под рельеф грунта, но с ростом глубины d эта адаптация монотонно ослабевает.

Список использованных источников

  1. I. Amidror. Scattered data interpolation methods for electronic imaging systems: a survey. J. Electron. Imaging, 11(2), 157–176 (апрель 2002).

  2. Авторы Википедии. Gaussian blur. Википедия, свободная энциклопедия, https://en.wikipedia.org/wiki/Gaussian_blur (дата обращения: 3 августа, 2015).

  3. F. N. Fritsch and R. E. Carlson. Monotone Piecewise Cubic Interpolation. SIAM J. Numer. Anal., 17, 238–246 (1980).

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *


*

Можно использовать следующие HTML-теги и атрибуты: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>