Содержание к диссертации
Введение
1 . Моделирование сигналов лазерной триангуляционной системы
1.1. Вводные замечания 15
1.2. Построение модели сигнала с фотоприёмника 16
1.3. Анализ эффективности моделирования сигнала 22
1.4. Выводы 30
2. Синтез и анализ алгоритма оценивания положения одиночного видеоимпульса
2.1. Вводные замечания 32
2.2. Обзор методов получения эталонных оценок 35
2.3. Синтез косвенного метода оценки положения импульса 41
2.4. Анализ косвенного метода оценки положения импульса 50
2.5. Параметрическая оптимизация косвенного алгоритма 59
2.6. Выводы 63
3. Адаптивная калибровка лазерных триангуляционных систем
3.1. Вводные замечания 65
3.2. Двухпараметрическая калибровка триангуляционного сенсора 66
3.3. Интерполяция двумерной калибровочной зависимости на неравномерной сетке 73
3.4. Экстраполяция калибровочной зависимости 81
3.5. Выводы 81
4. Вычислительная оптимизация алгоритмов первичной обработки в лазерных триангуляционных системах
4.1. Вводные замечания 83
4.2. Комбинированный метод нахождения центра тяжести 83
4.3. Рекурсивное преобразование косвенного алгоритма нахождения положения импульса 87
4.4. Оптимизация логического блока косвенного метода 93
4.5. Параметрическая оптимизация цифровых фильтров с квантованными коэффициентами 96
4.6. Выводы 104
5. Применение алгоритмов первичной обработки в лазерных триангуляционных системах для задач технической диагностики .
5.1. Вводные замечания 106
5.2. Повышение точности оценки параметров подвижных объектов 109
5.3. Оптический виброметр с высокой разрешающей способностью 119
5.4. Визуальный конструктор триангуляционного сенсора 124
5.5. Выводы 128
Заключение 130
Список использованных источников 133
Приложения 147
- Построение модели сигнала с фотоприёмника
- Синтез косвенного метода оценки положения импульса
- Интерполяция двумерной калибровочной зависимости на неравномерной сетке
- Рекурсивное преобразование косвенного алгоритма нахождения положения импульса
Построение модели сигнала с фотоприёмника
Моделирование сигнала с фотоприёмника триангуляционного сенсора начинается с разделения сигнала на компоненты - различные составляющие шумов и полезный импульс. Внешняя засветка обладает высокой межкадровой корреляцией и может быть исключена в результате её маскирования предыдущим кадром с фотоприёмника с привлечением техники блинкования [64]. В дальнейшем под сигналом s с триангуляционного сенсора будем подразумевать аддитивную смесь шумов и полезного импульса х, образованного распределением интенсивности в сечении луча лазера.
Построение модели импульса и сигнала в целом происходит в статическом режиме, т.е. неограниченное время и произвольным образом. Поэтому выделение границ импульса в составе сигнала не представляет труда и может быть произведено в том числе и визуально по реализации. В результате выбирается диапазон индексов jx...(jx + А/) сигнала s в котором находится импульс х = {0,... ,x[jx],...,x[jx+Aj],... ,0}, где А/ - ширина основания импульса (вне указанного диапазона считается, что х[/] = 0). При подборе статистической модели шумов целесообразно применять теоретическое обоснование применения рассматриваемых распределений, т.е. физико-статистический подход к выбору распределения [65]. Это предполагает выбор предельного распределения на основе физических и технологических предпосылок его формирования. Такой подход для конкретной задачи позволяет создать наиболее адекватную и простую модель. Рассмотрим источники шумов. При моделировании шума следует разделять фоновые шумы f= {/[/]} и шумы р ={/?[/]} накладываемые на импульс х= [x\j]}.
В состав фоновых шумов, распределенных на протяжении всей длительности N сигнала s= {s\j]}, неизменно входят некоррелированные шумы оптического и электрического происхождения. Появление этих компонент связано с факторами внешней среды и трактом прохождения сигнала от фотоприёмника до устройства обработки. Также обязательно присутствует равномерно распределенный шум квантования. Отметим, что амплитуда шума квантования не превышает 0,5 кванта, что обычно на порядок меньше среднеквадратического значения прочих шумовых компонент. Поэтому основной вклад в фоновый шум вносит некоррелированная составляющая, обладающая большей интенсивностью. В этой связи в большинстве практических случаев качество оценки местоположения центра импульса незначительно зависит от полноты описания фоновых шумов, т.к. часто имеется одна преобладающая по интенсивности компонента. Можно предположить, что фоновый шум, как и большинство случайных процессов в природе, будет иметь нормальное распределение Л„, вида [66]: ЛД/) = 0,5Р„. ==ехр{- 1, (1.1) где / = 0...2Г-1, г - разрядность представления данных, а„ -среднеквадратическое значение шума, тп - среднее значение шума, Рп -площадь шума (индекс п показывает отношение распределения и его параметров к модели шума п = {«[/]}). Площадь шума сводится к выражению: N J l, где n\J] - отсчёты шума, N - длительность сигнала s, J - длительность анализируемого участка сигнала s вне диапазона импульса/ .. .(jx + А/). Средний уровень шума f устраняется при маскировании, т.е. ш/= 0, а дисперсия a2f находится численно как второй момент [66]. На рис. 1.1 приведена эпюра сигнала s с фотоприёмника реального триангуляционного измерителя после исключения коррелированной помехи. На рис. 1.2 «квадратами» показано распределение А/ фрагмента сигнала s длиной N= 1024, соответствующего фоновому шуму f, а сплошной линией -нормальное аналитическое распределение вида (1.1) с параметрами Pf, т/= 0 и j2f (индекс/показывает отношение данных параметров к фоновому шуму f).
Как легко увидеть из рис. 1.2, распределение А/ фонового шума с достаточной точностью может быть описано нормальной плотностью вероятности. Поскольку природа появления фонового шума для всех триангуляционных сенсоров одинакова и не зависит от их конструктивных особенностей, то можно утверждать, что компонента f всегда будет иметь распределение близкое к нормальному.
Синтез косвенного метода оценки положения импульса
Дисперсия оценки положения одиночного видеоимпульса преимущественно вызвана наличием аддитивного шума вследствие ошибки измерения, регистрации и влияния прочих факторов формирования сигнала. Поэтому требуется снизить влияние шума, т.е. произвести процедуру фильтрации низкочастотного по сравнению с шумом сигнала, которая записывается в виде линейной свёртки сигнала и импульсной характеристики (ИХ) фильтра: JV-I 4j] = Yg[k]-s[j + kl =0 где у = -iV.. .N; g = {g[k]} - вектор отсчётов ИХ; s = {s[j + к]} - входной сигнал представленный N значащими отсчётами в пределах 0 (/ + к) N, т.е. при (J + k) N или (j + k) 0 полагается, что s[j + k]=0, Чтобы в результате фильтрации не получить нерегулярного смещения, вызванного фазовыми искажениями сигнала при прохождении фильтра, фазочастотная характеристика (ФЧХ) фильтра должна быть линейной. Такая ФЧХ характерна для фильтров с нечётной ИХ. В матричном виде линейная свёртка (2.15) записывается в виде w = Sg, где w= {w[/]}, S - Как видно из выражения (2.16), сигнал на выходе фильтра имеет 2N отсчётов, т.е. N отсчётов при потактовом вхождении значащих отсчётов сигнала s в линию задержки фильтра и TV отсчётов при выходе сигнала s из фильтра. На рис. 2.4 штриховкой показана область индексов (/ + к) в которой осуществляется вычисление свёртки (2.15); сплошной линией показан входной сигнал s в различных областях: на входе линии задержки фильтра, полностью в её пределах и на её выходе; пунктиром показана произвольная нечётная ИХ g (опустим пока рассуждения о форме ИХ, ограничившись тем условием, что она является нечётной функцией к). Сигнал на выходе ЛЗ фильтра
Как видно из анализа структуры матриц Si и S2 их элементы дополняют друг друга, т.е. отсутствующим (нулевым) элементам матрицы Si соответствуют ненулевые элементы матрицы S2 и наоборот. Поэтому без потери информации можно перейти к более компактной записи в виде циклической свёртки: v = Cg = (S, + S2)g. (2.17) или в традиционном виде для циклической свёртки, которая находится как поэлементная сумма вида (w\j] + w\j + Щ) двух частей одной линейной свёртки (2.15) длиной TV: v[j] = Y,g[k]-s[(j + k)modN], (2.18) =0 гдеу = О...ЛЧ,у={у[/]}. В выражении (2.18) функция g[k] имеет смысл дискриминационной характеристики и, в частности, свёртка (2.18) для у = 0 характеризует отклонение сигнала s от среднего положения N12. Представим выражение (2.18) как два слагаемых и заменим пределы суммирования с учётом нечётности g: NI2-\ N-\ v[j]= g[k]-s[(j + k)modN]+ %g[k]-s[ j + k)modN] = i=0 k=NI2 = Nfg[k]-s[(j + k)modN]-Yg[k]-s[(N/2-U + k))modN]= (2.19) = Ng[k]is[U + k)modN]-s[{N/2-(j + k))modN]). Выражение (2.19) показывает, что принципиального значения форма g не имеет, и влияет она только на вид функциональной зависимости рассогласования площадей S+ и S_ сигнала до и после индекса N12 от положения центра сигнала (см. рис. 2.5): M(M) S+-S_, W/2-1 W/2-1 где S+ = YJSU] S- = YJS№ /2 + j], М- центр сигнала s (мода импульса x). ;=o ;=o Справедливость проведённых рассуждений подтверждается частным случаем - совпадением значения свёртки (2.18) при/ = 0 для случая линейной g с центром тяжести: M = fJg-s[g] (2.20) при соблюдении условия нормировки: 5M ] = 1- (2-21) Для иллюстрации частного случая введем в выражение (2.18) смещение ИХ в виде g [k]=g[k] + с: М = [к] [к] = к-5[к]- -8[к] + с [к], (2.22) =0 к=0 к=0 2 к=0 где g[k] = к-N12 - нечётная линейная ИХ, с = N12 - смещение линейной ИХ для приведения её в область значений g\k] є [0.. JV-1], т.е. g [k] = к. Выражение (2.22), в силу условия нормировки (2.21) и с учётом с = N12, сводится к равенству M = Y k-s[k], (2.23) что в отличие от (2.18) даёт не отклонение центра сигнала от среднего положения, а непосредственно оценку местоположения центра сигнала, являясь полностью эквивалентной оценке (2.20) центра тяжести. Поэтому введение постоянного смещения ИХ является корректным, т.к. сохраняется условие нечётности ИХ - нечётность наблюдается относительно точки (v[c]; с).
Проведем исследование обобщенной записи (2.18) и рассмотрим выходную последовательность (рис. 2.6) свёртки (2.18) для индексов сдвига j = 0...N-\, где последовательность s представлена в виде сигнала гауссовой формы (у! = у2 = 0, см. гл. 1) с центром в точке М= \i\ = N/2 (рис. 2.7).
Значение v[0] соответствует оценке положения импульса по методу центра тяжести. В то же время, по оси абсцисс можно видеть индекс сдвига j = N/2, которому соответствует значение vnop = v[N/2] = N/2. Отметим большую крутизну участка кривой v = {v[/]} в области j = N/2 по сравнению с участком вблизи 0-го индекса. Очевидно, что если находить оценку положения импульса косвенно - как координату точки пересечения кривой v с её уровнем N12 (средним), то мы получим меньшее рассеяние оценки [78] по сравнению с разбросом значений v[0]. Свёртка гауссова импульса Гауссова модель импульса N12 Рис. 2.6 N12 Рис. 2. v/J N12 N-\J м 0,5 1 ЛЧ,; На рис. 2.8 штриховкой показана область значений свёртки (2.18) при наличии аддитивного гауссова шума фиксированной интенсивности. Также на рис. 2.8 показаны величины Av и ДМ, характеризующие рассеяние оценок по методу центра тяжести и по косвенному методу соответственно. Область значений свёртки при аддитивном белом шуме N12 N12 Рис. 2.8 N-IJ Как видно из рис. 2.8, косвенная оценка имеет заметно меньший разброс значений - AM, по сравнению с разбросом Av. Так, при отношении сигнал/шум, составляющем единицы (3...4), выигрыш Av/AM достигает десятков раз в отсутствие предварительной пороговой обработки (1.5).
Интерполяция двумерной калибровочной зависимости на неравномерной сетке
Выше было отмечено, что даже аналитический расчёт калибровочной зависимости может не обеспечить покрытие с заданной дискретностью всех значений дальности и ширины пятна. Поэтому в таблице будут присутствовать пропуски. Более того, при практических измерениях они могут охватывать значительные области калибровочной зависимости, ввиду ограниченного количества замеров. Очевидно, также, что в результате флюктуации, оценки М и L могут не соответствовать точно узлам таблицы.
В результате полученная двумерная калибровочная зависимость реализована на неравномерной сетке. Непосредственную интерполяцию в этих условиях можно произвести с помощью методов сплайн-интерполяции [89] или линейной интерполяции. Основной недостаток сплайн-интерполяции заключается в её неустойчивости к ошибке вычисления производных в узловых точках, что существенно затрудняет её применение. Прочие методы, например, регрессия [65], могут привести к менее качественным результатам, т.к. не гарантируется сохранение узловых точек. Поскольку калибровочная зависимость характеризуется малой кривизной на всех своих участках (см. рис. 3.3), то именно линейная интерполяция оказывает предпочтительной в силу своей устойчивости к ошибкам оценивания параметров калибровки.
Линейная интерполяция производится внутри треугольников, на которые разбивается калибровочная зависимость. Каждый треугольник представлен тремя вершинами, по которым рассчитываются коэффициенты уравнения плоскости вида H(M,L) = m-M+l-L + bp из системы линейных уравнений: Нх = т М, + / Lx + bp H2=m-M2+l-L2+bp, (3.4) Н3=т-М3+1-Ь3+Ьр где т, I и Ьр - искомые коэффициенты уравнения плоскости, а Н„, Мп и Ln -известные координаты узловых точек. Точность решения системы ограничена только разрядностью представления чисел, что для рассматриваемой задачи несущественно. После получения коэффициентов уравнения плоскости производится расчёт пропущенных точек с желаемой дискретностью. Фактическим ограничением линейной интерполяции (как, впрочем, и любого другого метода интерполяции) является требование получения достаточного числа узловых точек на границах диапазонов по дальности и ширине пятна.
Таким образом, основой линейной интерполяции, обеспечивающей переход к равномерной сетке с заданной дискретностью, является процедура триангуляции. В данном контексте триангуляцией называется разбиение поверхности, представленной набором точек в трехмерном пространстве, на треугольники [90]. Применительно к калибровке измерителей, процедура триангуляции отличается от строгого определения тем, что внешний контур поверхности является конечным - он ограничен диапазонами 0...N изменения оценок М и L.
Вопрос преобразования набора точек в треугольники широко исследуется в компьютерной графике и ему посвящено множество работ. Среди наиболее цитируемых российскими разработчиками работ по машинной графике и триангуляции, в частности, кроме [90] можно назвать [91-93].
Традиционные операции триангуляции приведены в таблице 3.1 [90]. В тех или иных алгоритмах некоторые из этих операций могут не использоваться. Известные методы триангуляции ориентированы на специфические особенности поверхности, представленной набором точек, и оптимизируются по скорости. Последнее для задачи калибровки является несущественным, а особенность калибровочной зависимости позволяет разработать упрощенный алгоритм, ориентированный на обработку трехмерной поверхности, не имеющей экстремумов и «окон», т.е. незакрашиваемых областей [94].
Упрощенный алгоритм может быть представлен меньшим числом операций, которые приведены в таблице 3.2. Операция 1 служит для обозначения начала поиска и выполняется только один раз. Операции 2 и 3 далее повторяются для вновь найденных треугольников.
Рекурсивное преобразование косвенного алгоритма нахождения положения импульса
Основой предлагаемого в п. 2.2 косвенного алгоритма нахождения положения полезного импульса является процедура вычисления циклической свёртки. Известно три основных метода, различающихся точностью и объемом вычислительных затрат [80]. 1) Первый метод основан при применении быстрого преобразования Фурье [101] или быстрого преобразования Хартли [102] и обеспечивает существенное сокращение арифметических операций, длина N которых больше 32. К недостаткам следует отнести высокие ошибки округления и большой объем памяти. 2) Второй метод, использующий теоретико-числовые преобразования, является точным, т.к. служит для преобразования последовательностей в кольце целых чисел. Существенный недостаток метода - зависимость между длиной последовательности N и длиной кодового слова, что приводит к длинным кодовым словам при больших N. 3) Более высокую эффективность и точность обеспечивает метод модульной арифметики. Этот метод также сложен для реализации на ПЛИС, а достижение его потенциального быстродействия на различных процессорах затруднительно. Однако при необходимости достижения высокого быстродействия можно применить подход, который основан на представлении циклической свёртки (2.18) суммой двух линейных сверток vnp[/] и vo6p[/] условно названных соответственно прямой и обратной [78,103]: v[/] = vnp[/] + vo6p[/]. Каждое слагаемое представляет собой различные участки одной и той же линейной свёртки вида: vnp[j] = t,g-4g + j], (4.2) т.е. vo6p[/] = vnp[/ + N], где; = 0.. JV-1. На рис. 4.5 приведена функциональная схема устройства, реализующего косвенный алгоритм. Функциональная схема устройства реализации косвенного алгоритма Рис. 4.5 Устройство состоит из четырех основных блоков: блока вычисления циклической свёртки БВЦС, описываемого (4.2); блока накопления БН, вычисляющего vnop; блока компараторов БК; логического блока ЛБ, л обеспечивающего поиск М по перепаду состояний выходов компараторов 0-Я. Очевидно, этот подход по сравнению с другими методами требует больших вычислительных ресурсов, т.к. использует не менее 2N умножений для вычисления JV-точечной циклической свёртки (например, против (N log2(A0+ ) комплексных умножений для быстрой свёртки на основе БПФ) и выполняется за N тактов после получения всех входных данных. Однако он может быть преобразован к рекурсивной форме. Представим линейную свёртку (4.2) выходом цифрового нерекурсивного фильтра с импульсной характеристикой (ИХ) вида g[k] = {0,\,...,N-l}. (4.3) Обычно при проектировании цифровых рекурсивных фильтров (ЦРФ) по заданной импульсной характеристике применяется прямой метод расчёта фильтров в Z-плоскости [104]. Этот метод обладает тем недостатком, что порядок цифрового рекурсивного фильтра (ЦРФ) равен числу аппроксимируемых отсчётов g[k]. Для сокращения порядка фильтра, используем модифицированную методику расчёта, позволяющую снизить порядок при заданном уровне потерь в точности аппроксимации. Запишем передаточную функцию ЦРФ в виде: H{z) = — = ±h[kyk, ±a[k]z k " где а[к] и р - соответственно коэффициенты и порядок ЦРФ (полагается, что а[0] = 1), h[k] - отсчёты бесконечной импульсной характеристики ЦРФ. Требуется наилучшим образом аппроксимировать заданную ИХ g[k] в пределах первых N отсчётов. Аппроксимируемая ИХ g[k] и передаточная функция ЦРФ связаны соотношением: H(z) = \/a[k]z-k = Ntg[k]z k, которое может быть переписано в виде: G[z\A[z] = \ (4.4) где G[z] = "8[кук, A[z] = a[k]z-k. Обратное Z-преобразование от (4.4) при условии p = N-\ приводит к свёртке между коэффициентами фильтра а[к] и отсчётами аппроксимируемой ИХД ]: р \ 1, т = 0 Ї8[т-к]а[к] = \ (4.5) -о [0, т ї 0 где g[&] = 0 при А; 0. В матричном виде (4.5) может быть записано как: Ga = e, (4.6) где G - нижняя треугольная матрица, составленная из отсчётов g[k], а -искомый вектор коэффициентов ЦРФ, е =[1,0, ... , 0]. Вектор а можно выразить, умножив левую и правую части (4.6) на G" [105]: а = G_1e. (4.7) Поскольку N достаточно велико ( 1024), необходимо снижать порядок фильтра р. Для соответствия числа столбцов матрицы G и размерности вектора а, последний дополняется (N-p) нулями. На основании стробирующего свойства нулей это эквивалентно вычеркиванию в матрице G столбцов с номерами (р+\).. .(ЛЧ) при сохранении числа элементов вектора а равнымр+1, где/? N. Уравнение (4.6) дляр N преобразуется к виду: Смй1 = е, (4.8) где GM - прямоугольная (/?+1)х7У-мерная матрица, полученная из G путем вычеркивания последних N-p столбцов, а і - (р+1)-мерный вектор коэффициентов прир N. Система уравнений (4.8) не имеет точного решения в виде (4.7) и может быть найдена только приблизительно [105], т.к. она является переопределённой. Однако возможно найти решение наилучшее в смысле минимума СКО ИХ фильтра от заданной формы. Введем JV-мерный вектор-столбец невязок w = Ga - GMab характеризующий точность аппроксимации ИХ при уменьшении порядка фильтра. С учётом (4.6) можно записать: w = e-GMai. (4.9) Критерий качества, основанный на (4.9), соответствует минимизации квадрата нормы вектора невязок: wHw-»min, (4.10) где н - знак транспонирования и комплексного сопряжения. Глобальный экстремум ai_opt целевой функции (4.10) является решением системы нормальных уравнений: »,_.,=(с№м)"с;е (4Л1) Вид экстремума характеризуется определенностью матрицы (G GM). Эксперименты подтвердили, что полученное решение является глобальным минимумом, т.к. (G GM) 0 [105]. Заданной ИХ (4.3) в пределах первых N отсчётов точно соответствует реакция на 5-функцию ЦРФ порядка р- 2 с коэффициентами (4.11). Очевидно, что выбор столь малого порядка фильтра позволяет осуществлять вычисление прямой свёртки vnp[/] в темпе поступления данных, что исключает влияние этой части процедуры вычисления циклической свёртки на величину вычислительных затрат.
