Способ идентификации загрязнений морской поверхности

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

 

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

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

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

Известен «Способ идентификации типов растительности», патент RU №2242716, A01G, 23/00, 2004 г., - аналог. В способе-аналоге получают изображение подстилающей поверхности в виде зависимости спектральной яркости I (х, у) от пространственных координат, разбивают изображение на мозаику участков, вычисляют фрактальную размерность и автокорреляционную функцию сигнала изображения каждого участка и сравнивают их с эталонными, изображение регистрируют в зеленой полосе видимого спектра по двум взаимно ортогональным по поляризации каналам приема, рассчитывают средний уровень сигнала изображения в каждом канале, вычисляют попиксельные отношения I₁ (x, у), I₂ (x, у) изображений с большим средним к меньшему и формируют синтезированную матрицу изображений из этих соотношений, методами пространственного дифференцирования выделяют контуры на синтезированном изображении, рассчитывают числовые характеристики сигнала фрагментов изображений внутри выделенных контуров и по значениям коэффициентов фрактальной размерности и ширине автокорреляционной функции сигнала судят о принадлежности фрагмента изображения к данному типу растительности на нем.

К недостаткам способа-аналога следует отнести:

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

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

Ближайшим аналогом по технической сущности к заявленному является «Способ обнаружения аномалий морской поверхности», патент RU №2109304, G01S, 11/06, 13/89, 1997 г.

Способ ближайшего аналога включает получение изображения морской поверхности в виде матрицы цифровых отсчетов |m×n| элементов функции яркости I (х, у) от пространственных координат, обработку матрицы путем разбиения на мозаику фрактальных участков, вычисление огибающей пространственного спектра и автокорреляционной функции сигнала каждого участка, расчет интегрального признака z=R/В и сравнение его с фоновым z₀=R₀/В₀, вывод на отображение участков, для которых z/z₀>2, синтезирование из последовательности проанализированных участков мозаичной картины аномалии, где В, В₀ - максимальные значения автокорреляционных функций электрического сигнала матриц аномалии и фона соответственно, R, R₀ - ширина автокорреляционных функций на уровне 0,1 их максимального значения для аномалии и фона.

Недостатками способа ближайшего аналога являются:

- недостаточная достоверность результата из-за неиспользования при обработке сигнала всех координат взволнованной морской поверхности: х, у, z;

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

- не все из возможных параметров сигнала учитываются при расчете интегрального критерия.

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

Поставленная задача решается тем, что способ идентификации загрязнений морской поверхности включает зондирование прибрежных акваторий, содержащих эталонные участки, средствами, установленными на воздушно-космическом носителе, с получением изображений в ультрафиолетовом и красном участках солнечного спектра, привязку изображений по координатам системой позиционирования «ГЛОНАСС», формирование синтезированной матрицы из попиксельных отношений этих изображений, выделение контуров областей загрязнения программным расчетом градиента функции яркости I(х, у) синтезированного изображения, вычисление следующих параметров внутри выделенных контуров: среднего значения F частоты пространственного спектра функции яркости I(х, у) изображения, фрактальной размерности Ω изображения, площади S_p рельефа для анализируемого и соответствующего ему эталонного участка, определение параметра идентификации для анализируемого и соответствующего ему эталонного участка как , где S₀ геометрическая площадь участка, определение разности ΔП для анализируемого и эталонного участков, оценку уровня загрязнения в процентах через отношение ΔП к П эталонного участка.

Изобретение поясняется чертежами, где:

фиг. 1 - зависимость коэффициента отражения от длины волны падающего светового потока;

фиг. 2 - выделенный контур загрязнения поверхности;

фиг. 3 - амплитудно-частотные характеристики эталонной а) и загрязненной б) поверхностей;

фиг. 4 - функции фрактальной размерности эталонной в) и загрязненной г) поверхностей;

фиг. 5 - относительные S_p/S₀ площади рельефов эталонной д) и загрязненной е) поверхностей в зависимости от балльности волнения;

фиг. 6 - функциональная схема устройства, реализующая способ.

Техническая сущность изобретения состоит в следующем. При загрязнении воды взвесями, пленками органических веществ, планктоном изменяется ее коэффициент поверхностного натяжения в пределах от (0,071…0,03) Н/м. Изменение коэффициента поверхностного натяжения приводит к демпфированию мелкодисперсной ветровой ряби морского волнения на участках загрязнения и, как следствие, к изменению спектра пространственных волн. Одновременно с изменением спектра изменяется и коэффициент преломления морской воды, что приводит к изменению формы волнения и коэффициента отражения падающего светового потока.

Основным методом идентификации объектов является сравнение с эталоном. В качестве идентифицируемых параметров в заявленном способе измеряются:

- пространственный спектр волнения, F;

- форма волнения, технической характеристикой которого является фрактальная размерность, Ω;

- относительная площадь рельефа пикселей изображения, S_p/S₀.

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

Коэффициент отражения К электромагнитного поля от подстилающей поверхности определяется соотношениями Френеля:

где - комплексная диэлектрическая проницаемость отражающей поверхности;

γ - угол скольжения, для зондирования в надир γ→π/2.

Для зондирования в надир коэффициент отражения в первом приближении будет равен:

где n - коэффициент преломления.

Коэффициент преломления n существенно зависит от длины волны λ светового потока. Эта зависимость аппроксимируется числовым рядом Коши вида:

где А, В, С - коэффициенты аппроксимации.

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

Психологически восприятие образа объекта человеком-оператором происходит на уровне контуров. Последнее достигается путем выделения контуров (контурного рисунка) на изображениях, осуществляемого методами пространственного дифференцирования [см., например, Дуда P.O., Харт П.Е. «Распознавание образов и анализ сцен», перевод с англ., изд. Мир, М., 1976 г. «Пространственное дифференцирование», стр. 287-288]. Существует несколько стандартных операторов (Робертса, Лапласа, Собела), позволяющих вычислить контуры на двумерных изображениях. В частности, например, перекрестный оператор Робертса рассчитывают в окне (маске) 2×2 элемента, для каждой дискретной точки изображения I (х, у) из отношения:

Вычисленное значение оператора сравнивают с установленным порогом R (i, j)>П, выводят на экран точки, для которых R (i, j) превышает порог. Выделение контуров на изображении с использованием масок различных операторов представляется стандартной математической операцией [см., например, П.А. Минько. «Обработка графики Photoshop CS2», из-во Эксмо, 2007 г., стр. 47-56]. После выделения контура на синтезированном изображении рассчитывают количественные значения идентифицирующих параметров.

Падающий световой поток на взволнованную морскую поверхность по-разному отражается от него. Гребни волн отражают падающий поток практически зеркально, в то время как скаты волн - диффузно. Поэтому текстура изображения морской поверхности повторяет ее геометрию, т.е. изображение содержит информацию о спектре волнения. В соответствии с теоремой Котельникова-Шеннона, непрерывная функция однозначно определяется своими дискретными отсчетами через интервал где F_max - максимальная частота спектра функции [см., например, Теоретические основы радиолокации, под ред. В.Е. Дулевича, Сов. Радио, М., 1964 г., стр. 212]. В соответствии со шкалой Бофорта наименьшая длина волны (рябь) составляет 0,3 м. Следовательно, для неискаженного восстановления спектра ветрового волнения по его изображению видимого диапазона пространственное разрешение цифровой видеокамеры должно составлять порядка 15 см на пиксель. Существующие цифровые видеокамеры обеспечивают требуемое пространственное разрешение.

Амплитудно-частотный спектр G (F_x, F_y) волнения может быть восстановлен по его изображению расчетом Фурье-преобразования матрицы отсчетов из |m×n| элементов в соответствии с зависимостью:

где F_x, F_y - спектр волнения по координатам х, у;

I (х, у) - функция яркости изображения;

m, n - число строк, столбцов матрицы |m×n|.

Расчет спектра осуществляют алгоритмами быстрого Фурье-преобразования (БФП) по стандартным программам, входящим в комплект специализированного программного обеспечения ПЭВМ типа MATH CAD, ER MAPPER [см., например, НТИ «Специализированное программное обеспечение MATH CAD 6.0 PLUS», издание 2-е стереотипное, М.: Информационно-издательский дом «Филинъ», 1997 г., стр. 441]. АЧХ пространственных спектров двух изображений а) эталонного участка и б) загрязненного участка иллюстрируются графиками фиг. 3.

Одновременно установлено, что наибольший объем информации об объекте содержит его форма. Элементом формы объекта по Мандельброту является его фрактал [см., например, Mandelbrot В. Fractals, Forms, Chance and Dimensions, Freeman, San Francisco, 1977].

Фрактальная размерность является числовым параметром, характеризующим структуру природных образований, в частности, для изображения этот параметр заключается в интервале [2…3]. Для вычисления фрактальной размерности используют метод осцилляций.

Пусть (x₁, y₁) и (х₂, у₂) - двумерные координаты точек, а третья координата, яркость, задана в виде функции координат I (х, у).

Тогда ε - осцилляцией значений (I) будет разность наибольшего и наименьшего значения (I) в (ε) в окрестности (х, у).

После этого ε - вариацию значений I вычисляют как:

где a, b - пределы, в которых изменяется переменная х;

с, d - пределы, в которых изменяется переменная у.

Фрактальная размерность матрицы вычисляется как размерность Хаусдорфа:

Вычисление фрактальной размерности изображений объектов осуществляют по специализированной программе. Текст программы приведен ниже в примере реализации. Расчетные значения функций фрактальной размерности эталонного в) и загрязненного г) участков иллюстрируются графиками фиг. 4.

Кроме длины волны в координатах х, у ветровое волнение по шкале Бофорта характеризуется высотой волн по координате Z. Именно по этой координате в областях аномалий наблюдаются наибольшие изменения. В пространстве данная координата определяет площадь рельефа взволнованной поверхности. Площадь рельефа Sp вычисляют как поверхностный интеграл функции z (х, у) [см., например, Н.С. Пискунов. Дифференциальное и интегральное исчисления для ВТУЗов, учебник, 5-е издание, М.: Наука, 1964 г., § 7. Вычисление площади поверхности, стр. 73-74]:

Поскольку аналитическую зависимость z (х, у) получить невозможно, для вычисления S_p используют численные методы программного расчета [см., например, «Способ определения площади рельефа», Патенты RU №2251075, №2255357, 2005 г.]. Текст программы расчета площади раельефа пикселей изображения приведен ниже в примере реализации. Результат расчета относительной площади рельефа S_p/S₀ иллюстрируются графиками фиг. 5.

Пример реализации способа

Заявленный способ может быть реализован по схеме фиг. 6. Функциональная схема содержит авиационный носитель (1) (типа самолетной лаборатории, созданной в рамках международной программы «Открытое небо») с установленными на нем цифровой видеокамерой (2) инфракрасного диапазона (типа «Фиалка МВ-Космос») и спектрозональной камеры (3) (типа «MOMS-2P», Германия). Трассовую покадровую съемку запланированных акваторий в полосе сканирования (4) осуществляют от бортового комплекса управления (БКУ) (5) на основе заложенных в БКУ программ включения бортовой аппаратуры. Результаты покадровой съемки участков акваторий записывают в бортовое запоминающее устройство (6) с одновременной привязкой по координатам изображений от аппаратуры потребителей (7) космической системы позиционирования «ГЛОНАСС». После посадки самолетного носителя массивы полученных измерений помещают на сервер хранения данных (8). Тематическую обработку изображений осуществляют в центре обработки (9), где через устройство ввода и передачи (10) информация из сервера хранения поступает в электронно-вычислительную машину (11) со стандартным набором периферийных устройств: процессор (12), оперативное запоминающее устройство (13), винчестер (14), дисплей (15), принтер (16), клавиатура (17), сервер Internet (18). Предварительно в оперативное запоминающее устройство (13) записывают программу специализированного программного обеспечения MATH CAD. Затем формируют кадры синтезированных матриц из попиксельных отношений ультрафиолетового изображения и зонального красного изображения.

Данная операция реализуется использованием специализированного программного обеспечения [см., например, MATH CAD. 7.0. PLVS, издание 3-е стереотипное, Информационно-издательский дом «Филинъ», 1998 г., стр. 211, Векторизация элементов матрицы]. После этого осуществляют нормирование функции сигнала синтезированной матрицы в стандартной шкале 0…255 уровней квантования.

За счет контрастирования изображений достигается увеличение чувствительности тракта обработки, что обеспечивает достоверное выделение контуров на синтезированном изображении программным методом [см., например, П.А. Минько. «Обработка графики Photoshop CS2», изд-во «Эксмо», 2007 г., стр. 47-56, Глава 3. Методы выделения областей].

Результат выделения области загрязнения иллюстрируется фиг. 2.

Затем алгоритмами быстрого Фурье-преобразования по программам специализированного программного обеспечения MATH CAD вычисляют пространственный спектр функции яркости I (х, у) внутри выделенного контура. За среднее значение частоты пространственного спектра принимают значение, делящее площадь под графиками фиг. 3 пополам. Эти значения соответственно составили: для загрязненного участка F₁=0,3; для эталонного - F₂=0,5.

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

Текс программы вычисления фрактальной размерности изображений

Фрактальная размерность изображения, по Хаусдорфу, занимает интервал [2…3].

Результат программного расчета иллюстрируется графиками фиг. 4. Фрактальная размерность эталонного участка Ω_в=2,4; загрязненного - Ω_г=2,7.

Для численного расчета площади рельефа пикселей изображения последовательно, от начала матрицы, изображение разбивают на окна |2×2| элемента. Аппроксимируют площадь рельефа в окне мозаикой треугольников. При известном пространственном разрешении одного пикселя (Δх, Δy) площади треугольников в окне находят по формуле Герона:

где а, b, с - длины сторон треугольника; р - полупериметр. По теореме Пифагора рассчитывают длины сторон а, b, с [см. патент RU №2255357]. Поверхность рельефа S_p находят как сумму площадей треугольников каждого окна. Геометрическую площадь S₀ участка определяют как сумму площадей пикселей обрабатываемого участка.

Расчет осуществляют программным методом на ПЭВМ.

Текст программы

Результат программного расчета площади S_p/S₀ иллюстрируются графиками фиг. 5.

Для балльности ветрового волнения по шкале Бофорта, равной 3, эти значения составили: эталонного участка S_p/S₀=2; загрязненного - S_p/S₀=1,6. В целом для расчетных графиков фиг. 3, 4, 5 показатель идентификации эталонного участка: текущего участка разница ΔП=0,237; показатель загрязненности:

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

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

Способ идентификации загрязнений морской поверхности, включающий зондирование прибрежных акваторий, содержащих эталонные участки, средствами, установленными на воздушно-космическом носителе, с получением изображений в ультрафиолетовом и красном участках солнечного спектра, привязку изображений по координатам системой позиционирования ГЛОНАСС, формирование синтезированной матрицы из попиксельных отношений этих изображений, выделение контуров областей загрязнения программным расчетом градиента функции яркости I(x, y) синтезированного изображения, вычисление следующих параметров внутри выделенных контуров: среднего значения F частоты пространственного спектра функции яркости I(x, y) изображения, фрактальной размерности Ω изображения, площади S_p рельефа для анализируемого и соответствующего ему эталонного участка, определение параметра идентификации для анализируемого и соответствующего ему эталонного участка как , где S₀ - геометрическая площадь участка, определение разности ΔП для анализируемого и эталонного участков, оценку уровня загрязнения в процентах через отношение ΔП к П эталонного участка.