О погрешности расчета импульсной характеристики рассеяния лоцируемого объекта при аппроксимации его поверхности конечными элементами Текст научной статьи по специальности «Физика»

УДК 629.735.05

Ю. Л. Смирнова

О ПОГРЕШНОСТИ РАСЧЕТА ИМПУЛЬСНОЙ ХАРАКТЕРИСТИКИ РАССЕЯНИЯ ЛОЦИРУЕМОГО ОБЪЕКТА ПРИ АППРОКСИМАЦИИ ЕГО ПОВЕРХНОСТИ КОНЕЧНЫМИ ЭЛЕМЕНТАМИ

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

Пусть некоторый объект с оптически однородной диффузно рассеивающей поверхностью облучается и наблюдается в направлении оси OZ из точки, удаленной от объекта на расстояние Ь0, существенно превышающее размеры видимой части объекта. При зондирующем импульсе излучения с длиной волны А и мощностью I(£*) = Е8(Ь*), где Е — энергия импульса, 5(Ь*) — дельта-функция аргумента £* = £ + + 2Ь0/с, £ — время, с — скорость распространения электромагнитных волн в атмосфере для мощности излучения, рассеянного от элементарного участка ¿Б поверхности объекта в направлении источника и приходящегося на единицу телесного угла, получим [1]

где R\ = const — коэффициент рассеяния излучения с длиной волны А*; р — угол между внешней нормалью к dS и направлением на источник излучения; Q — телесный угол выхода излучения источника; z — координата, отсчитываемая в направлении наблюдения от точки видимой части поверхности объекта, ближайшей к источнику излучения.

Поскольку Lo » z, то E Rxf (Q(Lq + z)2) w E Rxf (ftLg) = A = = const. Тогда с учетом выражения (1) для полной мощности излучения, рассеянного всей видимой поверхностью объекта, получим

(1)

s

Согласно данным работы [1], импульсная характеристика рассеяния (ИХР) объекта, имеющая размерность м2/с, определяется как #(£) = п I(£)/А. Поэтому

,2 , f ■

g(t) = J (cos ^)2 - - dS. (2)

s

В том случае, когда видимая часть поверхности объекта является поверхностью вращения вокруг оси OZ, угол ^ можно представить функцией координаты z, т.е. F(z) dz = (cos ^(z))2 ds(z). Тогда, используя замену т = 2z/c, вместо уравнения (2) запишем

ZL

g(t) = f F (z) ¿ft - f)dz =

= 2 J F (?) ' (t - T) dT = |f(|) =£F (z),

0

где тl = 2zL/c, 0 ^ t ^ tl; zL — протяженность видимой части поверхности в направлении оси OZ. Отсюда следует, что функцию F(z) можно считать аналогом импульсной характеристики рассеяния, но с размерностью в метрах, и поэтому для вычисления значения F(z) будем использовать формулу:

F(z) = I; J cos2 p(z) dS(z). (3)

s

В дальнейшем удобнее в качестве положительного направления оси OZ выбрать направление на источник излучения. Рассмотрим кольцевой элемент поверхности, образованный вращением вокруг оси OZ выпуклой дуги с переменным (в общем случае), но конечным радиусом кривизны Ri(z), где z — аппликата точки M этой дуги (рис. 1). Аппликату центра O1 кривизны дуги в точке M обозначим через z0 (z). Радиус кривизны рассматриваемого элемента в окружном направлении будет равен отрезку нормали к дуге, заключенному между осью вращения и этой дугой, т.е. R2(z) = O2M для сечения поверхности плоскостью, перпендикулярной оси OZ и проходящей через точку M. Тогда для облучаемого участка такой поверхности равенство (3) может быть представлено в следующем виде:

ф4

d

F(z) = — 2nR1(z)R2(z) sinф cos2 (4)

dz

ф(г)

Рис. 1. Кольцевой элемент поверхности

где ф(г) — угол между внешней нормалью к дуге в точке М и осью причем ф(г) < ф* < п/2 и

'г - ^о(г)

= arccos

Ri(z)

(5)

Так как для выпуклой дуги функция ф(г), согласно уравнению (5), в пределах облучаемой части поверхности вращения (см. рис. 1) имеет однозначную обратную функцию г(ф) при ф € [0, п/2], то радиусы кривизны и Л2 являются функциями угла ф, тогда перейдем к рассмотрению функций г^ф) = ^(¿(ф)) и г2(ф) = Д2(г(ф)). Таким образом, подынтегральная функция в правой части равенства (4) будет явно зависеть лишь от переменного интегрирования и, как следствие,

d

Ф*

F(-) = 2nd- J ^(ф)г2(ф^т ф cos2 ф^ф.

Используя правило Лейбница дифференцирования интеграла по параметру [2] с учетом уравнения (5) и очевидного равенства

\Jl — (z — z0(z))2/R2(z) = sinф, приходим к следующей цепочке равенств для определения ИХР лоцируемого объекта:

r1 (ф)г2(ф) вЬф cos2 ф

F (-) = —, X

X

y/l - (z - *,(z))2/R2(z) 1 — dz0(z)/dz z — z0(z) dR1(z)

Ri(z)

R2(z)

dz

= 2nR (z) 1 —

dz0(z) dR1(z)

dz

dz

cos ф^ ) ) cos2 ф(z) =

^ , ч ^ , ч d cos ф(z) 2 , . .

= 2nR1(z)R2(z)-d^^ cos2 ф(z). (6)

Для дуги окружности длиной s выражение для определения ИХР ло-цируемого объекта относительно принятой точки наблюдения следует из выражения (6) при замене множителя 2п на значение углового размера as = s/(R2(ф)sinф этой дуги. Заметим, что приведенный результат позволяет рассчитать значения ИХР для участков поверхности вращения с выпуклой криволинейной образующей. Эти значения целесообразно использовать в качестве эталона при оценке погрешности расчета ИХР, вносимой аппроксимацией такой поверхности более простыми в геометрическом отношении участками. В частном случае сферической поверхности R1(z) = R2(z) = R = const функция (5) имеет однозначную обратную функцию z = R cos ф, а из (6) следует известное выражение для ИХР этой поверхности [3]:

F°(z) = 2nR cos2 ф^) = 2nz2/R.

которое далее будет использовано в качестве эталонного.

Поверхность, полученная вращением вокруг оси OZ эллипса с полуосями а и b, задается уравнением (ж2 + y2)/a2 + z 2/b2 = 1 ив точке M(ж*, 0, z*) имеет касательную плоскость, определяемую уравнением жж*/а2 + zz*/b2 = 1 [4] (рис.2). Нормаль к поверхности рассматриваемого эллипсоида вращения в точке M определяется уравнением прямой ж = ж* + (z — z*)(b/a)2x*/z* с угловым коэффициентом

tg фЫ = . (7)

a

Радиус кривизны поверхности в окружном направлении в точке M равен _

I z2a4

Рис. 2. Поверхность вращения эллипса

а в меридиональном направлении он совпадает с радиусом кривизны дуги эллипса в этой же точке и равен [4]

Ri(z.) = «2b2( + ■ (9)

-21 N 2 \ 3/2

a4 + b4

Используя выражения (7)-(9) и равенство xt = a^1 — z2/b2, позволяющее выразить первую координату точки M через ее третью координату из уравнения рассматриваемого эллипсоида вращения, согласно выражению (6), находим ИХР этого эллипсоида при z = z*:

„f ч о (x*(z*)b4 +z2a^2 x*/z* — dx*/dz* 2 i( \

F(Zi) = -a^-1 + (b/a)4(x*/z*)2 Sinф(^ COS ф(^) =

( Z2a4 \ = 2п I X*(z*) +--(* 4 I sin 0(z*) cos2 0(z*). (10)

При z* = b в правой части соотношения (10) возникает неопределенность типа 0/0, так как x*(b) = 0 и sin 0(b) = 0. Формальное раскрытие этой неопределенности позволяет найти F(b) = 2na2/b, что геометрически означает равенство обоих радиусов кривизны значению a2/b и возможность рассмотрения поверхности эллипсоида в малой окрестности точки с аппликатой z* = b как сферической. Действительно, подстановка в равенство F°(R) = 2nR для сферы значения R = a2/b приводит к указанному выше значению для F(b).

Если кольцевой элемент поверхности образован вращением вокруг оси OZ дуги с нулевой кривизной, то Ri ^ то. Так, например, при вращении вокруг оси OZ отрезка прямой MH, где H(0, 0,h) £ OZ, имеем элемент поверхности прямого кругового конуса (рис. 3). При его облучении вдоль оси OZ для расчета ИХР непосредственное применение равенств (6) невозможно, поскольку угол ф между нормалью к его боковой поверхности и осью OZ не зависит от z. Для преодоления возникших трудностей рассмотрим элементарный участок поверхности кольцевого элемента с аппликатой z и площадью dS = 2nR2(z) dz/ sin ф = 2n(H — z)ctgфdz/ sin ф, где H — высота рассматриваемой части конуса (см. рис. 3). Таким образом,

z z

S(z) = i dS = 2n[(H — z) dz = 2n(Hz — z2/2)

sin2 ф sin2 ф

0 0

и согласно формуле (3)

F (z) = 2n(H — z)^. (11)

sin2 ф

X к

И 2

Рис. 3. Поверхность конуса

В случае сложной формы поверхности лоцируемого объекта установление зависимости Б (г) в формуле (3) и нахождение соответствующего интеграла связаны с преодолением значительных трудностей. Поэтому для моделирования ИХР такого объекта целесообразно аппроксимировать его поверхность совокупностью плоских треугольных конечных элементов (КЭ) [3, 5]. Такие КЭ просты и удобны в применении. Из семейства двумерных КЭ они имеют наименьшее число сторон, что облегчает нахождение их ИХР, но вносит неизбежную погрешность при аппроксимации криволинейной поверхности лоци-руемого объекта, особенно в ее областях с малыми значениями главных радиусов кривизны. Аналогичными свойствами обладают также и плоские прямоугольные КЭ, каждый из которых можно рассматривать как совокупность двух плоских треугольных КЭ. Отметим, что в общем случае для уменьшения погрешности аппроксимации можно применить более сложные по форме изопараметрические КЭ, имеющие криволинейную поверхность [6,7].

Оценим погрешность вычисления ИХР сферической поверхности при ее аппроксимации плоскими КЭ. Пусть в прямоугольной декартовой системе координат определена сфера Ф радиусом Я с центром в начале координат (рис. 4). Найдем ИХР плоского квадратного КЭ ш с диагональю I и вершинами А, В, С и Б, принадлежащими Ф, т.е.

А : Р п Рз п Ф, В : Р п Рз п Ф, С : Р п Р4 п Ф, Б : Р п Р4 п Ф,

где Р,, г = 1, 4 — плоскости, определяемые в выбранной системе координат уравнениями

Pi : У

2^2

P2 : У

2^2

P3 : z — zi, P4 : z — Z2,

Рис. 4. Прямоугольный конечный элемент на поверхности сферы

причем 0 < z1 < z2 < R. Площадь такого КЭ равна Бш = /2/2, а аппликата его центра (рис. 5)

z, = v4R2 — /2 ^ = Rv4-a2cos w, 22

где A = //R, w — угол между единичным вектором n внешней нормали к w и осью OZ (см. рис. 4). Аппликаты сторон КЭ, перпендикулярных плоскости XOZ, определяются следующими равенствами:

z1 = R | л/4 — A2 cos w--^ sin w

2 V v^

z2 = R I — A2 cos ю + sin ю

2\ V2

В сечении сферы плоскостью Р1 (или плоскостью Р2) лежит окружность радиуса

r(l) = R cos ю* = \ / R2 —

l

2^

= R

V8—Ä2 2^2 :

2

где a* = arcsin(A/v^8)- Таким образом, согласно формуле (3) и предположению о плоском волновом фронте падающего на сферу излучения, для ИХР любого перпендикулярного оси OZ отрезка в пределах рассматриваемого КЭ ш имеем

Fl = cos2 A = R acos!^. (12)

z2 — z1 v 2 sin a

В силу конечности значения величины l угол a изменяется в пределах от amin до Amax = п/2 - Amin (см. рис.5), где Amin = = arcsin (l/(2v^2R cos a*)) = arcsin(l/V8R2 — l2). При этом наибольшее возможное значение l будет соответствовать значению a = п/4. Для такого КЭ Zi = 0, Z2 = r(l) и lmax/v^ = v/2r(l) = v/8R2—lmaX/2. Отсюда следует, что lmax = 2v/2R/v/3-

Сравним уравнение (12) с ИХР участка П поверхности сферы Ф, ограниченного линиями ее пересечения с плоскостями Д, i = 1, 4 (см. рис. 4). Отметим, что для всех точек элемента сферы, площадью dS = 2R2(sin a) arcsin (l/(2\/2R sin a)) dA, заключенного между плоскостями P3 и P4, угол a между внешней нормалью и осью OZ остается неизменным, причем cos a = z/R. Таким образом, для участка Q*(z) поверхности сферы, заключенного между плоскостями P3, P4 и плоскостями, содержащими ось OY и составляющими с осью OZ углы a(z) = arccos(z/R) и a = п/2, получим

п/2 п

cos2 AdS = 2R2 arcsin I ———:- I cos2 a sin AdA-

2\f2R sin a

А так как в рассматриваемом случае cos a = z/R и

dA(z) d arccos(z/R) 1

dz dz R^/1 — z2/R2 R sin a:

1

то, воспользовавшись равенством (3) и правилом дифференцирования интеграла по переменному нижнему пределу [2], для ИХР множества точек рассматриваемого участка поверхности, принадлежащего плоскости, перпендикулярной оси О^, находим

Относительное отличие ИХР КЭ ш от ИХР соответствующего ему участка Q поверхности сферы можно характеризовать величиной относительной погрешности А = (Fq — Fw)/Fq, где в качестве Fq выберем значение Fq* (z°) при аппликате z° = R cos ф точки пересечения сферической поверхности с нормалью к КЭ, проходящей через его центр.

На рис. 6 приведены результаты расчета относительной погрешности А в зависимости от угла ф £ [0, п/2 — ф*] для различных значений Л £ (0; /max/R), полученные с использованием равенств (12) и (13). Существенное увеличение А с уменьшением ф даже для сравнительно малых значений Л свидетельствует о непригодности квадратных КЭ для аппроксимации сферической поверхности в окрестности ее точки N с аппликатой z = R (см. рис. 4). При приближении к этой точке растет и абсолютная погрешность AFW = (Fq* (z°) — Fw)/(2nR) (рис. 7),

(13)

Д,отн.ед. 0,4

0 0,2 0,4 0,6 0,8 1,0 #(яг/2),отн;<р

Рис. 6. Относительная погрешность вычисления ИХР

Д/^отн.ед. 0,20

0,05

0,15

0,10

0 0,1 0,2 0,3 0,4 0,5 Ф!(х12),отн.ед.

Рис. 7. Абсолютная погрешность вычисления ИХР

нормированная по значению Р°(Я) = 2пЯ ИХР для сферической поверхности в точке N.

Оценим погрешность вычисления ИХР в окрестности точки N с аппликатой г = Я (см. рис.4) при ее аппроксимации вместо квадратных КЭ равными по площади КЭ в виде равнобедренных треугольников с общей вершиной в точке N и остальными вершинами, принадлежащими сферической поверхности и имеющими одинаковую аппликату ¿д. Примем для удобства сопоставления основание каждого треугольника равным //л/2, т.е. стороне рассмотренного ранее квадратного КЭ с диагональю, равной /. Тогда каждому целому числу п > 3 треугольников будет соответствовать вполне определенные значения ¿д = л/Я2 — /2/(8 зт2(п/п)) и угла фд = агсэт ((Я — ¿д)/^д),

где кд = у (Я — ¿д)2 + Я2 — ¿Д — /2/8 — высота треугольника. Но в этом случае его ИХР будет линейной функцией аппликаты г:

На рис.8 приведены графики зависимости ДРД = (Р°(г) — —пРд(г)) / (2пЯ), где г = (Я + ¿д)/2, от Л при фиксированных значениях п.

Ясно, что применение совокупности любых плоских треугольных КЭ с общей вершиной в точке N для аппроксимации сферической поверхности в окрестности этой точки не может устранить погрешность вычисления ИХР именно в этой точке, поскольку ИХР треугольников

/Лд cos2 Фд R — z

z е [R — ¿д, R].

Д/^отн.ед,

»4

уП=6

\ х-чи=12 • >\ л=16 \ \

О ОД 0,2 0,3 0,4 0,5 Л,отн.ед.

Рис. 8. Зависимость абсолютной погрешности вычисления ИХР от Л и n

при z = R обращаются в нуль. Ликвидировать эту погрешность можно путем применения криволинейных треугольных КЭ, получаемых путем изопараметрического отображения рассмотренных треугольников на сферическую поверхность или на параболоид вращения [8]. Другой путь связан с расположением в точке N центра квадратного КЭ с достаточно малой стороной h, нормаль к плоскости которого составляет с осью OZ также малый угол фн, удовлетворяющий условию (h/ sin фн)о^2 фн = 2nR. Однако применение такого КЭ приведет к нарушению осевой симметрии при аппроксимации сферической поверхности в окрестности точки N.

На рис. 9 представлены результаты расчета при различных значениях А нормированных по точному для сферы значению F°(R) = 2nR зависимостей ИХР от относительной аппликаты z/R при аппроксимации сферической поверхности совокупностью КЭ в виде прямоугольных треугольников с различными значениями А = 1/R (здесь I — длина наибольшего катета). Для А = 0,01 в пределах выбранного масштаба графика рассчитанная зависимость ИХР вплоть до малой окрестности точки сферы с аппликатой z = R практически совпадает с точной зависимостью F°(z)/(2nR) = (z/R)2. Для больших значений А награ-фике четко видны скачки значений ИХР при относительных аппликатах z/R, соответствующих границам соседних КЭ, перпендикулярным оси OZ.

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

Рис. 9. Сравнение приближенных значений ИХР с точным значением для сферы

в предположении равномерного распределения случайной величины ф по той части S* = 2nR2(cos фт1п — sin 0min) облучаемой поверхности сферы, в пределах которой изменяется угол ф при фиксированном значении Л. Тогда, с учетом формулы (12) для осредненного значения ИХР КЭ получим

Fш = — í F, dS =

Ш S* J 2nR2(cos 0min - sin 0min)

St

п/2-фтЫ

f RA cos2 ф . П - 40min

sin фаф = RA

J у^тф 4v^2(cos фтщ — sin фтт)

фшт

На рис. 10 с учетом выражения для фтт= агсвт(Л/\/8—Л2) сплошной кривой представлена зависимость Fw/F°(R) от Л для КЭ с усредненным значением ф. Для сравнения штриховой линией отмечена рассчитанная в соответствии с выражением (13) зависимость отнесенного к F°(R) = 2nR точного значения ИХР участка сферической поверхности, покрывающего указанный КЭ. Штрихпунктирной линией показана зависимость от Л относительной погрешности е = (Fq — F/Fq, которая возникает при аппроксимации участков сферической поверхности КЭ с усредненным значением ф. Таким образом, задавшись предельно допустимым значением относительной погрешности е* можно найти значение Л* = е-1(е*), из которого затем определить длину I* = RЛ* стороны КЭ.

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

FjF^Fjh", е. отн.ед.

0,10 0,08 0,06 0,04 0,02

0 0,2 0,4 0,6 0,8 1,0 Л,отн.ед.

Рис. 10. Сравнение относительных погрешностей вычисления ИХР сферы

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

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

СПИСОК ЛИТЕРАТУРЫ

1.Лебедько Е. Г., Порфирьев Л. Р., Хайтун Ф. И. Теория и расчет импульсных и цифровых оптико-электронных систем. - Л.: Машиностроение, 1984. - 192 с.

2. К о р н Г., Корн Т. Справочник по математике для научных работников и инженеров / Пер. с англ. - М.: Наука, 1968. 720 с.

3. Б у р ы й Е. В. Синтез системы распознавания объектов по форме огибающей импульса лазерного излучения при импульсно-периодической локации // Квантовая электроника. - 1998. - Т. 25, № 5. - С. 471-475.

4. Бронштейн И. Н., Семендяев К. А. Справочник по математике для инженеров и учащихся втузов. - М.: Наука, 1986. - 544 с.

5. Б у р ы й Е. В., Зубцов С. А., Петров В. А. Использование ультракоротких импульсов лазерного излучения для получения информации о форме и ориентации лоцируемых объектов // Вестник МГТУ им. Н.Э. Баумана. Сер. "Приборостроение". - 1991. - № 3. - С. 100-113.

6. Зенкевич О., Морган К. Конечные элементы и аппроксимация / Пер. с англ. - М.: Мир, 1986. - 318 с.

7. Сегерлинд Л. Применение метода конечных элементов / Пер. с англ. -М.: Мир, 1979.-392 с.

8. С т р е н г Г., Фикс Д ж. Теория метода конечных элементов / Пер. с англ. - М.: Мир, 1977. - 350 с.

Статья поступила в редакцию 13.02.2006

Смирнова Юлия Леонидовна родилась в 1980 г., окончила МГТУ им. Н.Э. Баумана. Мл. науч. сотр. лаборатории лазерных информационных систем НИИ РЛ МГТУ им. Н.Э. Баумана. Автор 13 научных работ в области математического моделирования лазерных систем.

Yu.L. Smirnova (b. 1980) graduated from the Bauman Moscow State Technical University. Junior researcher of the laboratory of laser information systems of the Research Institute for Radio-electronics of the Bauman Moscow State Technical University. Author of 13 publications in the field of mathematical simulation of laser systems.

В издательстве МГТУ им. Н.Э. Баумана в 2006 г. вышла в свет книга

Козинцев В.И.

Основы импульсной лазерной локации. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2006. - 512 с.

Изложены физические основы импульсной лазерной локации. Приведены сведения об оптических свойствах земной атмосферы, отражающих свойствах земной и морской поверхностей и объектов локации. Описаны эффекты, возникающие при распространении лазерных пучков в атмосфере. Рассмотрены методы расчета лазерных сигналов на трассе с отражением от неровной земной и взволнованной морской поверхностей, от светоотража-телей и от объектов сложной формы. Описаны помехи в системах лазерной локации. Изложены теоретические основы приема лазерных сигналов. Приведены примеры лазерных локационных систем различного назначения и описаны их основные элементы. Содержание учебного пособия соответствует курсу лекций, который читают авторы в МГТУ им. Н.Э. Баумана.

Для студентов технических вузов, обучающихся по направлению "Опто-техника", а также для научных работников и инженеров приборостроительного профиля.

По вопросам приобретения обращаться по тел. 263-60-45; e-mail: [email protected]