宇宙からの地球観測第７章 観測対象物の位置の決定

宇宙からの地球観測　第７章の章末問題は位置合わせの入門編です。衛星から撮影したデータが地上のどこにあたるのか。スマホと違い自分自身のGPSでポイントを決めれば良いわけではないのでちょっと想像するだけでも大変そうです。

問題7.1.1

ある時刻の人工衛星は $x_0 = 0 \text{km}$、$y_0 = 0 \text{km}$、$H = 500 \text{km}$ を地面に平行に 5 $\text{km/s}$ の速さで飛行している。

1. 光学センサを搭載し、 $\theta = 30 \text{度}$、 $\phi = 10 \text{度}$ で地上点を観測した。この位置 $(x, y)$ を計算しなさい。

基本公式理解

光学衛星では基本的に向いている方向が対象・観測点なのでシンプルに高度と角度から三角関数を計算するのみです。

本文図7.4 で衛星の位置と対象物の位置　それにたいする方位角$\phi$ 直下からの仰角$\theta$としたときの計算式が下記のように説明されています。
式 (7.6)

$$\mathbf{a} = \begin{pmatrix} \sin \theta \sin \phi \\ \sin \theta \cos \phi \\ -\cos \theta \end{pmatrix} r=\frac{H}{cos \theta}$$

式 (7.7)

$$\mathbf{r}'_p = \begin{pmatrix} H \sin \phi \tan \theta \\ H \cos \phi \tan \theta \end{pmatrix}$$

回答7.1.1

これは実際にPythonで3Dを描画してみるとよく分かります。

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

# パラメータ問題設定より
H = 500  # 人工衛星の高度（km）
theta = np.radians(30)  # 仰角（度）
phi = np.radians(10)  # 方位角（度）
#r = H / np.cos(theta)

# 地上点の座標
x = H * np.tan(theta) * np.sin(phi)
y = H * np.tan(theta) * np.cos(phi)
print(f"GroundTarget x:{x:.2f}, y:{y:.2f}")

# 3Dプロットの設定
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 人工衛星の位置
satellite_pos = np.array([0, 0, H])
ground_pos = np.array([x, y, 0])
origin = np.array([0, 0, 0])

# 人工衛星の位置をプロット
ax.scatter(*satellite_pos, color='red', label='Satellite (H)')
ax.text(*satellite_pos, '  Satellite', color='red')

# 地上点の位置をプロット
ax.scatter(*ground_pos, color='blue', label='Ground Point (r_p)')
ax.text(*ground_pos, '  Ground Point', color='blue')

# 原点をプロット
ax.scatter(*origin, color='green', label='Origin')
ax.text(*origin, 'Origin', color='green')

# 衛星の飛行方向を矢印で
ax.quiver(0,-0.5*y, H, 0, 1.8*y, 0, color='black', arrow_length_ratio=0.1)
ax.quiver(0,-0.5*y, 0, 0, 1.8*y, 0, color='black', arrow_length_ratio=0.1)


# 線を描画
ax.plot([0, 0], [0, 0], [0, H], color='black', linestyle='-')
#ax.plot([0, x], [0, 0], [H, H], color='black', linestyle='--')
ax.plot([0, x], [0, 0], [0, 0], color='black', linestyle='-')

# 衛星から、地上の点
ax.plot([0, x], [0, y], [H, 0], color='black', linestyle='--')
ax.plot([0, x], [0, y], [0, 0], color='black', linestyle='--')

# ラベルとタイトル
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Plot of Satellite and Ground Point')

# 凡例を追加
ax.legend()

# プロットを表示
plt.show()

回答はx:50.13km, y:284.29km です。

問題7.1.2

ある時刻の人工衛星は $x_0 = 0 \text{km}$、$y_0 = 0 \text{km}$、$H = 500 \text{km}$ を地面に平行に 5 $\text{km/s}$ の速さで飛行している。
2. レーダを搭載し、 $r = 600 \text{km}$、 $f_d = 1000 \text{Hz}$ を計測した。この位置 $(x, y)$ を計算しなさい。なお、 $f_0 = 1.0 \times 10^9 \text{Hz}$、 $c = 300000 \text{km/s}$ とする。

SAR衛星の場合は第六章で勉強したドップラー周波数と、対象までの距離と光速の関係を使います。
基本的には２つの式で得られる２つの曲線の交点を求めればよいということになります。

ドップラー周波数と、距離の公式

式 (7.13)

$$f_d = \frac{2}{\lambda} \mathbf{v} \cdot \frac{\mathbf{r_s} - \mathbf{r_p}}{|\mathbf{r_s} - \mathbf{r_p}|}$$

式 (7.14)

$$r=\frac{c}{2}t$$

上記２つの交点が地上の点になり、さらに衛星の直下の地上の点を原点とした直交座標系でかんがえると以下のように変形できる。

式 (7.15)

$$f_d = \frac{2 v y}{\lambda \sqrt{x^2 + y^2 + H^2}}$$

式 (7.16)

$$r^2 = x^2 + y^2 + H^2$$

組み合わせると
式 (7.17)

$$f_d = \frac{2 v y}{\lambda r}$$(7.16)

高度、波長、ドップラー周波数、速度がわかれば座標が計算できることが分かりました。
式 (7.18)

$$r'_p = \begin{pmatrix} \sqrt{ r^2 \left( 1 - \frac{f_d \lambda}{2v} \right)^2 - H^2 } \\ \frac{f_d \lambda r}{2v} \end{pmatrix}$$

ドップラー周波数を描画するための変形

ドップラー周波数と高度がわかった時
本文式 (7.13) と (7.15) をもとに、 Y を求める式に変形する手順を示します。

$$f_d\sqrt{x^2 + y^2 + H^2} = \frac{2 v y}{\lambda }$$

$${f_d^2x^2 + f_d^2y^2 + f_d^2H^2} = \frac{4 v^2 y^2}{\lambda ^2 }$$

$${\lambda ^2f_d^2x^2 + \lambda ^2f_d^2y^2 + \lambda ^2f_d^2H^2} = 4 v^2 y^2$$

$$y^2 (4v^2 - \lambda ^2f_d^2) = {\lambda ^2f_d^2x^2 + \lambda ^2f_d^2H^2}$$

$$y = \sqrt{\frac{ \lambda ^2f_d^2 (x^2 + H^2)}{4v^2 - \lambda ^2f_d^2}}$$

回答7.1.2

SARの場合も3Dで描画しつつ、回答を計算します。

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

# パラメータ問題設定より
H = 500  # 人工衛星の高度（km）
v = 5  # 人工衛星の速度（km/s）
f_d = 1000  # ドップラー周波数（Hz）　
f0 = 1e9 # (Hz)　 
c = 300000  # 光速（km/s）
r=600 # 対象物までの距離 (km)問題設定より

lambda_ = c/f0  # 波長（km）
rg = np.sqrt((r**2) - (H**2)) # 地上の半径


# 位置
r_s = np.array([0, 0, H])
r_0 = np.array([0, 0, 0])

#式7.18 より
x = np.sqrt((r**2)*(1-f_d*lambda_/(2*v))**2-H**2)
y = (f_d * lambda_ * r) / (2*v)
ground_pos = np.array([x, y, 0])
print(f"GroundTarget x:{x:.2f}, y:{y:.2f}")


# 3Dプロットの設定
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 人工衛星の位置をプロット
ax.scatter(*r_s, color='red', label='Satellite (r_s)')
ax.text(*r_s, '  Satellite', color='red')

# 原点をプロット
ax.scatter(*r_0, color='green', label='Origin')
#ax.text(*r_0, '  Origin', color='green')

# 地上点の位置をプロット
ax.scatter(*ground_pos, color='blue', label='Ground Point (r_p)')
ax.text(*ground_pos, '  Ground Point', color='blue')


# 衛星の飛行方向を矢印で
ax.quiver(0,-50, H, 0, 150, 0, color='black', arrow_length_ratio=0.1)
ax.quiver(0,-50, 0, 0, 150, 0, color='black', arrow_length_ratio=0.1)


# rの円を描画
u = np.linspace(0, 2 * np.pi, 100)
x_circle = rg * np.cos(u)
y_circle = rg * np.sin(u)
z_circle = np.zeros_like(u)
ax.plot(x_circle, y_circle, z_circle, color='purple', label='r = (c * t)/2')

# ドップラー周波数fdの曲線を描画
      
fdlambda = (f_d * lambda_) ** 2
x_circle_fd = np.linspace(-r, r, 100)  
y_circle_fd = np.sqrt(fdlambda * (x_circle_fd ** 2 + H ** 2) / (4 * v ** 2 - fdlambda))
z_circle_fd = np.zeros_like(x_circle_fd)
print(y_circle_fd[90])
ax.plot(x_circle_fd, y_circle_fd, z_circle_fd, label=f'f_d = {f_d} Hz')



# 線を描画
ax.plot([0, 0], [0, 0], [0, H], color='black', linestyle='-')
#ax.plot([0, x], [0, 0], [H, H], color='black', linestyle='--')
ax.plot([0, x], [0, 0], [0, 0], color='black', linestyle='-')

# 衛星から、地上の点
ax.plot([0, x], [0, y], [H, 0], color='black', linestyle='--')
ax.plot([0, x], [0, y], [0, 0], color='black', linestyle='--')


# ラベルとタイトル
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_ylim(-50, 100)
ax.set_zlabel('Z')
ax.set_title('3D Plot of SAR Satellite and Ground Point with Flight Direction')

# 凡例を追加
ax.legend()

# プロットを表示
plt.show()

回答は x:297.87km, y:18.00km です
3Dを描きながら試行錯誤すると自分の数式理解もはかどります。

問題7.2

1. 式 (7.11) から標準偏差の式 (7.12) を求めなさい。
2. $r = 700 \text{km}$、 $dr = 0.01 \text{km}$、 $d = 0.01 \text{度}$、 $e = 0.01 \text{度}$ のときの地上位置決定誤差を求めなさい。

基本公式の理解

式 (7.11) から式 (7.12) を導くためにまず各変数の意味を理解します。

式 (7.11): 位置誤差

$$\Delta r'_p = dr \bar{a} + r (f \delta + g \epsilon)$$

この式は、測定誤差 ($dr$、$\delta$、$\epsilon$) が位置決定にどのように影響するかを示しています。それぞれの変数は以下のような意味を持ちます。

• $\Delta r'_p$: 位置誤差ベクトル
• $dr$: 距離の誤差
• $\bar{a}$: 視線ベクトルの真値
• $r$: 地上点から衛星までの距離
• $f \delta$ と $g \epsilon$: 仰角 ($\delta$) と方位角 ($\epsilon$) の誤差に対応する補正ベクトル

式 (7.12): 誤差の標準偏差

$\sigma_{\Delta}$ は、位置決定誤差の標準偏差です。これは、実際の位置と測定された位置の間の誤差のばらつきを示します。具体的には、以下のような意味があります。

1. 位置決定の信頼性: $\sigma_{\Delta}$ が小さいほど、位置決定の信頼性が高いことを示します。これは、測定された位置が実際の位置に近いことを意味します。
2. 測定の精度: $\sigma_{\Delta}$ は測定の精度を示します。標準偏差が小さいほど、測定の精度が高いことを意味します。

$$\sigma_{\Delta} = \sqrt{\sigma_{dr}^2 + r^2 (\sigma_{\delta}^2 + \sigma_{\epsilon}^2)}$$

• $\sigma_{dr}$: 距離誤差 $dr$ の標準偏差。距離測定の誤差のばらつきを示します。
• $\sigma_{\delta}$: 仰角誤差 $\delta$ の標準偏差。仰角測定の誤差のばらつきを示します。
• $\sigma_{\epsilon}$: 方位角誤差 $\epsilon$ の標準偏差。方位角測定の誤差のばらつきを示します。
• $r$: 地上点から衛星までの距離。角度誤差が距離に比例して影響を与えることを示しています。

1. 距離誤差成分 ($\sigma_{dr}^2$):

   • 距離誤差 $dr$ の標準偏差の2乗です。距離測定の誤差の影響を示します。

2. 角度誤差成分 ($r^2 (\sigma_{\delta}^2 + \sigma_{\epsilon}^2)$):

   • 角度誤差 $\delta$ と $\epsilon$ の標準偏差の2乗に、距離 $r$ の2乗を掛け合わせたものです。これは、角度誤差が距離に比例して位置決定に影響を与えることを示します。

回答7.2.1

1. 平均値周りの展開

式 (7.11) を微小項に展開し、まず各誤差成分について2乗項を取り、それぞれの期待値を計算します。

2. 誤差の標準偏差の計算

誤差の標準偏差を計算するために、誤差成分の2乗平均 ($\sigma^2$) を考慮します。

まず、各誤差成分について2乗平均を考えます。

$$(dr \bar{a})^2 = (dr)^2 (\bar{a})^2$$

ここで、$\bar{a}$ は視線ベクトルの真値なので、$(\bar{a})^2$ は1として扱います。
したがって、距離誤差成分の2乗平均は以下のようになります。

$$(dr \bar{a})^2 = (dr)^2$$

補足説明：視線ベクトルの真値の２乗がなぜ１として扱えるのかについて

単位ベクトルの定義

• 長さ（ノルム）が1であるベクトル。

視線ベクトルの真値（$\bar{a}$）が単位ベクトルであるということは、$\bar{a}$ の各成分の2乗和が1であることを意味します。これは、ベクトルの大きさ（長さ）が1であることを示します。

視線ベクトル $\bar{a}$ は次のように定義されています。

$$\bar{a} = \begin{pmatrix} \sin \theta \sin \phi \\ \sin \theta \cos \phi \\ -\cos \theta \end{pmatrix}$$

このベクトルの長さ（ノルム）は次のように計算されます。

$$|\bar{a}| = \sqrt{(\sin \theta \sin \phi)^2 + (\sin \theta \cos \phi)^2 + (-\cos \theta)^2}$$

これを展開すると、

$$|\bar{a}| = \sqrt{\sin^2 \theta \sin^2 \phi + \sin^2 \theta \cos^2 \phi + \cos^2 \theta}$$

ここで、$\sin^2 \phi + \cos^2 \phi = 1$ という三角関数の基本的な恒等式を使用します。

$$|\bar{a}| = \sqrt{\sin^2 \theta ( \sin^2 \phi + \cos^2 \phi ) + \cos^2 \theta}$$

$$|\bar{a}| = \sqrt{\sin^2 \theta \cdot 1 + \cos^2 \theta}$$

$$|\bar{a}| = \sqrt{\sin^2 \theta + \cos^2 \theta}$$

ここで、再び三角関数の基本的な恒等式 $\sin^2 \theta + \cos^2 \theta = 1$ を使用します。

$$|\bar{a}| = \sqrt{1}$$

$$|\bar{a}| = 1$$

したがって、視線ベクトルの真値 $\bar{a}$ のノルムは1です。これは、$\bar{a}$ が単位ベクトルであることを意味し、その2乗和が1になる理由です。

物理的な意味
物理的には、視線ベクトルの真値 $\bar{a}$ は方向を示すために使用されます。単位ベクトルであることで、その方向だけが考慮され、長さ（距離）は1に正規化されます。これにより、計算が簡略化され、方向の情報だけに集中することができます。

次に、角度誤差成分について考えます。

$$r (f \delta + g \epsilon)$$

これを2乗すると、

$$[r (f \delta + g \epsilon)]^2 = r^2 (f \delta + g \epsilon)^2$$

角度誤差成分の2乗平均は、$f \delta$ と $g \epsilon$ の期待値を考えます。

$$(f \delta + g \epsilon)^2 = f^2 \delta^2 + 2 f g \delta \epsilon + g^2 \epsilon^2$$

期待値を取ると、いまガウス分布を想定しているので、$\delta \epsilon$ の期待値はゼロと仮定します。（誤差は＋ー同じ分だけ正規部分するはずという想定なの、それをすべて足し合わせると期待値は平均=0となる）

$$\mathbb{E}[(f \delta + g \epsilon)^2] = \mathbb{E}[f^2 \delta^2] + \mathbb{E}[g^2 \epsilon^2]$$

$$\mathbb{E}[f^2 \delta^2] = f^2 \sigma_{\delta}^2$$

$$\mathbb{E}[g^2 \epsilon^2] = g^2 \sigma_{\epsilon}^2$$

補足説明：分散 ($\sigma^2$) の定義

ランダム変数 $X$ の分散 $\sigma_X^2$ は次のように定義されます：

$$\sigma_X^2 = \mathbb{E}[X^2] - (\mathbb{E}[X])^2$$

しかし、ここでは $\delta$ の期待値が（ガウス分布によると本文で設定があるので） 0 なので、分散は単に次のようになります：

$$\sigma_{\delta}^2 = \mathbb{E}[\delta^2]$$

したがって、角度誤差成分の期待値は以下のようになります。

$$\mathbb{E}[(f \delta + g \epsilon)^2] = f^2 \sigma_{\delta}^2 + g^2 \sigma_{\epsilon}^2$$

さらに $f,g$の下記定義

$${f} = \begin{pmatrix} \cos \theta \sin \phi \\ \cos \theta \cos \phi \\ -\sin \theta \end{pmatrix}$$

$${g} = \begin{pmatrix} \sin \theta \cos \phi \\ -\sin \theta \sin \phi \\ 0 \end{pmatrix}$$

ここで少し困ってしまいました。
$f^2$ もaと同様に方向を示す単位ベクトルのノルムとして捉えられるので1になります。
ただ、 $g^2$はノルムを取ると $\sin ^2 \theta$が残ります。

このあままだと、下記になってしまうため何かが違うと思うのですが、
わからないのでPendingとしておこうとおもいます。。

$$\sigma_{\Delta} = \sqrt{\sigma_{dr}^2 + r^2 (\sigma_{\delta}^2 + \sin ^2 \theta \sigma_{\epsilon}^2)}$$

回答7.2.2

式7.12 から、問題の条件のときの誤差の標準偏差をもとめたいとおもいます。

import numpy as np

# パラメータ
r = 700  # 距離 (km)
dr = 0.01  # 距離誤差 (km)
d_deg = 0.01  # 仰角誤差 (度)
e_deg = 0.01  # 方位角誤差 (度)

# 度からラジアンへの変換
d_rad = np.deg2rad(d_deg)
e_rad = np.deg2rad(e_deg)

# 誤差の標準偏差
sigma_dr = dr
sigma_delta = d_rad
sigma_epsilon = e_rad

# 位置決定誤差の標準偏差の計算
sigma_delta_rp = np.sqrt(sigma_dr**2 + r**2 * (sigma_delta**2 + sigma_epsilon**2))

# 結果の表示
print(f"地上位置決定誤差の標準偏差: {sigma_delta_rp:.6f} km")

地上位置決定誤差の標準偏差: 0.173068 km　です

問題7.3

1. 式 (7.19) から標準偏差の式 (7.20) を求めなさい。
2. $\delta_1, \delta_2 = 1 \text{Hz}$、 $0.01 \text{km}$、 $\lambda = 23.6 \text{cm}$、 $r = 700 \text{km}$、 $x = 400 \text{km}$、 $v = 7 \text{km/s}$ としたときに、地上位置決定誤差を計算しなさい。

回答7.3.1

式 (7.19) から式 (7.20) を導くために、問題7.2 と同様に期待値と標準偏差の関係を利用します。　個別の誤差がガウス分布によれば足していくと0になるという性質を利用します。

式 (7.19)

$$\Delta r'_p = \begin{pmatrix} \left(1 - \frac{f_d \lambda}{2v}\right)^2 \left( - \frac{r^2 \lambda}{2xv} \delta_1 + \frac{r}{x} \delta_2 \right) \\ \frac{\lambda}{2v} \left( r \delta_1 + f_d \delta_2 \right) \end{pmatrix}$$

1. 各項の展開

まず、上記の式を展開し、各誤差成分について2乗の期待値を求めます。

2. 二乗の期待値を取る

式 (7.19) の各項の二乗を考えます。

上の成分（x成分）

$$\left(1 - \frac{f_d \lambda}{2v}\right)^2 \left( - \frac{r^2 \lambda}{2xv} \delta_1 + \frac{r}{x} \delta_2 \right)$$

これを2乗すると、

$$\left[\left(1 - \frac{f_d \lambda}{2v}\right)^2 \left( - \frac{r^2 \lambda}{2xv} \delta_1 + \frac{r}{x} \delta_2 \right)\right]^2$$

$$= \left(1 - \frac{f_d \lambda}{2v}\right)^4 \left( \left(- \frac{r^2 \lambda}{2xv} \delta_1\right)^2 + \left(\frac{r}{x} \delta_2\right)^2 + 2 \left(- \frac{r^2 \lambda}{2xv} \delta_1\right) \left(\frac{r}{x} \delta_2\right) \right)$$

ここで、クロス項の期待値 $\mathbb{E}[\delta_1 \delta_2]$ は独立性、ガウス分布を仮定しているのでゼロとなる

$$= \left(1 - \frac{f_d \lambda}{2v}\right)^4 \left( \left(- \frac{r^2 \lambda}{2xv} \delta_1\right)^2 + \left(\frac{r}{x} \delta_2\right)^2 \right)$$

下の成分（y成分）

$$\frac{\lambda}{2v} \left( r \delta_1 + f_d \delta_2 \right)$$

これを2乗すると、

$$\left[ \frac{\lambda}{2v} \left( r \delta_1 + f_d \delta_2 \right) \right]^2$$

$$= \left( \frac{\lambda}{2v} \right)^2 \left( (r \delta_1)^2 + (f_d \delta_2)^2 + 2 r f_d \delta_1 \delta_2 \right)$$

ここでも、クロス項の期待値 $\mathbb{E}[\delta_1 \delta_2]$ は独立性によりゼロとなると仮定すると、

$$= \left( \frac{\lambda}{2v} \right)^2 \left( r^2 \delta_1^2 + f_d^2 \delta_2^2 \right)$$

3. 標準偏差の計算

上の式をまとめて標準偏差の式にすると

$$\sigma_{\Delta}^2 = \left(1 - \frac{f_d \lambda}{2v}\right)^4 \left( \frac{r^4 \lambda^2}{4x^2v^2} \sigma_{\delta_1}^2 + \frac{r^2}{x^2} \sigma_{\delta_2}^2 \right) + \frac{\lambda^2}{4v^2} \left( r^2 \sigma_{\delta_1}^2 + f_d^2 \sigma_{\delta_2}^2 \right)$$

展開して

$$\sigma_{\Delta}^2 = \left(1 - \frac{f_d \lambda}{2v}\right)^4 \frac{r^4 \lambda^2}{4x^2v^2} \sigma_{\delta_1}^2 + \left(1 - \frac{f_d \lambda}{2v}\right)^4 \frac{r^2}{x^2} \sigma_{\delta_2}^2 + \frac{\lambda^2}{4v^2} r^2 \sigma_{\delta_1}^2 + \frac{\lambda^2}{4v^2} f_d^2 \sigma_{\delta_2}^2$$

項を整理すると

$$\sigma_{\Delta}^2 = \left(\left(1 - \frac{f_d \lambda}{2v}\right)^4 (\frac{r^2 \lambda}{2xv})^2 + (\frac{\lambda r}{2v} )^2 \right) \sigma_{\delta_1}^2 + \left(\left(1 - \frac{f_d \lambda}{2v}\right)^4 (\frac{r}{x})^2 + (\frac{\lambda f_d}{2v} )^2 \right)\sigma_{\delta_2}^2$$

よって、下記の7.20 が導けました。

$$\sigma_ = \sqrt{ \left[ \left( 1 - \frac{f_d \lambda}{2v} \right)^4 \left( \frac{r^2 \lambda}{2xv} \right)^2 + \left( \frac{\lambda r}{2v} \right)^2 \right] \delta_1^2 + \left[ \left( 1 - \frac{f_d \lambda}{2v} \right)^4 \left( \frac{r}{x} \right)^2 + \left( \frac{\lambda f_d}{2v} \right)^2 \right] \delta_2^2 }$$

回答7.3.2

上記(7.20) のfd は （本文式7.17）$f_d=\frac{2vy}{\lambda r}$より、
真横を計測している時 y=0 と考えて以下の式で近似できます。

$$\sigma_ = \sqrt{ \left[ \left( \left( \frac{\lambda r}{2v} \right)^2 \left( 1 + \frac{r}{x} \right)^2 +\right) \right] \delta_1^2 + \left[ \left( \frac{r}{x} \right)^2 \right] \delta_2^2 }$$

import numpy as np

# 条件の設定
delta1 = 1  # Hz
delta2 = 0.01  # km
lambda_ = 23.6e-5  # cm to km
r = 700  # km
x = 400  # km
v = 7  # km/s


# 各項の計算
term1 = (lambda_ * r / (2 * v))**2 * (1 + (r / x)**2) * delta1**2
term2 = (r / x)**2 * delta2**2

# 地上位置決定誤差の計算
sigma = np.sqrt(term1 + term2)

# 結果の表示
print(f"地上位置決定誤差: {sigma:.6f} km")

地上位置決定誤差: 0.029528 km
（本書の解説だと、誤差は0.0238km となっているのでどこかfdの想定が違うと思うのですが分かりませんでした）

課題7.4

リモートセンシングでは移動体にセンサを搭載し観測対象物の位置を計測する。今、図7.Q1のように、地表高さHとX方向に速度Vで飛翔している航空機があり、それにセンサを搭載してレーダを搭載し、進行前方向の地上に電波を照射した。電波は即座に反射し、観測機で受信し、あわせてドップラー周波数を計測した。このドップラー周波数を用いて観測対象物の位置を計測することを考える。なお、観測はXの平面内で行われるものとする。

課題7.4.1

(1) レーダの送信周波数をf0としたとき、Pの対象物から反射した信号のドップラー周波数はどのように表現されるか？なお，光速をcとする（ドップラー周波数とは通信に使用した信号の周波数が移動体の視線方向の速度 / 送信信号速度分だけ偏移したもの）。

前の6章(式6.7)[https://zenn.dev/kozayupapa/articles/0879fb2ea61e2b#本文-式6.7の解説] でドップラー周波数が下記になると習ったのでそれをそのまま利用します。
$f_D = \frac{2}{\lambda}|v|\cos\theta$

光の速度とマイクロ波の波長、周波数には下記の関係があります。
$c=\lambda f_0$

よって、ドップラー周波数は下記のように表すことができます
$f_D = \frac{2f_0}{c}V\cos\theta$
また本文の図より
$\cos\theta = \frac{x}{\sqrt{x^2 + H^2}}$

よって
$f_D = \frac{2f_0}{c}V\frac{x}{\sqrt{x^2 + H^2} }$

課題7.4.2

(2) 次に、xはドップラー周波数でどう表現できるか？

こちらは解説を参考にすると、近似で考えていくようです。
上記7.4.1 の式を展開して考えて、

$$\frac{f_D c}{2 f_0 V} \sqrt{x^2 + H^2} = x$$

まず、以下の変数を定義します：

$$a = \frac{f_D c}{2 f_0 V}$$

これにより、式は次のように書き換えられます：

$$a \sqrt{x^2 + H^2} = x$$

式の変形

1. 両辺を2乗する

両辺を2乗して、平方根を取り除きます：

$$a^2 (x^2 + H^2) = x^2$$

2. 方程式を展開する

方程式を展開して整理します：

$$a^2 x^2 + a^2 H^2 = x^2$$

3. xについて解く

すべての $x^2$ の項を一方にまとめます：

$$a^2 H^2 = x^2 - a^2 x^2$$

これをさらに整理すると：

$$a^2 H^2 = x^2 (1 - a^2)$$

したがって、 $x^2$ は次のようになります：

$$x^2 = \frac{a^2 H^2}{1 - a^2}$$

したがって、 $x$ は次のように表されます：

$$x = \frac{a H}{\sqrt{1 - a^2}}$$

このあとですが、　結局aが十分に小さい数なのであれば
下記の近似がなりたちます。

$$x \approx a H$$

最終的には下記になります

$$x \approx \frac{f_D c H}{2 f_0 V}$$

課題7.4.3

(3) これでドップラー周波数を計測して，位置を見出すことが可能になった。計測には誤差がつきものである。今，ドップラー周波数の計測誤差を$\delta f_Dとすると，位置決定精度はどのように表現されるか？

上記の式の$f_D$が $f_D + \delta f_D$となるだけなので　誤差は単純に代入して下記になります。

$$\delta x \approx \frac{\delta f_D c H}{2 f_0 V}$$

課題7.4.4

(4) 最後に，$f_0$を1.27GHz，V=7 km/s,H=500Km, $\delta f_d$=1Hzとしたとき、位置決定精度はいくらくらいになるか？

import numpy as np

# 条件の設定
f0 = 1.27e9
delta = 1  # Hz
H = 500  # km
v = 7  # km/s
c = 300000

sigma = (delta * c * H )/(2*f0*v)
# 結果の表示
print(f"位置決定誤差: {sigma:.6f} km")

位置決定誤差: 0.008436 km

---

以上です。　なんとか、1週間かけて　第七章　解いていきました。
いやー、　SARは難しいですね。。