概要
レーダ断面積（Radar Cross-Section: RCS）は，レーダによる測定物の散乱特性を特徴づける重要なパラメータです．RCSは測定物の遠方界において定義される量であり，波長と比較して電気長の大きな物体のRCSを遠方界条件を満たす距離で直接測定することは，測定設備の空間的な制約により困難な場合が多くあります．このため，近傍界における測定データを信号処理的に遠方界のRCSへ変換する近傍界遠方界変換（Near-Field-to-Far-Field Transformation: NFFFT）が有効な手法となります．とくに，合成開口レーダ（Synthetic Aperture Radar: SAR）または逆合成開口レーダ（Inverse SAR: ISAR）に基づくNFFFT法 [1]〜[3] は，実装の簡便性と高い変換精度からよく利用される手法です．
画像に基づくNFFFT法では，まず近傍界におけるSARまたはISAR測定によりレーダ画像を再構成し，このレーダ画像から遠方界における散乱界を算出するという二段階の手順をとります．この手法は多くの状況において精度よく遠方界のRCSを計算可能ですが，これまでの研究で形状の非対称性が強い測定物の計測においてRCSの変換誤差が大きくなることが示されていました．そこで文献 [2] では，画像再構成に用いる積分変換において補正係数を導入し，RCSの計算精度の改善を図っています．この補正係数を用いたNFFFT法は，従来は円形走査や円筒走査について個別に定式化がなされてきました．
本稿ではこれらを一般化し，任意のアンテナ走査曲面に対応可能な手法を提案します [3]．提案手法の要点は，補正係数の計算に現れるアンテナ走査曲面の偏微分係数を数値微分により求めることで，一般的なアンテナ走査曲面に対応した補正係数の計算を可能としたことです．簡易的な計算機シミュレーションにより，提案するNFFFT法の有効性を示します．
問題の定式化
本節では問題の定式化を行います．任意曲面においてアンテナ走査を行うモノスタティック散乱測定を考え，測定した受信データの画像化に必要となる焦点化関数と補正係数を導出し，あわせて画像によるRCSの算出について述べます．
システムモデル
図1に，本稿で考える測定系のシステムモデルを示します．三次元空間領域 \((x, y, z)\) における測定を考え，空間領域内の一般的な位置を \(\boldsymbol{r}\) で表します．パラメータ \((u, v)\) で特徴づけられる任意の曲面 \(\boldsymbol{r}_0(u, v)\)（走査曲面）に配置されたアンテナから測定物に対して電磁波を照射し，その散乱波を測定します．しばらくの間は走査曲面を連続関数 \(\boldsymbol{r}_0(u,v)\) として表しますが，後に示すように離散的な三次元位置を用いることができます．これらのベクトル \(\boldsymbol{r}\) と \(\boldsymbol{r}_0(u,v)\) を以下のように定義しておきます．
$$\boldsymbol{r}_0(u,v) = x_0(u,v)\hat{\boldsymbol{x}} + y_0(u,v)\hat{\boldsymbol{y}} + z_0(u,v)\hat{\boldsymbol{z}} \tag{1a}$$
$$\boldsymbol{r} = x\hat{\boldsymbol{x}} + y\hat{\boldsymbol{y}} + z\hat{\boldsymbol{z}} \tag{1b}$$
ここで，\(\hat{\boldsymbol{x}}\)，\(\hat{\boldsymbol{y}}\)，\(\hat{\boldsymbol{z}}\) はそれぞれ \(x\) 方向，\(y\) 方向，\(z\) 方向の単位ベクトルです．また，位置 \(\boldsymbol{r}\) から走査曲面上のアンテナ位置 \(\boldsymbol{r}_0(u,v)\) へ向かうベクトル \(\boldsymbol{R}\) は次のように計算できます．
$$\boldsymbol{R}(u,v) = \boldsymbol{r}_0 – \boldsymbol{r} = R_x\hat{\boldsymbol{x}} + R_y\hat{\boldsymbol{y}} + R_z\hat{\boldsymbol{z}} \tag{2a}$$
$$R_x(u,v) = x_0(u,v) – x \tag{2b}$$
$$R_y(u,v) = y_0(u,v) – y \tag{2c}$$
$$R_z(u,v) = z_0(u,v) – z \tag{2d}$$
ベクトル \(\boldsymbol{R}\) を用いて，空間内の任意の位置 \(\boldsymbol{r}\) と走査曲面上のアンテナ位置 \(\boldsymbol{r}_0(u,v)\) との距離は次式のように表されます．
$$R(u,v) = |\boldsymbol{R}(u,v)| = |\boldsymbol{r}_0 – \boldsymbol{r}| = \sqrt{R_x^2 + R_y^2 + R_z^2} \tag{3}$$
受信信号モデル
測定に用いる角周波数を \(\omega\)，波動の伝搬速度を \(c\) として，対応する波数を \(k = \omega/c\) と表記します．位置 \(\boldsymbol{r}\) にある点散乱体の反射係数を \(C\) とすれば，アンテナ位置 \(\boldsymbol{r}_0\) で受信される散乱波 \(E^s(k,\boldsymbol{r}_0)\) は次のように表されます．
$$E^s(k,\boldsymbol{r}_0) = P^2(k,\boldsymbol{r}_0,\boldsymbol{r})\,\frac{k^2 C}{\sqrt{4\pi}}\,\frac{e^{-2jk|\boldsymbol{r}_0-\boldsymbol{r}|}}{|\boldsymbol{r}_0-\boldsymbol{r}|^2} \tag{4}$$
ここで，\(P(k,\boldsymbol{r}_0,\boldsymbol{r})\) は位置 \(\boldsymbol{r}_0\) から位置 \(\boldsymbol{r}\) へ向かう方向におけるアンテナの指向性を表す関数です．散乱体が半径 \(a\) の小さな導体球（\(ka < 0.4\)）の場合，\(C\) は次式で与えられます．
$$C = 3\sqrt{\pi}\,a^3 \tag{5}$$
なお，複数の散乱体が存在する場合の受信信号は，式(4)をそれぞれの散乱体位置について計算した結果の総和もしくは積分として表されます．式(4)の受信信号モデルは，後述する計算機シミュレーションで用います．
画像再構成
画像再構成の目的は，受信信号 \(E^s(k,\boldsymbol{r}_0)\) から散乱体の位置と反射係数を特定することです．画像再構成が理想的になされた場合，位置 \(\boldsymbol{r}_1 = (x_1, y_1, z_1)\) にある反射係数 \(C_1\) の点散乱体の三次元画像 \(\psi(\boldsymbol{r})\) は次のようになります．
$$\psi(\boldsymbol{r}) = C_1\,\delta(x-x_1)\,\delta(y-y_1)\,\delta(z-z_1) \tag{6}$$
ここで，\(\delta(\cdot)\) はデルタ関数を表します．画像再構成は，一般的に重み関数 \(F(k,\boldsymbol{r}_0,\boldsymbol{r})\) による受信信号 \(E^s(k,\boldsymbol{r}_0)\) の積分変換として次式のように表されます．
$$\psi(\boldsymbol{r}) = \int_0^{\infty}\!\!\int_{\Omega} E^s(k,\boldsymbol{r}_0)F(k,\boldsymbol{r}_0,\boldsymbol{r})\,d\boldsymbol{r}_0\,dk = \int_0^{\infty}\!\!\int_{\Omega_u}\!\!\int_{\Omega_v} E^s(k,\boldsymbol{r}_0)F(k,\boldsymbol{r}_0,\boldsymbol{r})\,du\,dv\,dk \tag{7}$$
ここで，\(\Omega\)，\(\Omega_u\)，\(\Omega_v\) はそれぞれ \(\boldsymbol{r}_0\)，\(u\)，\(v\) に関する積分範囲を表します．また，重み関数 \(F(k,\boldsymbol{r}_0,\boldsymbol{r})\) は焦点化関数（Focusing Factor）と呼ばれ，次式で与えられます [2], [3]．
$$F(k,\boldsymbol{r}_0,\boldsymbol{r}) = g(k,\boldsymbol{r}_0,\boldsymbol{r})\,\frac{|\boldsymbol{r}_0-\boldsymbol{r}|^2}{P^2(k,\boldsymbol{r}_0,\boldsymbol{r})}\,e^{2jk|\boldsymbol{r}_0-\boldsymbol{r}|} \tag{8}$$
ここで，\(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) は補正係数であり，次節で導出を行います．式(4)の受信信号モデルと式(8)を比較するとわかるように，焦点化関数は位置 \(\boldsymbol{r}\) に対応する位相を打ち消すように受信信号を補正します．したがって，焦点化関数を重みとすれば受信信号はコヒーレントに積分され，位置 \(\boldsymbol{r}\) の応答が得られることになります．なお，\(g(k,\boldsymbol{r}_0,\boldsymbol{r}) = 1\) とすれば，従来の合成開口レーダ信号処理で一般的に用いられる形式となります．
補正係数
本節では補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) を導出します．補正係数は，式(7)の計算結果が式(6)になるように決定する必要があります．式(4)に従い，位置 \(\boldsymbol{r}_1\) にある点散乱体の散乱波を次式のように考えます．
$$E_1^s(k,\boldsymbol{r}_0) = P^2(k,\boldsymbol{r}_0,\boldsymbol{r}_1)\,\frac{k^2 C_1}{\sqrt{4\pi}}\,\frac{e^{-2jk|\boldsymbol{r}_0-\boldsymbol{r}_1|}}{|\boldsymbol{r}_0-\boldsymbol{r}_1|^2} \tag{9}$$
式(9)の受信信号を式(7)に代入すると次式が得られます．
$$\psi(\boldsymbol{r}) = \frac{C_1}{\sqrt{4\pi}}\int_0^{\infty}\!\!\int_{\Omega_u}\!\!\int_{\Omega_v} g(k,\boldsymbol{r}_0,\boldsymbol{r})\,\frac{P^2(k,\boldsymbol{r}_0,\boldsymbol{r}_1)}{P^2(k,\boldsymbol{r}_0,\boldsymbol{r})}\,\frac{|\boldsymbol{r}_0-\boldsymbol{r}|^2}{|\boldsymbol{r}_0-\boldsymbol{r}_1|^2}\,e^{j2ks}\,k^2\,du\,dv\,dk \tag{10a}$$
$$s = |\boldsymbol{r}_0-\boldsymbol{r}| – |\boldsymbol{r}_0-\boldsymbol{r}_1| \tag{10b}$$
いま，図2に示すように，位置 \(\boldsymbol{r}\) を原点とする局所座標系 \((x’, y’, z’)\) を考えます．この局所座標系において，ベクトル \(\boldsymbol{R}\) の方位角 \(\alpha(u,v)\) と天頂角 \(\beta(u,v)\) は次式のように表されます．
$$\alpha(u,v) = \tan^{-1}\left[R_y(u,v)/R_x(u,v)\right] \tag{11a}$$
$$\beta(u,v) = \tan^{-1}\left[\rho/R_z(u,v)\right] \tag{11b}$$
$$\rho(u,v) = \sqrt{R_x^2(u,v) + R_y^2(u,v)} \tag{11c}$$
本稿の議論において，逆正接 \(\tan^{-1}(\cdot/\cdot)\) はその分子と分母の符号により \((-\pi, \pi]\) の範囲で計算するものとします．式(11)の角度 \(\alpha\) と \(\beta\) を用いると，散乱体位置 \(\boldsymbol{r}_1\) の近傍の位置 \(\boldsymbol{r}\)（\(\boldsymbol{r}\to\boldsymbol{r}_1\)）について，関数 \(s\) を次のように近似できます [2], [3]．
$$s \approx (x-x_1)\sin\beta\cos\alpha + (y-y_1)\sin\beta\sin\alpha + (z-z_1)\cos\beta \tag{12}$$
式(12)を式(10)に代入して \(\boldsymbol{r}\to\boldsymbol{r}_1\) とすれば次式を得ます．
$$\psi(\boldsymbol{r}) = \frac{C_1}{\sqrt{4\pi}}\int_0^{\infty}\!\!\int_{\Omega_u}\!\!\int_{\Omega_v} g(k,\boldsymbol{r}_0,\boldsymbol{r})\,e^{j2k(x-x_1)\sin\beta\cos\alpha}\,e^{j2k(y-y_1)\sin\beta\sin\alpha}\,e^{j2k(z-z_1)\cos\beta}\,k^2\,du\,dv\,dk \tag{13}$$
次に，以下の変数変換による置換積分を考えます．
$$k_x = 2k\sin\beta\cos\alpha \tag{14a}$$
$$k_y = 2k\sin\beta\sin\alpha \tag{14b}$$
$$k_z = 2k\cos\beta \tag{14c}$$
式(14)は点 \(\boldsymbol{r}\) を原点とする局所座標系における空間周波数です．測定曲面 \(\boldsymbol{r}_0(u,v)\) が測定物を含む画像化領域を完全にとり囲む閉曲面であるとき，局所座標系の方位角 \(\alpha\) の範囲は \([0, 2\pi]\)，天頂角 \(\beta\) の範囲は \([0, \pi]\) ですから，波数 \(k_x\)，\(k_y\)，\(k_z\) による積分区間はそれぞれについて \([-\infty, \infty]\) となります．測定曲面が閉曲面でない場合は，仮想的に測定曲面を延長して画像化領域をとり囲む閉曲面とし，延長した測定曲面における受信信号を零と考えれば同様の積分区間となります．式(14)を式(13)に代入して置換積分すれば次式を得ます．
$$\psi(\boldsymbol{r}) = \frac{C_1}{\sqrt{4\pi}}\int_{-\infty}^{\infty}\!\!\int_{-\infty}^{\infty}\!\!\int_{-\infty}^{\infty} g(k,\boldsymbol{r}_0,\boldsymbol{r})\,\frac{k^2}{|J|}\,e^{jk_x(x-x_1)+jk_y(y-y_1)+jk_z(z-z_1)}\,dk_x\,dk_y\,dk_z \tag{15}$$
ここで，\(J\) は次式で表されるヤコビアンです（本ページ末尾の付録を参照）．
$$J = 8k^2\sin\beta\,(\alpha_v\beta_u – \alpha_u\beta_v) \tag{16}$$
なお，式(16)において以下の定義を用いました．
$$\alpha_u = \frac{\partial\alpha}{\partial u} = \frac{1}{\rho^2}\left(R_x y_{0u} – R_y x_{0u}\right) \tag{17a}$$
$$\alpha_v = \frac{\partial\alpha}{\partial v} = \frac{1}{\rho^2}\left(R_x y_{0v} – R_y x_{0v}\right) \tag{17b}$$
$$\beta_u = \frac{\partial\beta}{\partial u} = \frac{1}{R^2}\left(R_z\rho_u – \rho z_{0u}\right) \tag{17c}$$
$$\beta_v = \frac{\partial\beta}{\partial v} = \frac{1}{R^2}\left(R_z\rho_v – \rho z_{0v}\right) \tag{17d}$$
$$\rho_u = \frac{\partial\rho}{\partial u} = \frac{1}{\rho}\left(R_x x_{0u} + R_y y_{0u}\right) \tag{17e}$$
$$\rho_v = \frac{\partial\rho}{\partial v} = \frac{1}{\rho}\left(R_x x_{0v} + R_y y_{0v}\right) \tag{17f}$$
$$\chi_\xi = \frac{\partial\chi(u,v)}{\partial\xi},\quad \chi\in\{x_0, y_0, z_0\},\quad \xi\in\{u, v\} \tag{17g}$$
ここで目的としているのは，式(15)の計算結果が式(6)となるように補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) を決定することです．したがって，デルタ関数の積分表示
$$\delta(x-x_1)\delta(y-y_1)\delta(z-z_1) = \frac{1}{(2\pi)^3}\int_{-\infty}^{\infty}\!\!\int_{-\infty}^{\infty}\!\!\int_{-\infty}^{\infty} e^{jk_x(x-x_1)+jk_y(y-y_1)+jk_z(z-z_1)}\,dk_x\,dk_y\,dk_z \tag{18}$$
と式(15)を比較すれば，補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) は次式のように選べばよいことがわかります．
$$g(k,\boldsymbol{r}_0,\boldsymbol{r}) = \frac{|J|}{2\pi^{\frac{5}{2}}k^2} = \frac{4}{\pi^{\frac{5}{2}}}\left|\sin\beta\right|\left|\alpha_v\beta_u – \alpha_u\beta_v\right| \tag{19}$$
円筒走査系や球面走査系では，式(17)における \(\chi(u,v)\) を方程式で与えることができるため，偏微分 \(\chi_\xi(u,v)\) を解析的に評価したうえで補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) を求めることができます．しかしながら，ここでは任意の曲面を考えているため，この偏微分は解析的に計算できない場合もあります．そこで次節では，この偏微分を数値微分に置き換える手法を提案します．
数値微分による補正係数の計算
前節で示したように，補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) の計算には，走査曲面 \(\boldsymbol{r}_0(u,v)\) の \(u\) および \(v\) に関する偏微分が必要となります．しかしながら，任意の測定曲面では解析式が得られない場合もあるため，数値微分による偏微分の評価が必要となります．数値微分の手法として最も簡易的なものは差分による近似ですが，より数値的に安定した手法として，以下では測定曲面 \(\boldsymbol{r}_0(u,v)\) を多項式で近似し，その近似多項式の微分係数を数値微分の値として用いる方法を述べます．
図3に示すように，測定曲面のある部分領域が離散的な座標値 \(\boldsymbol{r}_0(u_m, v_n)\) として与えられている場合を考えます．いま，画素位置 \(\boldsymbol{r}\) の補正係数 \(g[k, \boldsymbol{r}_0(u_m,v_n), \boldsymbol{r}]\) を計算するために，図3のように位置 \(\boldsymbol{r}_0(u_m,v_n)\) における測定曲面の微分係数を数値的に求めたいとします．ベクトル \(\boldsymbol{r}_0(u_m,v_n)\) の成分を \(\chi(u_m,v_n),\ \chi\in\{x_0,y_0,z_0\}\) で表記し，またその変数 \(\xi\in\{u,v\}\) に関する偏微分を \(\chi_\xi(u,v)\) で表記して議論を進めます．さらに，表記を簡潔にするため，走査曲面上の座標 \((u_m, v_n)\) におけるベクトル \(\boldsymbol{r}_0\) の成分 \(\chi(u_m,v_n)\) を \(\chi_{m,n}\) と略記します．まず，\(\chi_{m,n}\) を次のように \(u_m\)，\(v_n\) の多項式で展開します．
$$\chi_{m,n} = \sum_{p=0}^{K}\sum_{q=0}^{L} a^{(\chi)}_{p,q}\,u_m^p\,v_n^q \tag{20}$$
ここで，\(a^{(\chi)}_{p,q}\) は次数 \(p\) および \(q\) の展開係数であり，\([\cdot]^{(\chi)}\) はこの係数がベクトル \(\boldsymbol{r}_0\) の成分 \(\chi\) に関するものであることを表します．また，\(K\)，\(L\) はそれぞれ \(u_m\)，\(v_n\) に関する多項式の次数を表します．式(20)は以下のようにベクトルと行列で表記できます．
$$\chi_{m,n} = \boldsymbol{u}_m^T\boldsymbol{A}^{(\chi)}\boldsymbol{v}_n \tag{21a}$$
$$\boldsymbol{u}_m = \left[1, u_m, u_m^2, \ldots, u_m^p, \ldots, u_m^K\right]^T \tag{21b}$$
$$\boldsymbol{v}_n = \left[1, v_n, v_n^2, \ldots, v_n^q, \ldots, v_n^L\right]^T \tag{21c}$$
$$\boldsymbol{A}^{(\chi)} = \begin{bmatrix} a^{(\chi)}_{0,0} & a^{(\chi)}_{0,1} & \cdots & a^{(\chi)}_{0,q} & \cdots & a^{(\chi)}_{0,L}\\ a^{(\chi)}_{1,0} & a^{(\chi)}_{1,1} & \cdots & a^{(\chi)}_{1,q} & \cdots & a^{(\chi)}_{1,L}\\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ a^{(\chi)}_{p,0} & a^{(\chi)}_{p,1} & \cdots & a^{(\chi)}_{p,q} & \cdots & a^{(\chi)}_{p,L}\\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots\\ a^{(\chi)}_{K,0} & a^{(\chi)}_{K,1} & \cdots & a^{(\chi)}_{K,q} & \cdots & a^{(\chi)}_{K,L}\end{bmatrix} = \left[\boldsymbol{a}^{(\chi)}_0, \boldsymbol{a}^{(\chi)}_1, \ldots, \boldsymbol{a}^{(\chi)}_p, \ldots, \boldsymbol{a}^{(\chi)}_K\right]^T \tag{21d}$$
$$\boldsymbol{a}^{(\chi)}_p = \left[a^{(\chi)}_{p,0}, a^{(\chi)}_{p,1}, \ldots, a^{(\chi)}_{p,q}, \ldots, a^{(\chi)}_{p,L}\right]^T \tag{21e}$$
式(21)をさらに変形し，次式のようにベクトルの内積で表現することができます．
$$\chi_{m,n} = (\boldsymbol{u}_m\otimes\boldsymbol{v}_n)^T\left[\boldsymbol{a}^{(\chi)T}_0, \boldsymbol{a}^{(\chi)T}_1, \ldots, \boldsymbol{a}^{(\chi)T}_p, \ldots, \boldsymbol{a}^{(\chi)T}_K\right]^T = \boldsymbol{w}_{m,n}^T\boldsymbol{a}^{(\chi)} \tag{22a}$$
$$\boldsymbol{w}_{m,n} = \boldsymbol{u}_m\otimes\boldsymbol{v}_n \tag{22b}$$
$$\boldsymbol{a}^{(\chi)} = \left[\boldsymbol{a}^{(\chi)T}_0, \boldsymbol{a}^{(\chi)T}_1, \ldots, \boldsymbol{a}^{(\chi)T}_p, \ldots, \boldsymbol{a}^{(\chi)T}_K\right]^T \tag{22c}$$
ここで，\(\otimes\) はKronecker積を表します．これを部分領域内のすべてのサンプル点について考えると，次式の表現が得られます．
$$\boldsymbol{\chi} = \boldsymbol{W}\boldsymbol{a}^{(\chi)} \tag{23a}$$
$$\boldsymbol{\chi} = \left[\boldsymbol{\chi}_1^T, \boldsymbol{\chi}_2^T, \ldots, \boldsymbol{\chi}_m^T, \ldots, \boldsymbol{\chi}_M^T\right]^T \tag{23b}$$
$$\boldsymbol{\chi}_m = \left[\chi_{m,1}, \chi_{m,2}, \ldots, \chi_{m,n}, \ldots, \chi_{m,N}\right]^T \tag{23c}$$
$$\boldsymbol{W} = \left[\boldsymbol{W}_1^T, \boldsymbol{W}_2^T, \ldots, \boldsymbol{W}_m^T, \ldots, \boldsymbol{W}_M^T\right]^T \tag{23d}$$
$$\boldsymbol{W}_m = \left[\boldsymbol{w}_{m,1}, \boldsymbol{w}_{m,2}, \ldots, \boldsymbol{w}_{m,n}, \ldots, \boldsymbol{w}_{m,N}\right]^T \tag{23e}$$
式(23a)の係数ベクトル \(\boldsymbol{a}^{(\chi)}\) は，最小二乗法に基づいて次式から計算できます．
$$\boldsymbol{a}^{(\chi)} = \left(\boldsymbol{W}^T\boldsymbol{W}\right)^{-1}\boldsymbol{W}^T\boldsymbol{\chi} \tag{24}$$
この係数 \(\boldsymbol{a}^{(\chi)}\) を用いれば，小曲面内において次のような多項式展開が行えます．
$$\chi(u,v) = \boldsymbol{u}(u)\boldsymbol{A}^{(\chi)}\boldsymbol{v}^T(v) \tag{25a}$$
$$\boldsymbol{u}(u) = \left[1, u, u^2, \ldots, u^m, \ldots, u^K\right]^T \tag{25b}$$
$$\boldsymbol{v}(v) = \left[1, v, v^2, \ldots, v^n, \ldots, v^L\right]^T \tag{25c}$$
したがって，偏微分 \(\chi_\xi(u,v)\) は次のように求められます．
$$\chi_u(u,v) = \frac{\partial\chi(u,v)}{\partial u} = \boldsymbol{u}'(u)\boldsymbol{A}^{(\chi)}\boldsymbol{v}^T(v) \tag{26a}$$
$$\chi_v(u,v) = \frac{\partial\chi(u,v)}{\partial v} = \boldsymbol{u}(u)\boldsymbol{A}^{(\chi)}\boldsymbol{v}’^T(v) \tag{26b}$$
$$\boldsymbol{u}'(u) = \left[0, 1, 2u, \ldots, mu^{(m-1)}, \ldots, Ku^{(K-1)}\right]^T \tag{26c}$$
$$\boldsymbol{v}'(v) = \left[0, 1, 2v, \ldots, nv^{(n-1)}, \ldots, Lv^{(L-1)}\right]^T \tag{26d}$$
以上の手順をすべての成分 \(\chi\in\{x_0, y_0, z_0\}\) に対して繰り返して偏微分 \(\chi_\xi(u,v)\) を求めれば，式(17)の計算が行えることになり，最終的に式(19)の補正係数を求めることが可能になります．
RCSの算出
文献 [2], [3] で議論されているように，波数 \(k\)，方位角 \(\phi_0\)，天頂角 \(\theta_0\) のRCSは，再構成したレーダ画像 \(\psi(\boldsymbol{r})\) の空間逆Fourier変換として次式のように表せます．
$$\sigma\left[\boldsymbol{k}_r(k,\phi_0,\theta_0)\right] = k^4\left|\mathcal{F}^{-1}_{(\boldsymbol{r})}\left[\psi(\boldsymbol{r})\right]\right|^2 \tag{27}$$
ここで，\(\mathcal{F}^{-1}_{(\boldsymbol{r})}[\cdot]\) は空間変数 \(\boldsymbol{r}\) に関する逆Fourier変換であり，\(\boldsymbol{k}_r(k,\phi_0,\theta_0)\) は次式で定義される波数ベクトルです．
$$\boldsymbol{k}_r(k,\phi_0,\theta_0) = 2k\,\hat{\boldsymbol{r}}_0(\phi_0,\theta_0) \tag{28a}$$
$$\hat{\boldsymbol{r}}_0(\phi_0,\theta_0) = \boldsymbol{r}_0/|\boldsymbol{r}_0| = \sin\theta_0\cos\phi_0\,\hat{\boldsymbol{x}} + \sin\theta_0\sin\phi_0\,\hat{\boldsymbol{y}} + \cos\theta_0\,\hat{\boldsymbol{z}} \tag{28b}$$
式(27)を計算機で実装する場合は，単純に逆離散Fourier変換（Inverse Discrete Fourier Transform: IDFT）を用いるか，空間領域画像 \(\psi(\boldsymbol{r})\) の三次元IFFT（Inverse Fast Fourier Transform）により直交座標系の空間周波数領域 \((k_x, k_y, k_z)\) のスペクトラム画像を求め，この画像を極座標系 \((k, \phi_0, \theta_0)\) に補間して所望の周波数，方位角，仰角におけるRCSを算出すればよいことになります．
提案アルゴリズムのまとめ
これまでに述べたNFFFTのアルゴリズムを以下にまとめます．
- Step 1. 画像化対象とする画素位置 \(\boldsymbol{r} = (x, y, z)\) を決定する．
- Step 2. 決定した画素位置 \(\boldsymbol{r}\) について，式(2)から局所座標系におけるアンテナ位置 \((R_x, R_y, R_z)\) を求め，また式(11c)から \(\rho(u,v)\) を求める．
- Step 3. Step 2の計算結果を用い，式(17)に基づいて局所座標系におけるアンテナの方位角 \(\alpha\) および天頂角 \(\beta\) の偏微分 \(\alpha_u, \alpha_v, \beta_u, \beta_v\) を決定する．これらの偏微分が解析的に評価できない場合は，前述の数値微分に基づいた評価を行う．
- Step 4. Step 3の計算結果を用い，式(19)に基づいて画素位置 \(\boldsymbol{r}\) に対する補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) を計算し，これを基に式(8)の焦点化関数 \(F(k,\boldsymbol{r}_0,\boldsymbol{r})\) を計算する．
- Step 5. Step 4で求めた焦点化関数 \(F(k,\boldsymbol{r}_0,\boldsymbol{r})\) を用い，式(7)に基づいて画素位置 \(\boldsymbol{r}\) の画素値を再構成する．
- Step 6. 上記のStep 1からStep 5までを，すべての画素位置に対して繰り返す．
- Step 7. 式(27)に基づき，再構成した画像 \(\psi(\boldsymbol{r})\) の三次元逆Fourier変換により，所望の波数 \(k\)，方位角 \(\phi_0\)，天頂角 \(\theta_0\) に対するRCS \(\sigma(\boldsymbol{k}_r)\) を求める．
計算機シミュレーション
本節では，点散乱体モデルを用いた簡易的な計算機シミュレーションにより，提案する任意走査曲面に対する補正係数およびNFFFT法の有効性を検証します．
シミュレーション概要
本シミュレーションでは，式(4)および式(5)に基づく点散乱体モデルを用います．また，アンテナ走査曲面として，次式に示す原点を中心とする半径 \(r_0\) の球面を仮定します．
$$x_0(u,v) = r_0\cos u\sin v \tag{29a}$$
$$y_0(u,v) = r_0\sin u\sin v \tag{29b}$$
$$z_0(u,v) = r_0\cos v \tag{29c}$$
次式に示すように，式(29)に対する補正係数 \(g(k,\boldsymbol{r}_0,\boldsymbol{r})\) の計算に必要な偏微分 \(\alpha_u\)，\(\alpha_v\)，\(\beta_u\)，\(\beta_v\) は解析的に評価できます [3]．
$$\alpha_u = \frac{r_0\sin v}{\rho^2}\left(R_x\cos u + R_y\sin u\right) \tag{30a}$$
$$\alpha_v = \frac{r_0\cos v}{\rho^2}\left(R_x\sin u – R_y\cos u\right) \tag{30b}$$
$$\beta_u = \frac{R_z r_0\sin v}{\rho R^2}\left(R_y\cos u – R_x\sin u\right) \tag{30c}$$
$$\beta_v = \frac{r_0}{\rho R^2}\left[R_z\cos v\left(R_x\cos u + R_y\sin u\right) + \rho^2\sin v\right] \tag{30d}$$
なお本シミュレーションでは，式(30)に基づく厳密な偏微分の評価と，前述の数値的な偏微分の評価を両方行い，同一の結果が得られることを確認しています．よって以下では，数値微分による方法の結果のみを示します．
| 項目 | 諸元 |
|---|---|
| 送信アンテナ | 等方性アンテナ |
| 走査球の半径 \(r_0\) | 1 m |
| 中心周波数 | 1 GHz |
| 帯域幅 | 500 MHz |
| 周波数サンプル間隔 | 10 MHz |
| 方位角サンプル間隔 | 2° |
| 方位角範囲 | [0°, 360°] |
| 天頂角サンプル間隔 | 2° |
| 天頂角範囲 | [0°, 180°] |
| 導体球の位置 | (−0.3, 0, −0.3) m (0.3, 0, 0.3) m |
| 導体球の半径 \(a\) | 1.5 cm |
表1にシミュレーション諸元を示します．中心周波数は1 GHz，帯域幅は500 MHzであり，アンテナは全球走査を仮定しました．表1に示す位置に半径1.5 cmの導体球を配置しています．
シミュレーション結果
図4と図5は，それぞれ従来手法（補正係数なし）および提案手法（補正係数あり）を用いたシミュレーション結果を示します．これらの図において，(a)は三次元の空間領域画像，(b)は(a)に示す画像の空間周波数スペクトラム，(c)は \(\theta_0 = \pi/2\) の場合における方位角方向のRCS，(d)は \(\phi_0 = 0\) の場合における天頂角方向のRCSを示します．(a)の画像において，アンテナの走査半径を超える領域の画素値は零としています．式(27)に基づいて(a)に示す三次元画像の逆Fourier変換をIFFTにより計算したものが(b)の空間周波数スペクトラムであり，ここでは \(k_z = 0\) における \((k_x, k_y)\) 平面の断面図と，\(k_y = 0\) における \((k_x, k_z)\) 平面の断面図を表示しています．この空間周波数スペクトラムが(c)と(d)に示すRCSにそのまま対応しますが，ここでは補間誤差を避けるため，(a)の画像からIDFTにより式(27)を評価して(c)と(d)の結果を出力しています．また，(c)と(d)にはRCSの理論値を表示しています．
図4(a)に示す従来手法によるレーダ画像では散乱体の周囲にアーチファクトが生じているのに対し，図5(a)に示す提案手法のレーダ画像では散乱体の位置を中心として対称的な形状に応答が現れており，これらのレーダ画像から提案する補正係数の効果が明確に確認できます．
図4(b)に示す空間周波数スペクトラムの \((k_x, k_z)\) 平面の断面図を見てみると，明らかに \(\theta_0 = 0\) および \(\theta_0 = \pi\) に対応する位置の強度が強くなっていることがわかります．これは本シミュレーションにおいて方位角 \(\phi_0\) と天頂角 \(\theta_0\) のサンプリング間隔を一定とした球面走査を行っており，このため球面の北極点（\(\theta_0 = 0\)）と南極点（\(\theta_0 = \pi\)）の付近で空間サンプル点の密度が高くなるためです．補正係数を用いるとこのサンプル密度の差異が補正されるため，図5(b)の空間周波数スペクトラムでは，すべての方位角と天頂角方向において均一にスペクトラムが分布していることが確認できます．
最終的なRCSの計算結果である図4(c),(d)と図5(c),(d)とを比較すれば，補正係数を導入したことによるRCS算出の改善効果は明らかです．特に，図4(d)に示す補正係数無しの天頂角特性ではRCSの理論値とNFFFT結果が大きく異なっており，これは前述した空間周波数スペクトラムの分布に関する議論により説明できます．図5(c),(d)の補正係数を導入した場合のNFFFT結果はRCSの理論値と良好に一致しており，提案する補正係数の有効性が確認できます．
まとめと今後の予定
本稿では，これまで円筒走査系や球面走査系で定式化されてきたレーダ画像方式の近傍界遠方界変換を一般化し，これらの走査系を包含する任意の走査曲面に対して用いることができる統一的な手法を構築しました．提案手法の要点は，走査曲面の偏微分を数値微分により求めることで，任意の走査曲面に対する補正係数の導出を可能としたことです．簡易的な計算機シミュレーションにより，提案手法の有効性を確認しました．今後は，複雑形状の計算機シミュレーションによる提案理論の検証と，電波暗室内における実験検証を行う予定です．
付録：ヤコビアンの計算
式(15)中のヤコビアン \(J\) は次式で定義されます．
$$J = \frac{\partial(k_x, k_y, k_z)}{\partial(k, u, v)} = \begin{vmatrix} \dfrac{\partial k_x}{\partial k} & \dfrac{\partial k_x}{\partial u} & \dfrac{\partial k_x}{\partial v}\\[6pt] \dfrac{\partial k_y}{\partial k} & \dfrac{\partial k_y}{\partial u} & \dfrac{\partial k_y}{\partial v}\\[6pt] \dfrac{\partial k_z}{\partial k} & \dfrac{\partial k_z}{\partial u} & \dfrac{\partial k_z}{\partial v}\end{vmatrix} \tag{31}$$
上式の各要素は以下のように計算できます．
$$\frac{\partial k_x}{\partial k} = 2\cos\alpha\sin\beta \tag{32a}$$
$$\frac{\partial k_x}{\partial u} = 2k\left(\beta_u\cos\alpha\cos\beta – \alpha_u\sin\alpha\sin\beta\right) \tag{32b}$$
$$\frac{\partial k_x}{\partial v} = 2k\left(\beta_v\cos\alpha\cos\beta – \alpha_v\sin\alpha\sin\beta\right) \tag{32c}$$
$$\frac{\partial k_y}{\partial k} = 2\sin\alpha\sin\beta \tag{32d}$$
$$\frac{\partial k_y}{\partial u} = 2k\left(\beta_u\sin\alpha\cos\beta + \alpha_u\cos\alpha\sin\beta\right) \tag{32e}$$
$$\frac{\partial k_y}{\partial v} = 2k\left(\beta_v\sin\alpha\cos\beta + \alpha_v\cos\alpha\sin\beta\right) \tag{32f}$$
$$\frac{\partial k_z}{\partial k} = 2\cos\beta \tag{32g}$$
$$\frac{\partial k_z}{\partial u} = -2k\beta_u\sin\beta \tag{32h}$$
$$\frac{\partial k_z}{\partial v} = -2k\beta_v\sin\beta \tag{32i}$$
ここで，\(\alpha_u\)，\(\alpha_v\)，\(\beta_u\)，\(\beta_v\) は式(17)で与えられます．
参考文献
- T. Vaupel and T. F. Eibert, “Comparison and application of near-field ISAR imaging techniques for far-field radar cross section determination,” IEEE Trans. Antennas Propag., vol. 54, no. 1, pp. 144–151, Jan. 2006.
- A. Osipov, H. Kobayashi, and H. Suzuki, “An improved image-based circular near-field-to-far-field transformation,” IEEE Trans. Antennas Propag., vol. 61, no. 2, pp. 989–993, Feb. 2013.
- T. Watanabe and H. Yamada, “Far-field radar cross-section determination from near-field 3-D synthetic aperture imaging with arbitrary antenna scanning surfaces,” 2021. [Online]. Available: 10.36227/techrxiv.16802419（後に IEEE Trans. Antennas Propag., vol. 70, no. 7, pp. 5831–5840, July 2022 として出版）