Известия РАН. Теория и системы управления, 2019, № 6, стр. 103-118

ВЫДЕЛЕНИЕ И ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЦИЛИНДРИЧЕСКИХ И СФЕРИЧЕСКИХ ПОВЕРХНОСТЕЙ ПРИ СЕГМЕНТАЦИИ ДАЛЬНОСТНЫХ ЛАЗЕРНО-ЛОКАЦИОННЫХ ИЗОБРАЖЕНИЙ

В. М. Лисицын a*, К. В. Обросов a, В. А. Сафонов a

a ФГУП ГосНИИАС
Москва, Россия

* E-mail: lvm@gosniias.ru

Поступила в редакцию 28.05.2019
После доработки 25.06.2019
Принята к публикации 22.07.2019

Полный текст (PDF)

Аннотация

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

Введение. Эффективным средством дистанционного зондирования являются лазерно-локационные (ЛЛ) системы [1, 2]. Они представляют собой результат естественного слияния технологий радиолокационных и оптико-электронных датчиков и сочетают в себе радиолокационные методы обработки и анализа сигналов с оптическими методами формирования изображений. В зависимости от назначения и конструкции ЛЛ-системы могут обеспечивать при зондировании измерение дальности, скорости, степени деполяризации излучения и т.д. При этом достоверность измерений вероятностным образом связана с отношением сигнал/шум в тракте детектирования отраженного сигнала, которое также может быть измерено известными в радиолокации методами [3].

Особенности формируемой ЛЛ-информации обеспечивают возможность автоматической обработки при решении ряда задач в различных областях применения. С точки зрения систем технического зрения наибольший интерес представляют ЛЛ-датчики, формирующие дальностные 3-D-изображения. В перспективе роботизированные системы на основе информации, формируемой ЛЛ-датчиками, смогут анализировать наблюдаемые сцены, ориентироваться и перемещаться, распознавать образы объектов, принимать решения и т.д. Необходимым этапом обработки изображений для решения задач автоматического распознавания образов и анализа сцен является сегментация [49]. Сегментация – это разбиение изображения на непересекающиеся однородные по некоторому признаку области, которые имеют смысловую суть. Как правило, предполагается, что области соответствуют образам объектов или их частям, а границы областей – границам образов объектов. Сегментация изображения является первым шагом на пути автоматического анализа и “понимания” изображения. В частности, если говорить о дальностном ЛЛ-изображении, то сегментация является инструментом построения объемной модели сцены в виде совокупности фрагментов простейших поверхностей (плоскости, цилиндры, сферы, конуса и пр.). Сегментация позволяет перейти от дальностного ЛЛ-изображения (от совокупности не связанных между собой пространственно разнесенных замеров дальности) к геометрической модели сцены в виде совокупности непересекающихся областей с геометрически связанными между собой замерами дальности в пределах каждой области, т.е. к семантическому описанию форм на наблюдаемой сцене.

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

В [10] был предложен, а в [1113] развит метод выделения плоских поверхностей при сегментации доплеровских и дальностных ЛЛ-изображений. Метод состоит из двух этапов: поиска и последующего наращивания областей, содержащих признаки плоских поверхностей. Поиск плоскостей базируется на аппроксимации фрагмента сцены, попадающего в “скользящее окно”, методом наименьших квадратов исходя из гипотезы, что это фрагмент плоскости, и применении критерия согласия, основанного на распределении χ-квадрат, для принятия решения о верности данной гипотезы. Наращивание областей осуществляется итерационно с помощью критерия максимума апостериорной вероятности и использования марковской модели ЛЛ-изображения. Метод показал высокую эффективность при обработке реальных ЛЛ-изображений.

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

1. Постановка задачи. ЛЛ-система зондирует сцену под настильными углами, формируя растровое дальностное изображение, на которой расположен цилиндрический объект радиуса R с известной ориентацией оси. Центр объекта удален от ЛЛ-системы на расстояние Dc. Будем полагать, что ${{D}_{c}} \gg R$, зондирование производится в малом телесном угле и измерения дальности преобразованы из центральной проекции в Декартову систему координат (СК) OXYZ, начало которой совмещено с ЛЛ-системой, а ось OZ параллельна оси цилиндра (ось OY направлена на объект, ось OX перпендикулярна плоскости OYZ). Без потери общности для краткости будем предполагать, что ось цилиндра вертикальна, а ось OZ направлена вверх. Зондирования производятся только по наблюдаемой части цилиндрического объекта. В соответствии с [3] ошибки измерения дальности ξ имеют нормальное распределение с нулевым математическим ожиданием, а среднее квадратическое отклонение (СКО) оценок времени приема отраженного сигнала τс и дальности по i-му измерению определяются выражениями

(1.1)
${{\sigma }_{i}}({{\xi }_{\tau }}) = \frac{{{{\tau }_{с}}}}{{{{q}_{i}}\sqrt \pi }};\quad {{\sigma }_{i}}(\xi ) = {{\sigma }_{{\xi i}}} = \frac{{C{{\tau }_{с}}}}{{2{{q}_{i}}\sqrt \pi }},$
где τс – длительность импульса; С – скорость света; qi – отношение сигнал/шум. Причем $q_{i}^{2}$ = = S/N. Здесь S – энергия принимаемого сигнала; N – спектральная плотность мощности шума. В силу ${{D}_{c}} \gg R$ будем считать, что направления зондирования коллинеарны оси Y. Рассмотрим случай, когда эффективная ширина диаграммы направленности лазерного луча на местности относительно невелика, что характерно для лазерных локаторов с углом расходимости лазерного излучения менее 1 угл. мин на расстояниях 1–2 км. При этом СКО центра отражения от оси диаграммы направленности составляет единицы сантиметров, что существенно меньше СКО ошибок измерения дальности до отражающей поверхности. Поэтому в дальнейшем будем учитывать только ошибки измерения дальности.

Рассмотрим проекцию цилиндрической поверхности, ее оси, точек зондирования и измерений на горизонтальную плоскость (рис. 1). Примем гипотезу равноточности измерений дальности до всех точек цилиндрической поверхности. Направление зондирования совпадает с осью OY. Очевидно, что координаты точек зондирования Xi,Yi, i = $\overline {1,n} $, расположены на дуге окружности с центром, имеющим координаты Xс,Yс в соответствии с выражением:

(1.2)
${{({{Y}_{i}} - {{Y}_{c}})}^{2}} + {{({{X}_{i}} - {{X}_{c}})}^{2}} = R_{{}}^{2}.$
Рис. 1.

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

Измеряются координаты точек зондирования xi = Xi, yi = Yi + ξi c СКО σξ. Таким образом, i-м измерением в СК (XOY) является вектор (xi, yi), возникающий при координатах точки зондирования Xi, Yi.

2. Определение радиуса и координат оси цилиндрической поверхности. Рассмотрим результат однократного (i-го) зондирования цилиндрической поверхности как косвенное измерение ri радиуса R цилиндрической поверхности, имеющей координаты оси Xc, Yc на горизонтальной плоскости.

Отметим, что решение задачи оценки данных координат при традиционном способе оценки близости измерений к дуге окружности, т.е. путем расчета суммы квадратов отклонений самих косвенных измерений ri от оценки радиуса, приводит к появлению системы нелинейных уравнений, не имеющей аналитического решения. Поэтому при каждом зондировании будем оценивать косвенные измерения квадрата радиуса $r_{i}^{2}$ вместо ri. При координатах точки зондирования Xi,Yi прямое измерение производится с ошибкой ξi, т.е. xi = Xi, yi = Yi + ξi . Тогда косвенное измерение (ri)2 квадрата радиуса R2 имеет вид

(2.1)
$r_{{}}^{2} = {{\left( {{{X}_{c}} - {{X}_{i}}} \right)}^{2}} + {{\left( {{{Y}_{c}} - {{Y}_{i}} - {{\xi }_{i}}} \right)}^{2}}.$

Будем рассматривать i-е косвенное измерение (ri)2 квадрата радиуса как реализацию случайной величины r2. Учитывая (1.2), легко показать, что математическое ожидание случайной величины r2 при известных координатах Xс, Yс определяется выражением

(2.2)
$M(r_{{}}^{2}) = {{R}^{2}} + {{D}_{\xi }},$
где Dξ – дисперсия ошибки ξ.

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

(2.3)
$D(r_{i}^{2}) = M[{{(r_{i}^{2} - M(r_{i}^{2}))}^{2}}] = {{D}_{\xi }}[3{{D}_{\xi }} + 4{{\left( {{{Y}_{c}} - {{Y}_{i}}} \right)}^{2}}].$

Из (2.3) следует, что даже если принята гипотеза равноточности измерений дальности до всех точек цилиндрической поверхности, косвенные измерения квадрата радиуса не являются равноточными, поскольку зависят от Yi и неизвестной величины Yc. Это, вообще говоря, не позволяет для оценки квадрата радиуса использовать среднее арифметическое его измерений в различных точках зондирования. Однако, учитывая простоту вычисления среднего арифметического, целесообразно его использовать для оценки математического ожидания вспомогательной величины $R_{c}^{2}$, с помощью которой можно рассчитать искомую оценку радиуса R на основе процедуры, описанной ниже.

При известных координатах оси цилиндра (т.е. при рассмотрении их в качестве оптимизируемых параметров) выражение для оценки дисперсии случайной величины r2 можно использовать в качестве меры близости измерений xi, yi, i = $\overline {1,n} $, к дуге окружности:

(2.4)
$D(r_{{}}^{2}) = \frac{{Q({{x}_{1}},...,{{x}_{n}},{{y}_{1}},...,{{y}_{n}}\left| {{{x}_{c}},} \right.{{y}_{c}})}}{{n - 1}},$
где

$Q({{x}_{1}},...,{{x}_{n}},{{y}_{1}},...,{{y}_{n}}\left| {{{x}_{c}},} \right.{{y}_{c}}) = \sum\limits_{j = 1}^n {{{{\left[ {{{{\left( {{{x}_{c}} - {{x}_{j}}} \right)}}^{2}} + {{{\left( {{{y}_{c}} - {{y}_{j}}} \right)}}^{2}} - \frac{1}{n}\sum\limits_{j = 1}^n {[{{{\left( {{{x}_{c}} - {{x}_{i}}} \right)}}^{2}} + {{{\left( {{{y}_{c}} - {{y}_{i}}} \right)}}^{2}}]} } \right]}}^{2}}} .$

Для оценки координат центра окружности Xc, Yc определим оптимальные значения параметров xc, yc, минимизируя выражение (2.4), т.е. приравняем к нулю частные производные по xc и yc:

$\frac{{\partial Q}}{{\partial {{x}_{c}}}} = 0;\quad \frac{{\partial Q}}{{\partial {{y}_{c}}}} = 0,$
и решим получившуюся систему уравнений.

После преобразований выражение $\frac{{\partial Q}}{{\partial {{x}_{c}}}}$ примет вид

$\frac{{\partial Q}}{{\partial {{x}_{c}}}} = 4\left\{ {2{{x}_{c}}\sum\limits_{j = 1}^n {{{{(\bar {x} - {{x}_{j}})}}^{2}}} + 2{{y}_{c}}\sum\limits_{j = 1}^n {(\bar {y} - {{y}_{j}})(\bar {x} - {{x}_{j}})} + \sum\limits_{j = 1}^n {(\bar {x} - {{x}_{j}})\left[ {x_{j}^{2} + y_{j}^{2} - \frac{1}{n}\sum\limits_{j = 1}^n {(x_{i}^{2} + y_{i}^{2})} } \right]} } \right\},$
здесь $\bar {x} = \frac{1}{n}\sum\limits_{i = 1}^n {{{x}_{i}}} $, $\bar {y} = \frac{1}{n}\sum\limits_{i = 1}^n {{{y}_{i}}} $.

Аналогичные преобразования можно провести для $\frac{{\partial Q}}{{\partial {{y}_{c}}}}$. Приравняв к нулю частные производные, можно получить следующую систему линейных уравнений:

(2.5)
$\left\{ {\begin{array}{*{20}{c}} {{{a}_{{11}}}{{x}_{c}} + {{a}_{{12}}}{{y}_{c}} = {{b}_{1}},} \\ {{{a}_{{21}}}{{x}_{c}} + {{a}_{{22}}}{{y}_{c}} = {{b}_{2}},} \end{array}} \right.$
где

${{a}_{{11}}} = 2\sum\limits_{i = 1}^n {{{{\left( {\bar {x} - {{x}_{i}}} \right)}}^{2}}} ,\quad {{a}_{{12}}} = 2\sum\limits_{i = 1}^n {\left( {\bar {y} - {{y}_{i}}} \right)\left( {\bar {x} - {{x}_{i}}} \right)} ,\quad {{a}_{{21}}} = {{a}_{{12}}},\quad {{a}_{{22}}} = 2\sum\limits_{i = 1}^n {{{{\left( {\bar {y} - {{y}_{i}}} \right)}}^{2}}} ,$
${{b}_{1}} = \sum\limits_{j = 1}^n {\left( {\bar {x} - {{x}_{j}}} \right)} \left[ {\frac{1}{n}\sum\limits_{i = 1}^n {(x_{i}^{2} + y_{i}^{2})} - x_{j}^{2} - y_{j}^{2}} \right] = - \sum\limits_{j = 1}^n {\left( {\bar {x} - {{x}_{j}}} \right)(x_{j}^{2} + y_{j}^{2})} ,$
${{b}_{2}} = \sum\limits_{j = 1}^n {\left( {\bar {y} - {{y}_{j}}} \right)\left[ {\frac{1}{n}\sum\limits_{i = 1}^n {(x_{i}^{2} + y_{i}^{2})} - x_{j}^{2} - y_{j}^{2}} \right]} = - \sum\limits_{j = 1}^n {\left( {\bar {y} - {{y}_{j}}} \right)(x_{j}^{2} + y_{j}^{2})} .$

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

(2.6)
${{x}_{c}} = \frac{{{{\Delta }_{1}}}}{\Delta },\quad {{y}_{c}} = \frac{{{{\Delta }_{2}}}}{\Delta },$
где Δ1 = b1a22b2a12, Δ2 = b2a11b1a21, Δ = a11a22a21a12.

Можно показать, что полученная оценка положения центра окружности является смещенной. Однако величину смещения центра в свою очередь можно оценить. Чтобы получить аналитически оценки смещения центра окружности в условиях наличия ошибок измерения, рассмотрим асимптотический случай, когда количество измерений n → ∞. При этом в соответствии с законом больших чисел оценки математического ожидания и дисперсии оцениваемых величин, полученные по выборке измерений, приближаются по вероятности к соответствующим истинным статистическим характеристикам.

Без потери общности будем считать, что ось OY проходит через центр окружности. В выборке измерений (x1, …, xn, …, y1, … ,yn, …) при n → ∞ будем рассматривать координату xi как реализацию случайной величины x, равномерно распределенной в диапазоне [–R, R], где R – радиус окружности (общий случай, когда наблюдается или зондируется лишь часть дуги, будет рассмотрен далее). Величина yi в соответствии с принятым нормальным законом распределения может быть выражена как yi = Yi + ξi, где Yi – координаты точки зондирования, лежащей на окружности при i-м измерении, ξ – нормально распределенная случайная ошибка с характеристиками M(ξ(x)) = 0, $D(\xi (x)) = \sigma _{\xi }^{2}$ для ∀x; Yi определяется очевидным геометрическим соотношением Yi = = ${{D}_{c}} - \sqrt {R_{{}}^{2} - x_{i}^{2}} $, где Dc – дальность до центра окружности, отсчитываемая вдоль оси OY.

Допустим, что оценка положения центра окружности сделана с некоторой ошибкой, и предполагаемый центр смещен по оси OY относительно истинного на некоторую величину Δ (рис. 2). В силу симметрии смещение по оси OX стремится к нулю при увеличении количества измерений.

Рис. 2.

Связь величин Rc, радиуса R, измерения (xi, yi), величины смещения Δ и ошибки измерения ξi

Обозначим через Rc отрезок, связывающий полученное измерение (xi, yi) с найденным смещенным центром окружности. Введем отрезок RΔ, связывающий точку зондирования Yi, т.е. точку, находящуюся на рассматриваемой окружности с координатами (xi, Yi), со смещенным центром (xc, yc) с координатами (Xc, Yc – Δ). Выразим Rc через значение R, измерение (xi, yi), величины смещения Δ и ошибки измерения ξi. Для i-го измерения имеют место следующие соотношения:

$\begin{gathered} R_{c}^{2} = {{(\sqrt {R_{\Delta }^{2} - x_{i}^{2}} - {{\xi }_{i}})}^{2}} + x_{i}^{2} = R_{\Delta }^{2} - 2{{\xi }_{i}}\sqrt {R_{\Delta }^{2} - x_{i}^{2}} + \xi _{i}^{2} = \\ \, = R_{{}}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{i}^{2}} + {{\Delta }^{2}} - 2{{\xi }_{i}}\sqrt {R_{{}}^{2} - x_{i}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{i}^{2}} + {{\Delta }^{2}}} + \xi _{i}^{2}; \\ \end{gathered} $
где ${{R}_{\Delta }} = \sqrt {{{{(\sqrt {R_{{}}^{2} - x_{i}^{2}} - \Delta )}}^{2}} + x_{i}^{2}} = \sqrt {R_{{}}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{i}^{2}} + {{\Delta }^{2}}} .$

Используя всю выборку измерений, оценим математическое ожидание и дисперсию величины $R_{c}^{2}$ для некоторого значения Δ на основе полученных соотношений и сделанных допущений:

$M(R_{c}^{2}) = R_{{}}^{2} - 2\Delta M(\sqrt {R_{{}}^{2} - x_{{}}^{2}} ) + {{\Delta }^{2}} - 2M\{ \xi \sqrt {R_{{}}^{2} - x_{{}}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{{}}^{2}} + {{\Delta }^{2}}} \} + M(\xi _{{}}^{2}).$

Из того, что случайные величины ξ и x независимы, а M(ξ) = 0, следует

$M\{ \xi \sqrt {R_{{}}^{2} - x_{{}}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{{}}^{2}} + {{\Delta }^{2}}} \} = M(\xi )M\{ \sqrt {R_{{}}^{2} - x_{{}}^{2} - 2\Delta \sqrt {R_{{}}^{2} - x_{{}}^{2}} + {{\Delta }^{2}}} \} = 0.$

Таким образом,

(2.7)
$M(R_{c}^{2}) = R_{{}}^{2} - 2\Delta M(\sqrt {R_{{}}^{2} - {{x}^{2}}} ) + {{\Delta }^{2}} + M({{\xi }^{2}});$
$M{{[R_{c}^{2} - M(R_{c}^{2})]}^{2}}\, = \,M{{\{ - 2\Delta [\sqrt {R_{{}}^{2} - {{x}^{2}}} \, - \,M(\sqrt {R_{{}}^{2} - {{x}^{2}}} )]\, - \,2\xi \sqrt {R_{{}}^{2} - {{x}^{2}} - 2\Delta \sqrt {R_{{}}^{2} - {{x}^{2}}} + {{\Delta }^{2}}} + {{\xi }^{2}}\, - \,M({{\xi }^{2}})\} }^{2}}\, = $
$ = \,\;4{{\Delta }^{2}}\{ R_{{}}^{2}\, - \,M({{x}^{2}}) - {{[M(\sqrt {R_{{}}^{2} - {{x}^{2}}} )]}^{2}}\} \, + \,4M({{\xi }^{2}})[R_{{}}^{2} - M({{x}^{2}})\, - \,2\Delta M(\sqrt {R_{{}}^{2} - {{x}^{2}}} )\, + {{\Delta }^{2}}]\, + \,M{{[{{\xi }^{2}} - M({{\xi }^{2}})]}^{2}}.$

Далее для случайной величины x, распределенной равномерно в интервале [–R, R], имеем

(2.8)
$M({{x}^{2}}) = \frac{1}{{2R}}\int\limits_{ - R}^R {{{x}^{2}}dx} = \frac{1}{{2R}}\left. {\frac{{{{x}^{3}}}}{3}} \right|_{{ - R}}^{R} = \frac{{R_{{}}^{2}}}{3}.$

Величина $M(\sqrt {R_{{}}^{2} - {{x}^{2}}} )$ вычисляется с помощью интегрирования по частям:

(2.9)
$M(\sqrt {R_{{}}^{2} - {{x}^{2}}} ) = \frac{1}{{2R}}\int\limits_{ - R}^R {\sqrt {R_{{}}^{2} - {{x}^{2}}} dx = \frac{1}{{4R}}\left. {\left( {R_{{}}^{2}\arcsin \frac{x}{R} + x\sqrt {R_{{}}^{2} - {{x}^{2}}} } \right)} \right|_{{ - R}}^{R} = \frac{{\pi R}}{4}} .$

Имея в виду, что $M({{\xi }^{2}}) = \sigma _{\xi }^{2},$ $M({{\xi }^{4}}) = 3\sigma _{\xi }^{4},$получаем

(2.10)
$\begin{gathered} M{{[R_{c}^{2} - M(R_{c}^{2})]}^{2}} = 4{{\Delta }^{2}}\left( {R_{{}}^{2} - \frac{{{{R}^{2}}}}{3} - \frac{{{{\pi }^{2}}{{R}^{2}}}}{{16}}} \right) + 4\sigma _{\xi }^{2}\left( {R_{{}}^{2} - \frac{{{{R}^{2}}}}{3} - 2\Delta \frac{{\pi R}}{4} + {{\Delta }^{2}}} \right) + 3\sigma _{\xi }^{4} - \sigma _{\xi }^{4} = \\ \, = 4{{\Delta }^{2}}{{R}^{2}}\left( {\frac{2}{3} - \frac{{{{\pi }^{2}}}}{{16}}} \right) + 4\sigma _{\xi }^{2}\left( {\frac{2}{3}{{R}^{2}} - \Delta \frac{{\pi R}}{2} + {{\Delta }^{2}}} \right) + 2\sigma _{\xi }^{4} = 4{{\Delta }^{2}}\left[ {{{R}^{2}}\left( {\frac{2}{3} - \frac{{{{\pi }^{2}}}}{{16}}} \right) + \sigma _{\xi }^{2}} \right] - \\ \, - 2\sigma _{\xi }^{2}\Delta \pi R + 2\sigma _{\xi }^{2}\left( {\frac{4}{3}{{R}^{2}} + \sigma _{\xi }^{2}} \right). \\ \end{gathered} $

Значение Δ, которое доставляет минимум выражению (2.10), соответствует математическому ожиданию $R_{c}^{2}$ для положения центра окружности, найденного с помощью (2.5), (2.6). Чтобы определить значение Δ, доставляющее минимум (2.10), проведем дифференцирование, и полученное выражение приравняем нулю:

$\frac{{\partial M{{{[R_{c}^{2} - M(R_{c}^{2})]}}^{2}}}}{{\partial \Delta }} = 8\Delta \left[ {{{R}^{2}}\left( {\frac{2}{3} - \frac{{{{\pi }^{2}}}}{{16}}} \right) + \sigma _{\xi }^{2}} \right] - 2\sigma _{\xi }^{2}\pi R = 0.$

Отсюда

(2.11)
$\Delta = \frac{{\sigma _{\xi }^{2}\frac{{\pi R}}{4}}}{{R_{{}}^{2}\left( {\frac{2}{3} - \frac{{{{\pi }^{2}}}}{{16}}} \right) + \sigma _{\xi }^{2}}}.$

Подставив найденное значение смещения центра (2.11) с учетом (2.8) в выражение (2.7), можно получить следующее уравнение относительно радиуса окружности R:

(2.12)
${{R}^{6}} + {{a}_{1}}{{R}^{4}} + {{a}_{2}}{{R}^{2}} + {{a}_{3}} = 0,$
где

(2.13)
$\begin{gathered} {{a}_{1}} = \frac{{\sigma _{\xi }^{2}}}{k}\left( {2 + k - 2{{{\left( {\frac{\pi }{4}} \right)}}^{2}}} \right) - M(R_{c}^{2});\quad k = \frac{2}{3} - \frac{{{{\pi }^{2}}}}{{16}}; \\ {{a}_{2}} = \sigma _{\xi }^{2}\left[ {\frac{{\sigma _{\xi }^{2}}}{{{{k}^{2}}}}\left( {1 + 2k - {{{\left( {\frac{\pi }{4}} \right)}}^{2}}} \right) - \frac{2}{k}M(R_{c}^{2})} \right]; \\ {{a}_{3}} = \frac{{\sigma _{\xi }^{4}}}{{{{k}^{2}}}}[\sigma _{\xi }^{2} - M(R_{c}^{2})]. \\ \end{gathered} $

В (2.13) вместо величины $M(R_{c}^{2})$ необходимо подставить ее оценку:

(2.14)
$M(R_{c}^{2}) = {{\bar {r}}^{2}} = \frac{1}{n}\sum\limits_{i = 1}^n {r_{i}^{2}} .$

Решив бикубическое уравнение (2.12) и взяв положительный корень (моделирование показало, что уравнение имеет только один действительный положительный корень), получаем оценку радиуса окружности $\widehat R$. В свою очередь, подставив оценку радиуса $\widehat R$ в (2.11), определяем смещение центра по оси Y. Увеличив на Δ значение yc, получаем несмещенную оценку координат центра оси цилиндрической поверхности Yc = yc + Δ.

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

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

2. Оценивается величина начального радиуса R1 (j = 1) с помощью выражения

(2.15)
${{\hat {R}}_{j}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {[{{{\left( {{{x}_{c}} - {{x}_{i}}} \right)}}^{2}} + {{{\left( {{{y}_{c}} - {{y}_{i}}} \right)}}^{2}}]} - \sigma _{\xi }^{2}} ,$
которое вытекает из (2.2).

3. Используя известное значение σξ и полученное значение радиуса Rj (j = $\overline {1,k} $, где k – количество итераций), по формуле (2.11) вычисляется величина смещения центра по координате, коллинеарной направлению измерений Δj.

4. Центр окружности смещается на величину Δj в сторону увеличения дальности в точку с координатами (xц, yц) = (xн, yн+ Δj).

5. С помощью (2.15) рассчитывается уточненное значение радиуса Rj + 1 относительно нового положения центра.

Далее процедура повторяется с п. 3. Моделирование показало, что процедура быстро сходится к истинным значениям искомых величин.

Приведенные выше выражения справедливы для случая, когда проекции точек зондирования распределены по всей полуокружности. Аналогичные выражения можно получить для случая, когда проекции точек зондирования находятся только на фрагменте полуокружности. Для построения алгоритма выделения цилиндрических поверхностей на реальных ЛЛ-изображениях целесообразно рассмотреть ситуацию, в которой для выделения цилиндрической поверхности используется лишь часть зондирований, проекции которых соответствуют некоторой дуге окружности, симметрично расположенной относительно оси цилиндра. При этом в представленных выше преобразованиях меняются пределы интегрирования в выражениях (2.8) и (2.9). Введем коэффициент m ∈ (q, 1], где q из практических соображений может быть выбрано как q ≥ 0.1…0.2.

Для случайной величины x, распределенной равномерно в интервале [–mR, mR], имеем

(2.16)
$M({{x}^{2}}) = \frac{1}{{2kR}}\int\limits_{ - mR}^{mR} {{{x}^{2}}dx} = \frac{1}{{2mR}}\left. {\frac{{{{x}^{3}}}}{3}} \right|_{{ - mR}}^{{mR}} = \frac{{{{m}^{2}}{{R}^{2}}}}{3}.$

Величина $M(\sqrt {{{R}^{2}} - {{x}^{2}}} )$ аналогично (2.9) вычисляется с помощью интегрирования по частям:

(2.17)
$\begin{gathered} M(\sqrt {{{R}^{2}} - {{x}^{2}}} ) = \frac{1}{{2mR}}\int\limits_{ - mR}^{mR} {\sqrt {{{R}^{2}} - {{x}^{2}}} dx} = \frac{1}{{4mR}}\left. {\left( {{{R}^{2}}\arcsin \frac{x}{R} + x\sqrt {{{R}^{2}} - {{x}^{2}}} } \right)} \right|_{{ - mR}}^{{mR}} = \\ \, = \frac{R}{{2m}}(\arcsin m + m\sqrt {1 - {{m}^{2}}} ). \\ \end{gathered} $

Подставив полученные выражения (2.16), (2.17) в (2.10), проведя дифференцирование и приравняв нулю результат, можно получить величину смещения для рассматриваемого случая наблюдения части дуги окружности:

(2.18)
$\Delta = \frac{{{{C}_{1}}R\sigma _{\xi }^{2}}}{{{{C}_{2}}{{R}^{2}} + \sigma _{\xi }^{2}}},$
где

${{C}_{1}} = \frac{1}{{2m}}[arcsinm + m\sqrt {1 - {{m}^{2}}} ],\quad {{C}_{2}} = 1 - \frac{{{{m}^{2}}}}{3} - C_{1}^{2}.$

Тогда уравнение (2.12) примет вид

(2.19)
${{a}_{1}}{{R}^{6}} + {{a}_{2}}{{R}^{4}} + {{a}_{3}}{{R}^{2}} + {{a}_{4}} = 0,$
где ${{a}_{1}} = C_{2}^{2};$

${{a}_{2}} = {{C}_{2}}[\sigma _{\xi }^{2}(2 - 2C_{1}^{2}\, + {{C}_{2}}) - {{C}_{2}}M(R_{с}^{2})];$
${{a}_{3}} = \sigma _{\xi }^{2}[\sigma _{\xi }^{2} - C_{1}^{2}\sigma _{\xi }^{2} + 2{{C}_{2}}\sigma _{\xi }^{2} - 2{{C}_{2}}M(R_{c}^{2})];$
${{a}_{4}} = \sigma _{\xi }^{4}[\sigma _{\xi }^{2} - M(R_{c}^{2})].$

Здесь $M(R_{c}^{2})$ оценивается с помощью (2.14).

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

В качестве меры близости измерений xi, yi ,zi, i = $\overline {1,n} $, к поверхности сферы аналогично критерию (2.6) будем использовать оценку дисперсии косвенных измерений квадрата радиуса сферы:

$\widehat D({{r}^{2}}) = \frac{{Q({{x}_{1}},...,{{x}_{n}},{{y}_{1}},...,{{y}_{n}},{{z}_{1}},...,{{z}_{n}}\left| {{{x}_{c}},} \right.{{y}_{c}},{{z}_{c}})}}{{n - 1}},$
где

(3.1)
$\begin{gathered} Q({{x}_{1}},...,{{x}_{n}},{{y}_{1}},...,{{y}_{n}},{{z}_{1}},...,{{z}_{n}}\left| {{{x}_{c}},} \right.{{y}_{c}},{{z}_{c}}) = \\ \, = \sum\limits_{j = 1}^n {{{{\left[ {{{{\left( {{{x}_{c}} - {{x}_{j}}} \right)}}^{2}} + {{{\left( {{{y}_{c}} - {{y}_{j}}} \right)}}^{2}} + {{{\left( {{{z}_{c}} - {{z}_{j}}} \right)}}^{2}} - \frac{1}{n}\sum\limits_{i = 1}^n {[{{{({{x}_{c}} - {{x}_{i}})}}^{2}} + {{{\left( {{{y}_{c}} - {{y}_{i}}} \right)}}^{2}} + {{{\left( {{{z}_{c}} - {{z}_{j}}} \right)}}^{2}}]} } \right]}}^{2}}.} \\ \end{gathered} $

Для оценки координат центра сферы Xc, Yc, Zc определим оптимальные значения параметров xc, yc, zc, минимизируя выражение (3.1), т.е. приравняем к нулю частные производные по xc, yc и zc:

$\frac{{\partial Q}}{{\partial {{x}_{c}}}} = 0,\quad \frac{{\partial Q}}{{\partial {{y}_{c}}}} = 0,\quad \frac{{\partial Q}}{{\partial {{z}_{c}}}} = 0.$

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

(3.2)
$\left\{ {\begin{array}{*{20}{c}} {{{a}_{{11}}}{{x}_{c}} + {{a}_{{12}}}{{y}_{c}} + {{a}_{{13}}}{{z}_{c}} = {{b}_{1}},} \\ {{{a}_{{21}}}{{x}_{c}} + {{a}_{{22}}}{{y}_{c}} + {{a}_{{23}}}{{z}_{c}} = {{b}_{2}},} \\ {{{a}_{{31}}}{{x}_{c}} + {{a}_{{32}}}{{y}_{c}} + {{a}_{{33}}}{{z}_{c}} = {{b}_{3}},} \end{array}} \right.$
где

${{a}_{{11}}} = 2\sum\limits_{i = 1}^n {{{{\left( {{{x}_{i}} - \bar {x}} \right)}}^{2}}} ;\quad {{a}_{{12}}} = 2\sum\limits_{i = 1}^n {\left( {\bar {y} - {{y}_{i}}} \right)\left( {\bar {x} - {{x}_{i}}} \right)} ;\quad {{a}_{{13}}} = 2\sum\limits_{i = 1}^n {\left( {\bar {z} - {{z}_{i}}} \right)\left( {\bar {x} - {{x}_{i}}} \right)} ;\quad {{a}_{{21}}} = {{a}_{{12}}};$
${{a}_{{22}}} = 2\sum\limits_{i = 1}^n {{{{\left( {{{y}_{i}} - \bar {y}} \right)}}^{2}}} ;\quad {{a}_{{23}}} = 2\sum\limits_{i = 1}^n {\left( {\bar {z} - {{z}_{i}}} \right)\left( {\bar {y} - {{y}_{i}}} \right)} ;\quad {{a}_{{31}}} = {{a}_{{13}}};\quad {{a}_{{32}}} = {{a}_{{23}}};$
${{a}_{{33}}} = 2\sum\limits_{i = 1}^n {{{{\left( {{{z}_{i}} - \bar {z}} \right)}}^{2}}} ;\quad {{b}_{1}} = - \sum\limits_{j = 1}^n {\left( {\bar {x} - {{x}_{j}}} \right)(x_{j}^{2} + y_{j}^{2} + z_{j}^{2})} ;\quad {{b}_{2}} = - \sum\limits_{j = 1}^n {\left( {\bar {y} - {{y}_{j}}} \right)(x_{j}^{2} + y_{j}^{2} + z_{j}^{2})} ;$
${{b}_{3}} = - \sum\limits_{j = 1}^n {\left( {\bar {z} - {{z}_{j}}} \right)(x_{j}^{2} + y_{j}^{2} + z_{j}^{2}){\kern 1pt} } .$

Здесь

$\bar {x} = \frac{{\sum\limits_{i = 1}^n {{{x}_{i}}} }}{n},\quad \bar {y} = \frac{{\sum\limits_{i = 1}^n {{{y}_{i}}} }}{n},\quad \bar {z} = \frac{{\sum\limits_{i = 1}^n {{{z}_{i}}} }}{n}.$

Система (3.2) может быть решена по правилу Крамера.

Можно показать, что данная оценка положения центра сферы является смещенной. Чтобы получить оценку смещения центра, проанализируем асимптотический случай, когда количество измерений n → ∞. Без потери общности будем считать, что ось OY проходит через центр сферы. В выборке измерений (x1, …, xn, …, y1, …, yn, …, z1, … zn, …) при n → ∞ будем рассматривать координаты xi,zi i-го измерения как реализацию независимых случайных величин x, z, равномерно распределенных на отрезке [–R, R], где R – радиус сферы (предлагаемый подход позволяет получить оценки положения центра сфера и ее радиус и для случая, когда зондируется лишь часть полусферы). Величина yi может быть выражена как yi = Yi + ξi, где Yi – координаты точки зондирования, лежащей на поверхности сферы при i-м измерении, ξ – нормально распределенная случайная ошибка с характеристиками: M(ξ(x)) = 0, $D(\xi (x)) = \sigma _{\xi }^{2}$ для ∀x; Yi определяется очевидным геометрическим соотношением ${{Y}_{i}} = {{D}_{c}} - \sqrt {R_{{}}^{2} - x_{i}^{2} - z_{i}^{2}} $, где Dc – дальность до центра сферы, отсчитываемая вдоль оси OY. Предположим, что оценка положения центра сделана с некоторой ошибкой и предполагаемый центр смещен по оси OY относительно истинного на некоторую величину Δ. В силу симметрии смещение по осям OX и OZ стремится к нулю при увеличении количества измерений.

Обозначим через Rc отрезок, связывающий полученное измерение (xi, yi, zi) с найденным смещенным центром сферы. Для i-го измерения имеют место следующие соотношения:

(3.3)
$\begin{gathered} R_{c}^{2} = {{(\sqrt {{{R}^{2}} - x_{i}^{2} - z_{i}^{2}} - \Delta - {{\xi }_{i}})}^{2}} + x_{i}^{2} + z_{i}^{2} = \\ \, = {{R}^{2}} - 2\Delta \sqrt {{{R}^{2}} - x_{i}^{2} - z_{i}^{2}} + {{\Delta }^{2}} - 2{{\xi }_{i}}(\sqrt {{{R}^{2}} - x_{i}^{2} - z_{i}^{2}} - \Delta ) + \xi _{i}^{2}. \\ \end{gathered} $

Из (3.3) следует

(3.4)
$M(R_{c}^{2}) = {{R}^{2}} - 2\Delta M(\sqrt {{{R}^{2}} - {{x}^{2}} - {{z}^{2}}} ) + {{\Delta }^{2}} + M({{\xi }^{2}}),$
(3.5)
$\begin{gathered} M{{[R_{c}^{2} - M(R_{c}^{2})]}^{2}} = 4{{\delta }^{2}}\{ {{R}^{2}} - M({{x}^{2}} + {{z}^{2}}) - {{[M(\sqrt {{{R}^{2}} - {{x}^{2}} - {{z}^{2}}} )]}^{2}}\} + \\ \, + 4M({{\xi }^{2}})[{{R}^{2}} - M({{x}^{2}} + {{z}^{2}}) - 2\Delta M(\sqrt {{{R}^{2}} - {{x}^{2}} - {{z}^{2}}} ) + {{\Delta }^{2}}] + M{{[{{\xi }^{2}} - M({{\xi }^{2}})]}^{2}}. \\ \end{gathered} $

Плотности распределения случайных величин x, z задаются выражениями $p\left( x \right) = p\left( z \right) = 1{\text{/2}}R$. Их совместная плотность распределения с учетом их независимости $p(x,z) = p(x)p(z) = 1{\text{/}}4{{R}^{2}}$. Очевидно, что ${{x}^{2}} + {{z}^{2}} \leqslant {{R}^{2}}$. Вероятность этого события $P({{x}^{2}} + {{z}^{2}} \leqslant {{R}^{2}}) = \pi {{R}^{2}}{\text{/}}4{{R}^{2}} = \pi {\text{/}}4$. Условная плотность вероятности pc(x, z) распределения координат x,z в круге ${{x}^{2}} + {{z}^{2}} \leqslant {{R}^{2}}$определяется выражением

${{p}_{c}}\left( {x,z} \right) = \frac{{p\left( {x,z} \right)}}{{P({{x}^{2}} + {{z}^{2}} \leqslant R_{{}}^{2})}} = \frac{{1{\text{/}}(4R_{{}}^{2})}}{{\pi {\text{/}}4}} = \frac{1}{{\pi R_{{}}^{2}}}.$

Тогда

$M({{x}^{2}} + {{z}^{2}}) = \iint\limits_{{{x}^{2}} + {{z}^{2}} \leqslant R_{{}}^{2}} {({{x}^{2}} + {{z}^{2}}){{p}_{c}}\left( {x,z} \right)dxdz} = \frac{1}{{\pi {{R}^{2}}}}\iint\limits_{{{x}^{2}} + {{z}^{2}} \leqslant R_{{}}^{2}} {({{x}^{2}} + {{z}^{2}})dxdz}.$

Путем перехода в полярную СК с якобианом преобразования J = ρ можно получить

$\iint\limits_{{{x}^{2}} + {{z}^{2}} \leqslant R_{{}}^{2}} {({{x}^{2}} + {{z}^{2}})dxdz} = \iint\limits_{\begin{array}{*{20}{c}} {0 \leqslant \rho \leqslant R} \\ {0 \leqslant \phi \leqslant 2\pi } \end{array}} {{{\rho }^{2}}Jd\varphi d\rho } = \int\limits_0^R {\left( {\int\limits_0^{2\pi } {{{\rho }^{3}}d\varphi } } \right)} d\rho = 2\pi \int\limits_0^R {{{\rho }^{3}}d\rho } = 2\pi \frac{{{{R}^{4}}}}{4} = \tfrac{{\pi {{R}^{4}}}}{2}.$

Тогда

$M({{x}^{2}} + {{z}^{2}}) = \frac{{\pi {{R}^{4}}{\text{/}}2}}{{\pi {{R}^{2}}}} = \frac{{{{R}^{2}}}}{2}.$

Аналогично

$\begin{gathered} M(\sqrt {{{R}^{2}} - {{x}^{2}} - {{z}^{2}}} ) = \iint\limits_{{{x}^{2}} + {{z}^{2}} \leqslant {{R}^{2}}} {\sqrt {{{R}^{2}} - {{x}^{2}} - {{z}^{2}}} {{p}_{c}}\left( {x,z} \right)dxdz} = \\ \, = \frac{1}{{\pi {{R}^{2}}}}\int\limits_0^R {\left( {\int\limits_0^{2\pi } {\rho \sqrt {{{R}^{2}} - {{\rho }^{2}}} d\varphi } } \right)} d\rho = \frac{2}{3}R. \\ \end{gathered} $

Подставляя найденные выражения для математических ожиданий в (3.5), можно получить

(3.6)
$M{{[R_{c}^{2} - M(R_{c}^{2})]}^{2}} = 4{{\Delta }^{2}}\left[ {\frac{{{{R}^{2}}}}{{18}} + \sigma _{\xi }^{2}} \right] - \frac{{16\sigma _{\xi }^{2}\Delta R}}{3} + 2\sigma _{\xi }^{2}({{R}^{2}} + \sigma _{\xi }^{2}).$

Как и в разд. 2, значение Δ, которое доставляет минимум выражению (3.6), соответствует математическому ожиданию $R_{c}^{2}$ для положения центра окружности, найденного с помощью (3.2). Чтобы определить значение Δ, доставляющее минимум (3.6), проведем дифференцирование и результат приравняем нулю:

(3.7)
$\frac{{\partial M{{{[R_{c}^{2} - M(R_{c}^{2})]}}^{2}}}}{{\partial \Delta }} = 8\Delta \left[ {\frac{{{{R}^{2}}}}{{18}} + \sigma _{\xi }^{2}} \right] - \frac{{16\sigma _{\xi }^{2}R}}{3} = 0.$

Отсюда

(3.8)
$\Delta = \frac{{\frac{2}{3}\sigma _{\xi }^{2}R}}{{\frac{{{{R}^{2}}}}{{18}} + \sigma _{\xi }^{2}}}.$

Подставив определенное с помощью (3.7) значение смещения центра (3.8) с учетом (3.6) в выражение (3.4), можно получить следующее уравнение относительно радиуса окружности R:

${{R}^{6}} + {{a}_{1}}{{R}^{4}} + {{a}_{2}}{{R}^{2}} + {{a}_{3}} = 0,$
где

(3.9)
$\begin{gathered} {{a}_{1}} = 21\sigma _{\xi }^{2} - M(\tilde {R}_{c}^{2}); \\ {{a}_{2}} = 36\sigma _{\xi }^{2}[6\sigma _{\xi }^{2} - M(\tilde {R}_{c}^{2})]; \\ {{a}_{3}} = 324\sigma _{\xi }^{4}[\sigma _{\xi }^{2} - M(\tilde {R}_{c}^{2})]. \\ \end{gathered} $

Решив бикубическое уравнение с коэффициентами (3.9) и взяв положительный корень, получаем оценку радиуса сферической поверхности $\widehat R$. В свою очередь, подставив оценку радиуса $\widehat R$ в (3.8), определяем смещение центра по оси OY. Увеличив на Δ значение yc, находим несмещенную оценку координат центра сферической поверхности.

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

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

4. Моделирование по синтезированным ЛЛ-изображениям. Для проверки математических положений было проведено моделирование по синтезированным дальностным изображениям. На рис. 3 представлена теоретическая зависимость величины относительного смещения оси цилиндра (Δтеор/R) от ее истинного положения, от относительной величины СКО измерения дальности (σξ/R). Значение смещения Δтеор получено с помощью (2.12) для заданного постоянного радиуса R и вариаций значения σξ.

Рис. 3.

Теоретическая зависимость величины относительного смещения оси цилиндра от ее истинного положения, от относительной величины СКО измерения дальности

В таблице приведены результаты моделирования для следующих исходных данных: положение оси цилиндра Xc = 0, Yc = 2000 м, радиус цилиндра R = 3 м, количество измерений в одном эксперименте n = 1000, величина СКО измерения дальности варьировалась от σξ = 0.2 до σξ = 2 м, количество реализаций для каждого значения σξM = 500. В таблице представлены: величины смещения оси цилиндра Δтеор, рассчитанные с помощью (2.12); усредненные по 500 реализациям координаты оси цилиндра ${{\hat {X}}_{c}}$, ${{\hat {Y}}_{c}}$, полученные при решении системы уравнений (2.7); усредненная оценка радиуса цилиндра $\hat {R}$, найденная с помощью решения уравнения (2.13); усредненная экспериментальная оценка величины смещения Δэксп, полученная путем подстановки в (2.12) оценки радиуса цилиндра; скорректированное положение оси цилиндра ${{\bar {Y}}_{c}}$, т.е. ${{\bar {Y}}_{c}} = {{\hat {Y}}_{c}} + {{\Delta }_{{{\text{эксп}}}}}$.

Из таблицы видно, что даже при значительной величине СКО измерения дальности оценка радиуса и положения оси цилиндра производятся с большой точностью, следовательно, полученные выражения корректны. Было проведено моделирование предложенной в разд. 2 итерационной процедуры определения параметров цилиндрической поверхности по выборке измерений. Моделирование подтвердило, что процедура быстро сходится за 2–4 итерации (в зависимости от величины ошибок измерения дальности Dξ). В качестве примера на рис. 4 представлены результаты моделирования для σ(ξD) = 4 м. Штрихпунктирной линией обозначена проекция цилиндрической поверхности с радиусом R = 15 м на горизонтальную плоскость. Точки – проекции измерений на горизонтальную плоскость. Координаты начального положения оси цилиндра (отмечен крестиком) Xн = 0.13 м, Yн = 1994.63 м. Начальный радиус Rн = 12.21 м. Дуга показана сплошной линией. Начальное смещение положения оси составило Δ = 5.37 м.

Рис. 4.

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

На рис. 5 приведены результаты реализации четырех итераций (наложены друг на друга). Получены следующие значения: оценка радиуса цилиндра R4 = 15.03 м, ошибка определения положения оси цилиндра по продольной координате ΔY = 0.07 м.

Рис. 5.

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

Моделирование определения положения центра и радиуса сферической поверхности также подтвердило корректность выражений (3.2), (3.4), (3.5), (3.8) и (3.9).

5. Алгоритм выделения цилиндрических поверхностей при сегментации реальных ЛЛ-изображений. Полученные выражения для определения координат оси цилиндра и его радиуса были использованы при разработке алгоритма выделения цилиндрических поверхностей при сегментации реальных дальностных ЛЛ-изображений. Алгоритм заключается в следующем. ЛЛ-изображение обрабатывается “скользящим окном” m × n, горизонтальные размеры которого m не превосходят ожидаемую величину диаметра цилиндров (в пикселах), находящихся на сцене. При движении “скользящего окна” существуют положения, когда “окно” расположено на образе цилиндра симметрично относительно его оси. Поэтому для каждого положения “скользящего окна” выполняется проверка нескольких статистических гипотез. Каждая гипотеза заключается в том, что “окно” полностью находится на цилиндрической поверхности так, что ось цилиндра соответствует середине “окна” и горизонтальный размер “окна” составляет некоторую заранее заданную часть диаметра (например, k/10 при k = 1, …, 10). Альтернативой данным гипотезам рассматривается гипотеза, что “окно” расположено на поверхности общего вида или не симметрично относительно оси цилиндра. Для проведения такой проверки гипотез определяются координаты оси цилиндра (xc, yc) с помощью выражений (2.5), (2.6). Необходимо отметить, что если анализируемая поверхность не цилиндрическая, а является поверхностью произвольного вида, будут найдены координаты некоторой точки, не имеющей физического смысла. Затем для каждой гипотезы о соотношении горизонтального размера “окна” и диаметра цилиндра рассчитываются радиус Rk и смещение оси Δk от ее истинного положения с помощью выражений (2.11)–(2.13) и (2.18), (2.19), после чего координаты оси цилиндра корректируются на величину смещения Δk. Если какая-либо гипотеза верна, то положение оси и радиус цилиндра найдены корректно и, следовательно, математические ожидания измерений лежат на поверхности цилиндра, математическое описание которой получено. В соответствии с (1.1) измерения имеют относительно этой поверхности нормальное распределение с нулевым математическим ожиданием. Тогда нормированные ошибки измерений (в результате деления на СКО σξ) являются независимыми случайными величинами с единичной дисперсией и нулевым математическим ожиданием. Поэтому сумма по всем пикселам окна квадратов таких величин, т.е. Sk:

${{S}_{k}} = \sum\limits_{i = 1}^n {\frac{{{{{([{{y}_{c}} + \Delta - \sqrt {{{R}^{2}} - {{{({{x}_{i}} - {{x}_{c}})}}^{2}}} ] - {{y}_{i}})}}^{2}}}}{{\sigma _{\zeta }^{2}}}} ,$
подчинена закону распределения χ-квадрат и доверительная вероятность события ${{S}_{k}} < {{S}_{{{\text{пор}}}}}$ при известном числе измерений дальностей в “окне” n однозначно связана с пороговым значением ${{S}_{{{\text{пор}}}}}$, которое может быть найдено по таблицам для распределения χ-квадрат с n степенями свободы. Если же гипотеза не верна, то величина Sk имеет распределение, отличное от распределения χ-квадрат, причем доверительная вероятность $P({{S}_{k}} < {{S}_{{{\text{пор}}}}})$ будет существенно меньше, чем в случае цилиндрической поверхности. Для принятия или непринятия перечисленных гипотез можно использовать критерий согласия, основанный на использовании табулированных значений распределения χ-квадрат. В отличие от классического применения критериев согласия, например критерия Пирсона, нет необходимости разбивать выборку измерений на интервалы, определять частоты попаданий измерений в эти интервалы и сопоставлять их с теоретическим распределением. Величина Sk по построению имеет распределение χ-квадрат при верности k-й гипотезы и не имеет распределение χ-квадрат в противном случае. Задаваясь желаемой доверительной вероятностью α и сопоставляя величину Sk с пороговым значением ${{S}_{{{\text{пор}}}}}$ (т.е. с рассчитанной “процентной точкой” распределения χ-квадрат), принимается решение о верности одной из гипотез нахождения окна на цилиндрической поверхности с известным радиусом: при ${{S}_{k}} \leqslant \chi _{\alpha }^{2}$ k-я гипотеза принимается, а остальные отклоняются и всем пикселам фрагмента изображения размером 2Rk × n присваивается определенный индекс. Если указанный критерий выполняется для нескольких гипотез, то выбирается гипотеза, для которой Sk минимально. Если не для одной из гипотез не верно условие ${{S}_{k}} \leqslant \chi _{\alpha }^{2}$, то принимается гипотеза нахождения “окна” на поверхности общего вида. Вычисляется новое положение “скользящего окна”, и цикл повторяется. Предложенный подход позволяет провести выделение цилиндрических поверхностей при сегментации ЛЛ-изображений с заранее заданной доверительной вероятностью и определить положение оси и радиус каждого цилиндра. Последующее построение связных образов цилиндров на ЛЛ-изображении может быть произведено на базе эвристических подходов.

Для выделения сферических поверхностей был разработан аналогичный алгоритм.

6. Результаты моделирования с использованием реальных ЛЛ-изображений. На рис. 6–8 представлен пример исходного изображения и результаты обработки реального ЛЛ-изображения с применением алгоритма обнаружения фрагментов образов цилиндрических поверхностей на сцене по дальностному ЛЛ-изображению. На рис. 6 приведено изображение промышленно-городской сцены в видимом диапазоне, на рис. 7 – бинарное ЛЛ-изображение той же сцены: белый цвет соответствует принятому отраженному сигналу, черный – отраженный сигнал отсутствует. Изображения, показанные на рис. 6–8, сформированы с помощью активно-пассивного многоспектрального комплекса дистанционного зондирования [14]; на рис. 8 представлены результаты обработки ЛЛ-изображения предложенным алгоритмом.

Рис. 6.

Изображение промышленно-городской сцены в видимом диапазоне

Рис. 7.

Бинарное ЛЛ-изображение промышленно-городской сцены: белый цвет соответствует принятому отраженному сигналу, черный – отраженный сигнал отсутствует

Рис. 8.

Результаты обработки ЛЛ-изображения промышленно-городской сцены предложенным алгоритмом

Под каждой сегментированной областью (более темный тон серого) указано число пикселов N в ней, положение оси Xс, Yс и радиус цилиндра R. Отметим, что на рис. 8 выделены не только трубы с априори известной формой поверхности, но и цилиндрический бак для хранения горюче-смазочных материалов (в правом нижнем углу обработанного изображения), конфигурация которого не определяется по ТВ-изображению. Цилиндрическая форма данного объекта определена в результате анализа космоснимка. Таким образом, моделирование с использованием реальных ЛЛ-изображений подтвердило работоспособность предложенного алгоритма.

Заключение. Получено аналитическое решение задачи оценки местоположения фрагмента цилиндрической поверхности и его радиуса при обработке дальностного ЛЛ-изображения. Оценки являются оптимальными для случая, когда СКО центра отражения от оси диаграммы направленности существенно меньше СКО ошибок измерения дальности до отражающей поверхности. Предложенные методы и алгоритмы позволяют осуществлять поиск фрагментов таких поверхностей по ЛЛ-изображениям с заданной доверительной вероятностью. Метод обобщен для решения задачи выделения на ЛЛ-изображениях сферических поверхностей.

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

Работоспособность методов и алгоритмов проверена путем моделирования с использованием реальных ЛЛ-изображений промышленно-городских сцен.

Таблица
Параметры σξ, м
0.2 0.4 0.8 1.2 1.6 2.0
Смещение оси, полученное теоретически, Δтеор, м 0.193 0.620 1.386 1.797 2.005 2.119
Усредненные по 500 реализациям координаты оси цилиндра, м ${{\hat {X}}_{c}}$ 0.0003 0.0002 –0.00014 0.00007 0.0007 0.0053
${{\hat {Y}}_{c}}$ 1999.81 1999.388 1998.622 1998.212 1997.999 1997.885
Усредненная оценка радиуса, $\hat {R}$, м 3.0023 3.0052 3.0059 3.006 3.003 3.005
Экспериментальная усредненная оценка смещения, Δэксп, м 0.193 0.619 1.386 1.798 2.0059 2.121
Скорректированное положение оси, ${{\bar {Y}}_{c}}$, м 2000.003 2000.007 2000.008 2000.011 2000.005 2000.006

Список литературы

  1. Основы импульсной лазерной локации / Под ред. В.Н. Рождествина. М.: Изд-во МГТУ, 2006. 512 с.

  2. Дановский В.Н., Ким В.Я., Лисицын В.М., Обросов К.В., Тихонова С.В. Сравнение возможностей радиолокации и лазерной локации как методов информационного обеспечения безопасности маловысотного полета // Изв. РАН. ТиСУ. 2007. № 4. С. 153–165.

  3. Кузьмин С.З. Основы проектирования систем цифровой обработки радиолокационной информации. М.: Радио и связь, 1986. 352 с.

  4. Pal N.R., Pal S.K.A. Review on Image Segmentation Techniques // Pattern Recognition. 1993. V. 26. № 9.

  5. Gonzalez R.C., Woods R.E. Digital Image Processing. 3rd Edition / Pearson Education, Inc. Prentice-Hall. 2008. 976 p.

  6. Колдаев В.Д. Эвристические и квазитопологические алгоритмы контурной сегментации изображений // Изв. вузов. Электроника. 2008. № 6. С. 41–45.

  7. Фурман Я.А. Сегментация и описание трехмерных структур на базе кватернионных моделей // Наукоемкие технологии. 2007. № 9. С. 37–49.

  8. Zhang J., Lin X., Ning X. SVM-based Classification of Segmented Airborne LiDAR Point Clouds in Urban Areas // Remote Sensing. 2013. V. 5 (8). P. 3749–3775.

  9. Lin X., Zhang J. Segmentation-Based Filtering of Airborne LiDAR Point Clouds by Progressive Densification of Terrain Segments // Remote Sensing. 2014. V. 6 (2). P. 1294–1326.

  10. Лисицын В.М., Обросов К.В., Пасечный Н.Н., Стефанов В.А. Алгоритм сегментации доплеровских оптических изображений // Изв. АН СССР. Техн. кибернетика. 1990. № 2. С. 203–213.

  11. Lisitsyn V., Obrosov K., Pasechny N., Stefanov V. Iterative Algorithm for Multidimensional Image Analysis // International Archives of ISPRS. 1992. V. 29. Pt 3-1. P. 527–534.

  12. Лисицын В.М., Обросов К.В., Пасечный Н.Н., Стефанов В.А. Итерационные алгоритмы сегментации многомерных изображений // Изв. РАН. Техн. кибернетика. 1993. № 6. С. 103–113.

  13. Lisitsyn V., Pasechny N., Stefanov V., Samoilov V., Lukanidin D. Application of Random Decision Rules in Laser Locator Image Optimal Segmentation Algorithm // International Archives of ISPRS. 1996. V. 31. Pt B-3. P. 463–466.

  14. Лисицын В.М., Ким В.Я. Активно-пассивный многоспектральный комплекс дистанционного зондирования // Лазеры в науке, технике, медицине: сб. научн. тр. Т. 23 / Под ред. В.А. Петрова. М.: МНТОРЭС им. А.С. Попова, 2012. С. 16–21.

Дополнительные материалы отсутствуют.