Научная статья на тему 'Анализ корреляционной функции отраженного сигнала космического прецизионного радиовысотомера'

Анализ корреляционной функции отраженного сигнала космического прецизионного радиовысотомера Текст научной статьи по специальности «Физика»

CC BY
93
43
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ПРЕЦИЗИОННЫЙ РАДИОВЫСОТОМЕР / PRECISION RADAR ALTIMETER / МОРСКАЯ ПОВЕРХНОСТЬ / SEA SURFACE / КОРРЕЛЯЦИОННАЯ ФУНКЦИЯ / CORRELATION FUNCTION / АППРОКСИМАЦИЯ / APPROXIMATION / ИНТЕРВАЛ КОРРЕЛЯЦИИ / CORRELATION INTERVAL / СРЕДНЕКВАДРАТИЧЕСККАЯ ОРДИНАТА МОРСКИХ ВОЛН / ШИРОКОПОЛОСНЫЙ ЗОНДИРУЮЩИЙ СИГНАЛ / WIDEBAND SIGNAL / SIGNIFICANT WAVE HEIGHT

Аннотация научной статьи по физике, автор научной работы — Баскаков А.И., Мин-Хо Ка, Важенин Н.А., Гришечкин Б.Ю.

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Баскаков А.И., Мин-Хо Ка, Важенин Н.А., Гришечкин Б.Ю.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Analysis of the correlation function of reflected signals for a space borne precision radar altimeter

Analysis is presented in the article of the reflected signal correlation function for the precision radar altimeter at the use of probing signal with linear frequency modulation and possible small deviation of the antenna pattern axis from the vertical. It appears to be rather helpful for the synthesis and analysis of optimum units designed for the procession of precision radar altimeter signals reflected from the sea surface. Correlation interval for the fast fluctuations of reflected signal is defined.

Текст научной работы на тему «Анализ корреляционной функции отраженного сигнала космического прецизионного радиовысотомера»

Радиолокация и радионавигация

УДК 621.396.966

А. И. Баскаков

Московский энергетический институт (технический университет)

Мин-Хо Ка

Корейский политехнический университет

Н. А. Важенин Московский авиационный институт Б. Ю. Гришечкин

ФГУП "Российский научно-исследовательский институт космического

приборостроения"

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

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

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

При проведении геофизических исследований океанографический прецизионный радиовысотомер (ПРВ) космического базирования в сочетании с другими приборами дистанционного зондирования способен давать информацию для широкого круга задач [1], [2]: уточнения морского геоида; картирования гравитационных аномалий; контроля уровня поверхности (приливов, отливов, штормовых нагонов); контроля морских течений; определения высоты морских волн; контроля ледяного покрова; контроля скорости ветра и многих других процессов в Мировом океане, связанных с изменением уровня морской поверхности (МП). ПРВ способен также давать полезную информацию для изучения рельефа суши, уровня воды крупных внутренних водоемов и т. д. Для эффективного решения указанных задач необходима высокая точность измерения высоты, при которой средне-квадратическая ошибка не превышает единиц сантиметров. Такая высокая точность обусловлена тем, что регистрируемые перепады уровня поверхности чрезвычайно малы.

© Баскаков А. И., Мин-Хо Ка, Важенин Н. А., Гришечкин Б. Ю., 2007

53

Отсюда следует необходимость изучения статистических характеристик отраженных радиосигналов. Для их анализа прежде всего определяется корреляционная функция (КФ), необходимая для нахождения алгоритма оптимальной обработки отраженных сигналов и расчета точностных характеристик ПРВ. КФ при использовании в ПРВ зондирующих радиосигналов с линейной частотной модуляцией (ЛЧМ) определена в [3], однако без учета возможных отклонений оси диаграммы направленности антенны (ДНА) ПРВ от вертикали. Кроме того, полученное выражение очень громоздко, так что использовать его для синтеза и анализа оптимальных устройств обработки затруднительно. Поэтому задачи данной статьи -учесть возможные незначительные отклонения оси ДНА от вертикали, определить интервал корреляции (ИК) быстрых флуктуаций отраженных от морской поверхности сигналов ПРВ и проанализировать возможность удобной аппроксимации выражения для КФ при использовании широкополосных зондирующих ЛЧМ-сигналов с высокой разрешающей способностью. При этом желательно представить выражение для КФ в виде произведения двух функций с разделяющимися временными переменными, что позволит существенно упростить задачу синтеза и анализа оптимальных устройств обработки отраженных от МП сигналов ПРВ.

КФ радиосигналов ПРВ, отраженных от морской поверхности. Запишем импульсный сигнал, излученный передатчиком ПРВ по направлению нормали к среднему уровню

поверхности моря, как £(г) = 2//уУ(г)ехр ()}, где УУ (г) = Уа (г)ехр[(г)]

- комплексный закон модуляции; Уа (г), уа - амплитуда и фаза сигнала соответственно.

МП считаем крупношероховатой, т. е. нормально распределенные ординаты морских волн и радиусы их кривизны много больше длины радиоволны X. Используя феноменологическую модель морской поверхности в виде совокупности элементарных отражателей, парциальные сигналы от которых независимы, имеют случайную амплитуду, задержку и равномерно распределенную в пределах 0...2п случайную фазу, представим сигнал на входе приемника ПРВ суперпозицией парциальных сигналов от облучаемой поверхности моря. При этом отраженный сигнал представляет собой нормальный случайный процесс, его математическое ожидание (когерентная составляющая) согласно принятой модели равно нулю. Чтобы полностью охарактеризовать статистику этого сигнала, достаточно найти его КФ. Пусть начало координат совмещено с проекцией фазового центра антенны на средний уровень МП и движется вместе с ПРВ со скоростью КА. В настоящей статье учтем возможность отклонения оси ДНА от вертикали и запишем общее выражение для КФ отраженного от МП сигнала, прошедшего согласованный фильтр приемника [3]:

Р (8, Ф )

Я (12 ) = Яе

2До^Одао до

(4п)3 9ф

р [^ - 21 (8, Ф)/с] р[г2 - 21 (8, Ф)/с] ё8ёф

(1)

¿4 (8, ф )

где До - излучаемая передатчиком мощность; Од - коэффициент усиления антенны (передающую и приемную антенны считаем совмещенными); од - удельная эффективная

площадь отражения морской поверхности для нормально падающей радиоволны; 9 -угол, отсчитываемый от нормали к среднему уровню МП (от оси z); ф - угол на плоско-

======================================Известия вузов России. Радиоэлектроника. 2007. Вып. 1

сти, совпадающей со среднем уровнем МП, отсчитываемый от оси х; F (9, ф) - нормированная функция, учитывающая ДНА и диаграмму обратного рассеяния (ДОР) МП; L (9, ф)

- расстояние до элемента МП в зоне облучения; р (•) - сигнал на выходе согласованного фильтра приемника, совпадающий с точностью до постоянного коэффициента с автокорреляционной функцией (АКФ) зондирующего сигнала; c - скорость света.

Для случая зондирования вблизи вертикали нормированная диаграмма антенны, работающей на излучение и прием, аппроксимируется гауссовской кривой (9) =

= ехр (-5.5502/92), причем F (9, ф) можно представить выражением [4]:

F (9, ф) = ехр [-5.55 (э2ь-92 )/9д + 11.1-99^008 ф/9д ] ехр (-1§2 9/ а?), (2)

где 95ь - угол отклонения оси ДНА от вертикали; 9о - ширина ДНА по половинной

22

мощности, не превосходящая нескольких угловых градусов; аг = 2ад - параметр шероховатости поверхности, равный удвоенной дисперсии наклонов морских волн ад, определяющий ширину ДОР МП Л9дор ~ л/Паг.

Окончательное выражение для Я (д, ¿2) в аналитическом виде можно получить, используя зондирующий ЛЧМ-сигнал с гауссовской огибающей

им и) = ехр [- (%2 - ]Ъ) 12 ], (3)

где % = >/П/ Г5 (Г5 - длительность огибающей импульса на уровне 0.46); Ъ = лД/5/ Г5 (А/5

- ширина спектра сигнала). АКФ этого сигнала [3]:

р (т) = ,¡0-5411 + (Ъ2/ %4) ехр {-0.5%2 [1 + (ъV%4)] т2}. Введя замены переменных ¿1-Т0 = I и ¿2 -^0 = ^ + т (т0 = 2Н/с; Н - высота относительно среднего уровня морской поверхности), получим

R (t, т) = Re4

PqT^G^OQCD

64п2 H 3A/s

0.5 -Ф

exp ( - уюот ) exp

_ 2 ( Dg )2 + 2 ( Dga )2 ( 2t + т ) Dg

2 "I Г / 2 "|2 2 1 + ( Dg а ) т + v/ ( Dg ) ( Dg )

-v t

2 Dg

2^1 + 2 (Dga )2

exp (-5.559У )

(4)

где В = Д/5Г5 ж 1 + Ъ2/ %4 - коэффициент сжатия;

• 2 / „4

V = С

с/(н^2) (V2 = a2©2/|s.55ar2 [l -(З.ЗЗ©^/©2 )] + ©g});

* 2

Ф(х) = (l/Vn) je~y dy интеграл вероятности; а = 2аz/c (аz - среднеквадратическая ор-

дината морских волн).

2

х

4

х

Выражение (4) справедливо при 0s^ < 0.2590. Множитель exp(-5.556;?^/0°) в этом выражении учитывает уменьшение мощности при отклонении оси ДНА от вертикали.

Воспользуемся представлением (4) в виде R (t, т) = r (t, т) cos («от), где r (t, т) - огибающая КФ, которая при т = 0 определяется усредненной мощностью сигнала P (t) = r (t, 0), являющейся информационным сигналом ПРВ [4]. P(t) характеризует усредненную форму импульса, отраженного от МП, прошедшего согласованный фильтр и квадратичный детектор приемника:

P (t) = Po} 2G0*ocD exp 64п2 H 3A/s

( 2 2 \ - ( , \

V ц -Vt 0.5-Ф \>ц t

1 4 1 2 ц)

exp (-5.550°h/0? ), (5)

где p

4

2a2 +1/

n

(Afs )2

,

на форму фронта отраженного сигнала.

Как видно из (5), P (t) является функцией времени, что объясняется нестационарностью случайного процесса отражения зондирующего сигнала от МП. На рис. 1 показан рассчитанный по (5) участок переднего фронта, нормированный к величине P0\2G^a0c^(64п2H3Afs ) exp(-5.55e2h/00 ) при H = 800 км ; / = 320 МГц; 0о = = 0.035 рад, 0sh = 0 и различной степени взволнованности морской поверхности: 1 -

az = 0.5 м, = 0.044 рад2; 2 - az = 1 м, = 0.074 рад2; 3 - az = 4 м, = 0.141 рад2 .

Момент t = 0 соответствует отражению от среднего уровня поверхности.

Оценка высоты до среднего уровня морской поверхности и среднеквадратических ординат морских волн производится, соответственно, по положению и по форме фронта усредненной мощности сигнала P (t).

На рис. 2 показан рассчитанный по (5) участок переднего фронта нормированной усредненной огибающей мощности PH (t) при H = 800 км ; A/s = 320 МГц; 00 = 0.035 рад ,

22

слабом волнении моря (az = 0.5 м, ar = 0.044 рад ) и различной величине угла отклонения оси ДНА от вертикали: 1 - 0sh = 0; 2 - 0sh = 0.1500; 3 - 0sh = 0.2500. Из приведенных зависимостей следует, что кроме уменьшения мощности отраженного сигнала, определяемого коэффициентом exp(-5.550°^/0°), изменяется и крутизна среза P(t) .

0.5

0.5

-60 -40 -20 0 20 40 Рис. 1

60 t, нс -50 0

50 100 Рис. 2

150 200 t, нс

Зависимость г (т) от выбранных параметров ПРВ, режима облучения и степени взволнованности морской поверхности удобно анализировать при трехмерном представлении этой функции, а также при одновременном ее представлении в виде кривых равного уровня. На рис. 3-6 показаны трехмерные (а) и равноуровневые (б) представления начального участка (переднего фронта) огибающей КФ г (^, т) для Н = 800 км; 00 = 0.035 рад ; = 0, различной ширины спектра зондирующего сигнала Д/, различной степени взволнованности морской поверхности а2, аг : рис. 3 - А/^ = 80 МГц;

г

0.75 0.5 0.25

0

Л нс

нс

40

т, нс 40 20

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

0.848

0.727 0.606 0.485 0.364 0.242

г = 0.121

т, нс

Рис. 3

г

0.75 0.5 0.25 0

нс

¿ , нс

40

0.779 0.668 0.557 0.445 0.334 0.223

г = 0.111

т, нс 40 20

20

40

т, нс

Рис. 4

г

0.75 0.5 0.25 0

\ \ \ \ \0

т, нс 15 10 5

¿ , нс

40

0.858

0.735

0.613

0.49

0.368

0.245

г = 0.123

10

15

т, нс

Рис. 5

0

б

а

0

0

б

а

0

5

б

а

r

0.75 0.5 0.25

0

t, нс

40

0.845 0.724 0.604 0.483 0.362 0.241

r = 0.121

т, нс 15 10 5

10

15

т, нс

Рис. 6

2 2 2 2 az = 0.5 м, a2 = 0.044 рад2 ; рис. 4 - À/s = 80 МГц, az = 4 м, af = 0.141 рад2 ; рис. 5 -

A/g = 320 МГц ; a z = 0.5 м, af = 0.044 рад2; рис. 6 - A/s = 320 МГц, a z = 1 м,

a2 = 0.074 рад2 .

Анализ огибающей КФ r ( t, т ) показал следующее:

• функция r (t, т ) вытянута вдоль оси t и симметрична по оси т для всех t, соответствующих максимуму этой функции и ее спадающей в сторону увеличения t части;

• передний (нарастающий) участок функции r ( t, т ) может отклоняться от оси t, но это отклонение тем меньше, чем сильнее неравенство az > 0.5c/A/g ;

• при выполнении условия 0sh < 0.2500 влияние угла отклонения оси ДНА от вертикали 0S на огибающую КФ r (t, т) проявляется только в уменьшении ее амплитуды и крутизны среза ее спадающей вдоль оси t части;

• для широкополосных зондирующих ЛЧМ-сигналов с разрешающей способностью в 1.. .3 нс выражение для КФ можно существенно упростить, представив его в виде произведения двух отдельных функций с разделяющимися переменными t и т . Указанное упрощение полезно для синтеза и анализа оптимальных устройств обработки отраженных от МП сигналов ПРВ.

Последнее утверждение можно подтвердить аналитически. Из (1) с учетом (2) и (3), взяв интеграл по ф и перейдя к полярным координатам, получим

R (tb t2 ) = Re

nP\C

H1

OO

exp(-5.5502h/©2 )exp [ j&0 (t1 -12)] Jexp(-cy/hV2)x

xln

-, +00 * 4CYH(11.10g^/02)] J W(n)p(t1 -T0-Y-n)P(t2-T0-Y-n)dndY

-oo

(6)

где р = 2Р0Х20)ао/(4п)3 ; /0 (•) - функция Бесселя нулевого порядка от мнимого аргумента; Ж (•) - гауссовский закон распределения случайных ординат морских волн по оси 2 ; ц = -2г/с; у = ^2/сН (- текущий радиус облучаемой на поверхности области).

58

0

5

б

а

0

Заменив, как и ранее, ^ - Т0 = ¿2 - Т0 = 7 + т, осуществив усреднение по ансамблю парциальных отраженных сигналов (внутренний интеграл в (6)) и аппроксимируя функцию Бесселя экспонентой [5] с учетом выполнения условия 0^ < 0.2500, из (6) получим

Я (т ) = Яе {[п^св/(2Н 3 Д/5Ь )] ехр (-5.550.2h/0О ) ехр ( - Ую0 т ) И ( т ) х

х ) |^ехр ( —уу) и (7 — у) и (7 + т — у) ё у}, (7)

где И(т) = ехр -(я/2) /¡{Г2 (°2/Н-2)] - функция, отличающаяся от АКФ ЛЧМ-сигнала толь-

2 / 2

ко множителем а / ц в показателе экспоненты. Этот множитель

о2/ц2 = я/

_(5Ь/ о , )2 + 2я], (8)

где 5Ь = с/ (2Д/) - интервал разрешения по дальности зондирующего ЛЧМ-сигнала.

Обычно для широкополосных зондирующих ЛЧМ-сигналов условие а2 > 0.5с// выполняется, поэтому отношением 5^о2 в (8) можно пренебречь, и зависимость функции И (т) от состояния взволнованности морской поверхности практически не проявляется. При этом в (7) и (7) = ехр [-0.572/ц2 ] .

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

да

и(7 + т) = — [ и(ю)ехр[ую(7 + т)]ёю = — [ и(ю)ехр(ую7)(1 + уют)ёю = и (7) + ти'(7) .

—да —да

Приведенное соотношение справедливо, если с ростом ю величина |и(ю)| убывает быстрее, чем |1/ю|, и интеграл сходится. Тогда выражение (7) можно представить в виде

{2 2

^ О°?°СВ ехр(-5.5502ь/02 )ехр(-т)И(т)х 64п2 Н 3Д/5

|ехр (-уу) |и (7 - у)|2 ёу + т |ехр (-уу)и (7 - у)и' (7 - у)ёу 0 0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(9)

где функция ехр (-уу) изменяется очень медленно по сравнению с и (7) и второе слагаемое в (9) близко к нулю. Поэтому окончательно имеем

Я(7,т) = Яе{ехр(-ую0т)И(т)Р(7)}, (10)

где Р (7) определено в (5).

Соотношение (10) выполняется тем точнее, чем сильнее неравенство а2 > 0.5с/Д/5. На рис. 7 показано трехмерное и контурное представления начального участка (переднего фронта) огибающей КФ г (7, т), рассчитанные по (10) при тех же исходных данных, что и

г

0.75 0.5 0.25

0

t, нс

40

т, нс 15 10

а

_1_

-0.845 0.724 -0.604 -0.483 -0.362 -0.241

г = 0.121

10

15

т, нс

б

Рис. 7

для рис. 6. По оси t отложено 120 нс, по оси т отложено 20 нс. Видно, что полученные результаты практически совпадают.

ИК быстрых флуктуаций радиосигналов, отраженных от морской поверхности.

Используя огибающую КФ г (t, т), определим ИК быстрых флуктуаций отраженного сигнала. Поскольку в рассматриваемом случае процесс является нестационарным, ИК будет зависеть от времени. Определим его из выражения [6]:

00

Дтк 0) = [2г ^,0 )] -1 | г ^, т ) ёт .

(11)

Нетрудно показать, что для зондирующего ЛЧМ-сигнала (3) из (11) получим выра-

жение

Лтк ^) = т^ехр

АЛ

V

s у

0.5 + Ф< t ^ 2/ (пЛ/52 )_ + 2а2

)/[ 2/ (пЛ/2 )_ + 2а2 2

0.5 + Ф< t v^/ У (пЛ/2 )_ + 2а2

№ (пЛЛ2 )_ + о 2 2 + 2а

. (12)

На рис. 8 приведены зависимости нормированного ИК Л 0) = Дтк 0) Д/*5 для различных случаев взволнованности МП при Н = 800 км ; Д/5 = 320 МГц; 00 = 0.035 рад ;

05Ь < 0.2500 : 1 - а2 = 0.5 м, а2 = 0.044 рад2; 2 - а2 = 1 м, а2 = 0.074 рад2 ; 3 - а2 = 4 м,

а2 = 0.141 рад2.

Анализ полученных результатов показал, что в режиме с ограничением облучаемой на МП области интервалом разрешения зондирующего сигнала (что обычно имеет место в ПРВ) зависимость ИК от времени проявляется только в начальный момент формирования переднего фронта отраженного сигнала, а далее выполняется равенство Дтк = 1/А/8, что соответствует интервалу разрешения сжатого радиоимпульса на выходе оптимального приемного тракта ПРВ. При выполнении условия < 0.2500 ИК практически не зависит

5

0

5

—00

2

от угла отклонения оси ДНА . В режиме с ограничением облучаемой на МП области

л

'( 4пД/д2 )

- б0

- 30

0

Рис. S

30

t, нс

ДНА ПРВ коэффициент exp

в (12) характеризует увеличение ИК по сравнению с 1/ A/s из-за уменьшения облучаемой на МП области. Однако такой режим в современных ПРВ не используется из-за требуемых больших размеров антенной системы и жестких требований на ориентацию оси ДНА.

На основании полученных результатов можно сделать следующие выводы.

1. Полученное аналитическое выражение для КФ отраженных от МП радиосигналов ПРВ позволяет:

• исследовать зависимость КФ от режима облучения с учетом возможных отклонений оси ДНА от вертикали, выбранных параметров ПРВ и степени взволнованности МП;

• проанализировать вид усредненной мощностной огибающей, являющейся информационным сигналом ПРВ.

2. Анализ огибающей КФ дает возможность определить особенности поведения функции r (t, т ) (см. рис. 3-7).

3. Получено аналитическое выражение для ИК быстрых флуктуаций отраженного от МП сигнала радиовысотомера, и показано, что зависимость ИК от времени проявляется только в начальный момент формирования переднего фронта отраженного сигнала, а далее он определяется шириной спектра зондирующего сигнала Атк = 1/Afs, что соответствует интервалу разрешения зондирующего ЛЧМ-сигнала.

Библиографический список

1. Davis C. H. Satellite radar altimetry II IEEE Trans. Microwave Theory Techn. 1992. Vol. MTT-40, № б. P. 1070-107б.

2. Robert H. S. Introduction to physical oceanography I Texas A&M University. Texas, USA, 2003. 352 p.

3. Баскаков А. И. Исследование возможности использования сигналов с линейной частотной модуляцией для оценки взволнованности морской поверхности II Изв. вузов. Радиофизика. 1978. Т. XXI, № 5. С. 710-713.

4. Baskakov A. I., Vazhenin N. A., Morozov K. N. Comparison of information signals processed using time-and frequency-domain methods in oceanographic precision radar altimeters II Earth observation and remote sensing. 2000. Vol. 1б. P. 449-455.

5. Левин Б. Р. Теоретические основы статистической радиотехники. В 2 кн. Кн. 2-я. М.: Сов. радио, 19б8. 504 с.

6. Цветков Э. И. Основы теории статистических измерений. Л.: Энергия, 1979. 288 с.

3

A. I. Baskakov

Moscow power engineering institute

Ka Min-Ho

Korea polytechnic university

N. A. Vazhenin

Moscow aviation institute

B. Y. Grishechkin

FSUE "Russian institute of space device engineering"

Analysis of the correlation function of reflected signals for a space borne precision radar altimeter

Analysis is presented in the article of the reflected signal correlation function for the precision radar altimeter at the use ofprobing signal with linear frequency modulation and possible small deviation of the antenna pattern axis from the vertical. It appears to be rather helpful for the synthesis and analysis of optimum units designed for the procession of precision radar altimeter signals reflected from the sea surface. Correlation interval for the fast fluctuations of reflected signal is defined.

Precision radar altimeter, sea surface, correlation function, approximation, correlation interval, significant wave

height, wideband signal

Статья поступила в редакцию 26 мая 2006 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.