JP5699405B2

JP5699405B2 - レーダ受信信号処理装置とその方法

Info

Publication number: JP5699405B2
Application number: JP2010260632A
Authority: JP
Inventors: 英樹丸井; 聡久田中; 智志根; 文彦水谷; 将一和田; 石澤　寛; 寛石澤
Original assignee: NATIONAL UNIVERSITY CORPORATION TOKYO UNIVERSITY OF AGRICULUTURE & TECHNOLOGY; Toshiba Corp
Current assignee: NATIONAL UNIVERSITY CORPORATION TOKYO UNIVERSITY OF AGRICULUTURE & TECHNOLOGY; Toshiba Corp
Priority date: 2010-06-11
Filing date: 2010-11-22
Publication date: 2015-04-08
Anticipated expiration: 2030-11-22
Also published as: JP2012018154A; JP5699404B2; JP2012018155A

Description

本発明の実施形態は、レーダにおいて、目標信号と不要信号が混在するレーダ受信信号から目標信号のみを的確に検出するレーダ受信信号処理装置とその方法に関する。

周知のように、気象レーダは、気象目標のエコーを受信してその受信電力を解析することにより降水量を観測している。また、電磁波のドップラー効果を利用し、レーダ受信信号のドップラー周波数を解析することで風速も観測している。同様に、探知系レーダは飛翔体等のエコーを受信信号処理することにより飛翔体の探知を行っている。

但し、レーダ受信信号には、目標信号だけでなく、大地、山岳、海面のエコー等からなる不要信号が含まれる。不要信号は観測の妨げとなるため、レーダの受信信号処理装置には、その不要信号を除去するために移動目標指示装置（Moving Target Indicator：ＭＴＩ）が用いられる。

但し、ＭＴＩは、目標信号と不要信号の周波数変化の違いを利用して、予め遮断周波数が設定された低域遮断フィルタで不要信号を除去している。そのため、設定した遮断周波数以下の成分は目標信号、不要信号の区別をせずに除去してしまうことになる。

例えば、気象レーダにおいては、気象エコーと不要信号が同程度の周波数変化を持つ場合、的確に不要信号のみを除去することはできないため、降水量、風速を正しく観測することができない。

石原正仁，ドップラー気象レーダの原理と基礎，気象研究ノート第２００号，日本気象協会，pp.1-38，2001

上述の如く、従来のレーダ受信信号処理装置にあっては、目標信号と不要信号が同程度の周波数変化を持つ場合、両者を的確に区別することができず、目標信号まで除去してしまうことがあった。

実施形態によれば、一定周期で繰り返し得られ、気象エコー及び地形エコーを含む気象レーダのレーダ受信信号を周波数領域に変換して周波数スペクトルを取得し、前記一定周期で取得されたレーダ受信信号の周波数スペクトルについて、前記気象エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとするフォン・ミーゼス分布による第１の密度関数及び前記地形エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとする正規分布による第２の密度関数の和である混合密度関数で順次モデリングする処理を、前記混合密度関数のエコー混合比、前記第１の密度関数及び第２の密度関数それぞれのパラメータを変数として繰り返し実行して、前記周波数スペクトルに対する混合密度関数の当て嵌まり具合が最大となるときの前記エコー混合比及び前記第１及び第２の密度関数それぞれのパラメータを最適値として順次推定し、前記フォン・ミーゼス分布に基づいて、前記推定手段で順次推定された前記第１の密度関数のパラメータである電力値、ドップラー周波数、スペクトル幅の最適値から気象エコーの受信電力、平均ドップラー周波数、スペクトル幅を気象エコー情報として算出する。特に、前記パラメータの推定では、前記第２の密度関数を予め観測した地形エコーのデータから決定し、決定された第２の密度関数のパラメータを変数として、統計的性質から前記第２の密度関数が従う事前分布を決定し、決定された事前分布と前記レーダ受信信号の観測による尤度関数から事後分布を算出し、前記事後分布を最大化したときの前記混合密度関数のパラメータを最適値として推定するものとし、前記事後分布の最大化は、前記周波数スペクトルに対する前記混合密度関数の当て嵌まり具合をペナルティ付き尤度関数で評価し、前記ペナルティ付き尤度関数の期待値を最大化するときの前記混合密度関数の各パラメータを算出する、という手順を前記パラメータの値が収束するまで繰り返すものとし、前記ペナルティ付き尤度関数は、前記エコー群の持つ平均ドップラー周波数、スペクトル幅及び受信電力の統計的性質から混合密度関数の各パラメータが従う事前分布を決定し、前記事前分布を前記ペナルティとして尤度関数に付加することで決定するものとし、前記事前分布は、前記地形エコー成分である不要信号の電力値、ドップラー周波数、スペクトル幅それぞれに関するパラメータで表され、これらのパラメータを事前に観測したデータから決定するものとし、前記不要信号の電力値に関するパラメータはベータ分布に従うものとし、前記不要信号のドップラー周波数に関するパラメータは正規分布に従うものとし、前記不要信号のスペクトル幅に関するパラメータはガンマ分布に従うものとする。

実施形態のレーダ受信信号処理装置が適用される気象レーダ信号処理装置の構成を示すブロック図。実施形態のレーダ受信信号処理装置において、レーダ受信信号から目標信号の平均ドップラー周波数、スペクトル幅、受信電力を算出する処理の流れを示すフローチャートである。上記実施形態において、レーダ受信信号の周波数スペクトルを混合密度関数でモデリングする様子を示す波形図である。上記実施形態において、事前の知識を用いた最適な混合密度関数のパラメータ推定方法の処理の流れを示すフローチャートである。上記実施形態において、最適な混合密度関数のパラメータ推定方法の処理の流れを示すフローチャートである。上記実施形態において、レーダ受信信号から降水量、風速を算出する処理の流れの例を示すフローチャートである。

以下、実施の形態について、図面を参照して説明する。

図１は実施形態に係るレーダ受信信号処理装置が適用される気象レーダ信号処理装置の構成を示すブロック図である。図１において、パルス信号発生部１１で発生されたレーダパルス信号は送信装置１２で周波数変換され電力増幅されて、サーキュレータ１３を介してアンテナ１４から空間に向けて放射される。

上記アンテナ１４で受けたレーダ波の反射信号は、サーキュレータ１３を介して受信装置１５に送られる。この受信装置１５は、アンテナ１４で受けた信号を増幅し、ベースバンドに周波数変換するもので、その出力はレーダ受信信号処理装置１６に送られる。

上記レーダ受信信号処理装置１６は、入力されたレーダ受信信号をＡ／Ｄ変換器１６１でデジタルデータに変換し、ＩＱ検波器１６２で複素形式のＩＱデータに変換し、観測処理器１６３で地形エコー成分を除去して気象エコーによる観測データを得る。

図２は上記観測処理器１６３において、レーダ受信信号（ＩＱデータ）から目標信号（気象エコー成分）の平均ドップラー周波数、スペクトル幅、受信電力を算出する処理の流れを示すフローチャートである。

まず、レーダ受信信号を周波数領域に変換することで周波数スペクトルを得る（ステップＳ１１）。次に、得られた周波数スペクトルを混合密度関数でモデリングする（ステップＳ１２）。続いて、混合密度関数のパラメータの最適値を学習により推定し（ステップＳ１３）、推定したパラメータから目標信号の平均ドップラー周波数、スペクトル幅、受信電力を算出し、観測結果として出力する（ステップＳ１４）。

図３は、上記ステップＳ１２において、周波数スペクトルを混合密度関数でモデリングする処理を説明するための波形図であり、Ａは周波数スペクトル、Ｂは混合密度関数を示している。図３において、横軸Ｘはドップラー速度[m/s]、縦軸Ｙは振幅[dB]を示す。まず、周波数スペクトルをモデリングする混合密度関数を与える。混合密度関数は複数エコーそれぞれの密度関数の和で表され、そのパラメータはそれぞれのエコーの混合比とスペクトル形状に関係する値となっている。

上記周波数スペクトルに対応する最適な混合密度関数のパラメータの推定には、事前知識の利用が考えられる。図４に事前の知識を利用して最適な混合密度関数のパラメータを推定する場合の処理の流れを示す。

図４において、まず事前の知識から複数エコーそれぞれの密度関数のパラメータを変数として、統計的性質からそれらが従う事前分布を決定する（ステップＳ２１）。決定した事前分布と、観測による尤度から事後分布を算出し、事後分布を最大化することで最適な混合密度関数のパラメータを得る（ステップＳ２２）。

図５は、上記事後分布を最大化することで混合密度関数の最適なパラメータを推定する方法を示すフローチャートである。まず周波数スペクトルと事前の知識を用いてペナルティ付きの尤度関数を決定する（ステップＳ３１）。次に、決定されたペナルティ付きの尤度関数からペナルティ付き対数尤度関数の期待値を算出する（ステップＳ３２）。続いて、ペナルティ付き対数尤度関数の期待値を最大化するパラメータを算出する（ステップ３３）。ここでパラメータの値が収束したか否かを判断し（ステップＳ３４）、収束しない場合にはステップＳ３２に戻る。すなわち、期待値を最大化したパラメータを用いて、ペナルティ付き対数尤度関数の期待値算出からの手順を再び行い、パラメータの値が収束するまで繰り返す。この収束したパラメータの値をレーダ受信信号の周波数スペクトルをモデリングするのに最適な混合密度関数のパラメータとする。

（実施例１）
図６は上記実施形態を気象レーダ信号処理装置に適用した場合の具体例（実施例１とする）であり、レーダ受信信号から降水量、風速を算出する場合の処理の流れを示すフローチャートである。

図６において、まず、レーダ受信信号を周波数領域に変換し、周波数スペクトルの振幅成分y_iを得る（ステップＳ４１）。次に、得られた周波数スペクトルをモデリングする混合密度関数を与える（ステップ４２）。ここで、不要信号（グランドクラッタのエコー成分）をモデリングする分布として、次式（１）で表される正規分布f₁(x|θ₁) を用いる。式（１）において、x はドップラー速度を表す。また、θ = (τ₁,τ₂,θ₁₁,θ₁₂,θ₂₁,θ₂₂) は混合密度関数のパラメータであり、初期値は任意に設定する。
f₁(x|θ₁) = (2πθ₁₂)^-1/2exp{-(x-θ₁₁)²/(2θ₁₂)} …(1)
また、気象エコーをモデリングする分布として次式（２）で表されるフォン・ミーゼス分布f₂(x|θ₂) を用いる。
f₂(x|θ₂) = exp{θ₂₂ cos(x-θ₂₁)}/{2πI₀(θ₂₂)} …(2)
I₀は０次の修正ベッセル関数を表す。

これらの式(1)，(2)より、周波数スペクトルをモデリングする混合密度関数f(x|θ) を次式(3)とする。
f(x|θ) = τ₁f₁(x|θ₁) +τ₂f₂(x|θ₂) …(3)
続いて、与えた混合密度関数から混合比T_1,i、T_2,iを次式(4)，(5)で算出する（ステップＳ４３）。
T_1,i = τ₁f₁(x_i|θ₁)/{f(x|θ)} …(4)
T_2,i = τ₂f₂(x_i|θ₂)/{f(x|θ)} …(5)
次に、予め観測したデータから事前分布を決定する（ステップＳ４４）。まず、不要信号の電力値に関するパラメータは次式(6)で表されるベータ分布p₁(τ₁)に従うとする。
p₁(τ₁) = {τ₁ ^α-1(1-τ₁)^β-1}/B(α,β) …(6)
ここでB(α,β)はベータ関数を表す。

不要信号のドップラー周波数に関するパラメータは次式(7)で表される正規分布p₂(θ₁₁) に従うとする。
p₂(θ₁₁) = (2πφ₂)^-1/2exp{-(θ₁₁-φ₁)²/(2φ₂)} …(7)
また、不要信号のスペクトル幅に関するパラメータは次式(8)で表されるガンマ分布p₃(θ₁₂) に従うとする。
p₃(θ₁₂) = b^-aθ₁₂ ^a-1 exp(-θ₁₂/b)/Γ(a) …(8)
ここで、Γ(a) はガンマ関数を表す。これら事前分布のパラメータα、β、φ₁、φ₂、a 、b は事前に観測した不要信号のデータから決定する。

次に、式(9) に示すペナルティ付き対数尤度関数の期待値を各パラメータで偏微分することで、各パラメータQ(θ|θ^(t)) の最大化を行う（ステップＳ４５）。
Q(θ|θ^(t)) = Σⁿ _i=1[y_iT_1,i{logτ₁f₁(θ₁|x_i)} + y_iT_2,i{logτ₂f₂(θ₂|x_i)}]
+S₁log{p₁(τ₁)} +S₂log{p₂(θ₁₁)} + S₃log{p₃(θ₁₂)} …(9)
ここで、n はフーリエ変換のサンプル点数を表す。また、(t)はステップＳ４３，Ｓ４５の繰返し回数を表す。各パラメータの更新式を以下に示す。

τ₁ = {Σⁿ _i=1y_iT_1,i + S₁(α-1)}/{Σⁿ _i=1y_i + S₁(α+β-2)} …(10)
τ₂ = 1-τ₁ …(11)
θ₁₁ = {Σⁿ _i=1 y_iT_1,ix_i +S₂(φ₁θ₁₂/φ₂)}/{Σⁿ _i=1 y_iT_1,i + S₂(θ₁₂/φ₂)}
…(12)
θ₁₂=(1-S₃){Σⁿ _i=1 y_iT_1,i (x_i-θ₁₁)²/Σⁿ _i=1 y_iT_1,i }
+S₃{-Σⁿ _i=1 y_iT_1,i+2(a-1)+[{Σⁿ _i=1 y_iT_1,i-2(a-1)}²+8/bΣⁿ _i=1 y_iT_1,i(x_i-θ₁₁)²]^1/2}/{4/b}
…(13)
θ₂₁ = tan^-1{(Σⁿ _i=1y_iT_2,isinx_i)/(Σⁿ _i=1y_iT_2,icosx_i)} …(14)
I₁(θ₂₂)/I₀(θ₂₂) = Σⁿ _i=1 y_iT_2,icos(x_i − θ₂₁)/ Σⁿ _i=1 y_iT_2,i …(15)
以上のように、θ₂₂は式(15)を解くことで得る。また、S₁、S₂、S₃は1か0をとり、不要信号の電力値に関する事前知識を用いる場合はS₁=1、用いない場合はS₁=0である。また、不要信号のドップラー周波数に関する事前知識を用いる場合はS₂=1、用いない場合はS₂=0である。また、不要信号のスペクトル幅に関する事前知識を用いる場合はS₃=1、用いない場合はS₃=0である。

次に、上記の式に従って算出したパラメータが収束しているかどうか判定する（ステップＳ４６）。この判定処理において、混合比算出時のパラメータと、更新後のパラメータの差が1/1000以下の場合、パラメータが収束したと判定する。差がそれよりも大きい場合は、算出したパラメータを用いてステップＳ４３，Ｓ４５の処理を再び行う。

上記パラメータが収束したと判定された場合には、推定したフォン・ミーゼス分布のパラメータから、気象エコーの平均ドップラー周波数、スペクトル幅、受信電力を推定値として算出する（ステップＳ４７）。続いて、推定した気象エコーの平均ドップラー周波数、スペクトル幅、受信電力から、降水量、風速を算出する（ステップＳ４８）。

以上の処理を行うことにより、気象エコーに不要信号としてグランドクラッタによる地形エコーが重なる状況であっても、気象エコー、地形エコーそれぞれの成分を高精度に推定することができるので、気象エコーと地形エコーが混在するデータから降水量、風速を確実に推定することができる。

（実施例２）
ところで、上記実施例１では、受信信号の周波数スペクトルを、気象エコー、地形エコー（グランドクラッタ）をそれぞれ仮定した２つの密度関数の和でモデリングした場合について説明した。このような処理によれば、ある程度の効果が期待できるが、実際の受信信号には気象エコー、地形エコーの他に雑音等の不要成分（以下、不要エコーと称する）が含まれている。そのような受信信号の周波数スペクトルを気象エコー、地形エコーをそれぞれ仮定した２つの密度関数の和でモデリングすることは、気象エコー推定にて誤差を生むことが想定される。

そこで、実施例２として、さらに雑音等の不要エコーを含む３つの密度関数の和でモデリングした場合について、再度図６を参照して説明する。

図６において、まず、レーダ受信信号を周波数領域に変換し、周波数スペクトルの振幅成分y_iを得る（ステップＳ４１）。次に、得られた周波数スペクトルをモデリングする混合密度関数を与える（ステップ４２）。ここでは、不要信号（地形エコー、不要エコー）をモデリングする分布として、次式(16)で表される正規分布f₁(x|θ₁) を用いる。ここで、x はドップラー速度を表す。θ = (τ₁,τ₂,τ₃,θ₁₁,θ₁₂,θ₂₁,θ₂₂) は混合密度関数のパラメータであり、初期値は任意に設定する。
f₁(x|θ₁) = (2πθ₁₂)^-1/2exp{-(x-θ₁₁)²/(2θ₁₂)} …(16)
また、気象エコーをモデリングする分布として次式(17)で表されるフォン・ミーゼス分布f₂(x|θ₂) を用いる。
f₂(x|θ₂) = exp{θ₂₂ cos(x-θ₂₁)}/{2πI₀(θ₂₂)} …(17)
I₀は０次の修正ベッセル関数を表す。

また、上記以外の不要エコーをモデリングする分布として次式(18)で表される一様分布f₃(x|θ₃)を用いる。

f₃(x|θ₃) = 1/(2V_nyq) …(18)
ここでV_nyqはナイキスト速度を表す。

これらの式(16)，(17)，(18)より、周波数スペクトルをモデリングする混合密度関数f(x|θ) を次式(19)とする。
f(x|θ) = τ₁f₁(x|θ₁) +τ₂f₂(x|θ₂) +τ₃f₃(x|θ₃) …(19)
続いて、与えた混合密度関数から混合比T_1,i、T_2,i、T_3,iをそれぞれ次式(20)，(21)，(22)で算出する（ステップＳ４３）。
T_1,i = τ₁f₁(x_i|θ₁)/{f(x_i|θ)} …(20)
T_2,i = τ₂f₂(x_i|θ₂)/{f(x_i|θ)} …(21)
T_3,i = τ₃f₃(x_i|θ₃)/{f(x_i|θ)} …(22)
次に、予め観測したデータから事前分布を決定する（ステップＳ４４）。まず、不要信号の電力値に関するパラメータは次式で表されるベータ分布p₁(τ₁)に従うとする。
p₁(τ₁) = {τ₁ ^α-1(1-τ₁)^β-1}/B(α,β) …(23)
ここでB(α,β)はベータ関数を表す。

不要信号のドップラー周波数に関するパラメータは次式で表される正規分布p₂(θ₁₁) に従うとする。
p₂(θ₁₁) = (2πφ₂)^-1/2exp{-(θ₁₁-φ₁)²/(2φ₂)} …(24)
また、不要信号のスペクトル幅に関するパラメータは次式で表されるガンマ分布p₃(θ₁₂) に従うとする。
p₃(θ₁₂) = b^-aθ₁₂ ^a-1 exp(-θ₁₂/b)/Γ(a) …(25)
ここで、Γ(a) はガンマ関数を表す。これら事前分布のパラメータα、β、φ₁、φ₂、a 、b は事前に観測した不要信号のデータから決定する。

次に、式(26) に示すペナルティ付き対数尤度関数の期待値を各パラメータで偏微分することで、各パラメータQ(θ|θ^(t)) の最大化を行う（ステップＳ４５）。
Q(θ|θ^(t)) = Σⁿ _i=1[y_iT_1,i{logτ₁f₁(θ₁|x_i)}+ y_iT_2,i{logτ₂f₂(θ₂|x_i)}
+ y_iT_3,i{logτ₃f₃(θ₃|x_i)}]
+S₁log{p₁(τ₁)} +S₂log{p₂(θ₁₁)} + S₃log{p₃(θ₁₂)} …(26)
ここで、n はフーリエ変換のサンプル点数を表す。また、(t)はステップＳ４３，Ｓ４５の繰返し回数を表す。

各パラメータの更新式を以下に示す。
τ₁ = {Σⁿ _i=1y_iT_1,i + S₁(α-1)}/{Σⁿ _i=1y_i + S₁(α+β-2)} …(27)
τ₂ = {Σⁿ _i=1y_iT_2,i - S₁τ₁Σⁿ _i=1y_iT_2,i}/{Σⁿ _i=1y_i(T_2,i+S₁T_3,i)}
…(28)
τ₃ = 1-τ₁-τ₂ …(29)
θ₁₁ = {Σⁿ _i=1 y_iT_1,ix_i +S₂(φ₁θ₁₂/φ₂)}/{Σⁿ _i=1 y_iT_1,i + S₂(θ₁₂/φ₂)}
…(30)
θ₁₂=(1-S₃){Σⁿ _i=1 y_iT_1,i (x_i-θ₁₁)²/Σⁿ _i=1 y_iT_1,i }
+S₃{-Σⁿ _i=1 y_iT_1,i+2(a-1)+[{Σⁿ _i=1 y_iT_1,i-2(a-1)}²
+8/bΣⁿ _i=1 y_iT_1,i(x_i-θ₁₁)²]^1/2}/{4/b} …(31)
θ₂₁ = tan^-1{(Σⁿ _i=1y_iT_2,isinx_i)/(Σⁿ _i=1y_iT_2,icosx_i)} …(32)
I₁(θ₂₂)/I₀(θ₂₂) = Σⁿ _i=1 y_iT_2,icos(x_i − θ₂₁)/ Σⁿ _i=1 y_iT_2,i …(33)
ここでθ₂₂は式(33)を解くことで得られる。また、S₁ 、S₂ 、S₃ は1か0をとり、不要信号の電力値に関する事前知識を用いる場合はS₁ = 1、用いない場合はS₁ =0である。また、不要信号のドップラー周波数に関する事前知識を用いる場合はS₂ =1、用いない場合はS₂ =0である。また、不要信号のスペクトル幅に関する事前知識を用いる場合はS₃ =1であり、用いない場合はS₃ =0である。

以上の処理を行うことにより、気象エコーに不要信号としてグランドクラッタによる地形エコー、雑音等の不要エコーが重なる状況であっても、気象エコー、地形エコー、不要エコーそれぞれの成分を高精度に推定することができるので、気象エコーと地形エコー及び不要エコーが混在するデータから降水量、風速を確実に推定することができる。

上記実施形態は気象レーダの場合を想定して説明したが、通常のレーダ装置、すなわちレーダ受信信号から不要信号を除去して目標信号を得る場合にも同様に実施可能である。また、アクティブフェーズドアレイレーダのように、送信ビームを電子走査させるレーダ装置やデジタルビームフォーミング技術を用いた複数の受信ビーム形成を行うレーダ装置を用いる場合にも同様に実施可能である。したがって、上記の実施形態によれば、同程度の周波数変化を持つ目標信号と不要信号が混在するレーダ受信信号から、目標信号のみを的確に検出するレーダ受信信号処理装置を提供することができる。

尚、上記実施形態はそのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせでもよい。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。

１１… パルス信号発生部、１２…送信装置、１３…サーキュレータ、１４…アンテナ、１５…受信装置、１６…レーダ受信信号処理装置、１６１…Ａ／Ｄ変換器、１６２…ＩＱ検波器、１６３…観測処理器。

Claims

一定周期で繰り返し得られ、気象エコー及び地形エコーを含む気象レーダのレーダ受信信号を周波数領域に変換して周波数スペクトルを取得する取得手段と、
前記一定周期で取得されたレーダ受信信号の周波数スペクトルについて、前記気象エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとするフォン・ミーゼス分布による第１の密度関数及び前記地形エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとする正規分布による第２の密度関数の和である混合密度関数で順次モデリングする処理を、前記混合密度関数のエコー混合比、前記第１の密度関数及び第２の密度関数それぞれのパラメータを変数として繰り返し実行して、前記周波数スペクトルに対する混合密度関数の当て嵌まり具合が最大となるときの前記前記エコー混合比及び前記第１及び第２の密度関数それぞれのパラメータを最適値として順次推定する推定手段と、
前記フォン・ミーゼス分布に基づいて、前記推定手段で順次推定された前記第１の密度関数のパラメータである電力値、ドップラー周波数、スペクトル幅の最適値から前記レーダ受信信号に含まれる気象エコーの受信電力、平均ドップラー周波数、スペクトル幅を気象エコー情報として算出する気象エコー情報算出手段と
を具備し、
前記推定手段では、
前記第２の密度関数を予め観測した地形エコーのデータから決定し、決定された第２の密度関数のパラメータを変数として、統計的性質から前記第２の密度関数が従う事前分布を決定し、決定された事前分布と前記レーダ受信信号の観測による尤度関数から事後分布を算出し、前記事後分布を最大化したときの前記混合密度関数のパラメータを最適値として推定するものとし、
前記事後分布の最大化は、前記周波数スペクトルに対する前記混合密度関数の当て嵌まり具合をペナルティ付き尤度関数で評価し、前記ペナルティ付き尤度関数の期待値を最大化するときの前記混合密度関数の各パラメータを算出する、という手順を前記パラメータの値が収束するまで繰り返すものとし、
前記ペナルティ付き尤度関数は、前記エコー群の持つ平均ドップラー周波数、スペクトル幅及び受信電力の統計的性質から混合密度関数の各パラメータが従う事前分布を決定し、前記事前分布を前記ペナルティとして尤度関数に付加することで決定ものとし、
前記事前分布は、前記地形エコー成分である不要信号の電力値、ドップラー周波数、スペクトル幅それぞれに関するパラメータで表され、これらのパラメータを事前に観測したデータから決定するものとし、
前記不要信号の電力値に関するパラメータはベータ分布に従うものとし、
前記不要信号のドップラー周波数に関するパラメータは正規分布に従うものとし、
前記不要信号のスペクトル幅に関するパラメータはガンマ分布に従うものとすることを特徴とするレーダ受信信号処理装置。
さらに、前記気象エコー情報として算出された平均ドップラー周波数、スペクトル幅及び受信電力から降水量、風速の気象情報を算出する気象情報算出手段を具備することを特徴とする請求項１記載のレーダ受信信号処理装置。
一定周期で繰り返し得られ、気象エコー及び地形エコーを含む気象レーダのレーダ受信信号を周波数領域に変換して周波数スペクトルを取得し、
前記一定周期で取得されたレーダ受信信号の周波数スペクトルについて、前記気象エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとするフォン・ミーゼス分布による第１の密度関数及び前記地形エコーのスペクトル形状を表し、電力値、ドップラー周波数、スペクトル幅をパラメータとする正規分布による第２の密度関数の和である混合密度関数で順次モデリングする処理を、前記混合密度関数のエコー混合比、前記第１の密度関数及び第２の密度関数それぞれのパラメータを変数として繰り返し実行して、前記周波数スペクトルに対する混合密度関数の当て嵌まり具合が最大となるときの前記エコー混合比及び前記第１及び第２の密度関数それぞれのパラメータを最適値として順次推定し、
前記フォン・ミーゼス分布に基づいて、前記推定手段で順次推定された前記第１の密度関数のパラメータである電力値、ドップラー周波数、スペクトル幅の最適値から前記レーダ受信信号に含まれる気象エコーの受信電力、平均ドップラー周波数、スペクトル幅を気象エコー情報として算出する方法であって、
前記パラメータの推定では、
前記第２の密度関数を予め観測した地形エコーのデータから決定し、決定された第２の密度関数のパラメータを変数として、統計的性質から前記第２の密度関数が従う事前分布を決定し、決定された事前分布と前記レーダ受信信号の観測による尤度関数から事後分布を算出し、前記事後分布を最大化したときの前記混合密度関数のパラメータを最適値として推定するものとし、
前記事後分布の最大化は、前記周波数スペクトルに対する前記混合密度関数の当て嵌まり具合をペナルティ付き尤度関数で評価し、前記ペナルティ付き尤度関数の期待値を最大化するときの前記混合密度関数の各パラメータを算出する、という手順を前記パラメータの値が収束するまで繰り返すものとし、
前記ペナルティ付き尤度関数は、前記エコー群の持つ平均ドップラー周波数、スペクトル幅及び受信電力の統計的性質から混合密度関数の各パラメータが従う事前分布を決定し、前記事前分布を前記ペナルティとして尤度関数に付加することで決定するものとし、
前記事前分布は、前記地形エコー成分である不要信号の電力値、ドップラー周波数、スペクトル幅それぞれに関するパラメータで表され、これらのパラメータを事前に観測したデータから決定するものとし、
前記不要信号の電力値に関するパラメータはベータ分布に従うものとし、
前記不要信号のドップラー周波数に関するパラメータは正規分布に従うものとし、
前記不要信号のスペクトル幅に関するパラメータはガンマ分布に従うものとすることを特徴とするレーダ受信信号処理方法。
さらに、前記気象エコー情報として算出された平均ドップラー周波数、スペクトル幅及び受信電力から降水量、風速の気象情報を算出することを特徴とする請求項３記載のレーダ受信信号処理方法。