はじめに
スペックル(Speckle)をPython使って計算してみるシリーズの第2回です。実際の物質に近い表面状態に対してスペックル計算します。PSD(Power Spectral Density)を導入して、実際の高さ場を求めて位相スクリーン近似でスペックルを求めてみます。一般にスペックルは「粗面」で発生すると言われますが、ミラーのような平坦性の高い条件でもスペックルが発生します。
前回記事
スペックル計算の方針
実際の試料の高さを計算したあと、位相スクリーン近似(高さのばらつきが付加する位相に対応するという近似)によりスペックルを計算したいと思います。
試料は平坦度が高いことを前提とします。高さ場はAFMによる実測データがあるといいですが、実際はそういきません。そこでPSDの情報から高さ場を計算したいと思います。
PSDとは
表面の高さ場（プロファイル）を $h(x, y)$ とするとき、そのフーリエ変換の絶対値の2乗を、測定領域の面積で規格化したものがPSDです。
2次元の空間周波数を $(f_x, f_y)$、測定領域のサイズを $L_x \times L_y$ とすると、定義式は以下のようになります。
$$PSD(f_x, f_y) = \lim_{L_x, L_y \to \infty} \frac{1}{L_x L_y} \left| \int_{-L_y/2}^{L_y/2} \int_{-L_x/2}^{L_x/2} h(x, y) e^{-i 2\pi (f_x x + f_y y)} dx dy \right|^2$$
やや感覚的ですが、PSDのルートをとれば高さ場のフーリエ係数が求められるので、逆フーリエ変換で高さ場が求まります。ただし、PSDに位相情報はないのでこちらで用意します。
特定の光学応答を示す粗面や高さ場は無数にありますが、特定のPSDに対して位相の取り方を変えることでこれらの無数の粗面や高さ場を表現できるのです。
試料の高さ場の求め方
PSDはガウシアンと仮定して、相関長さは5μmとして計算してみます。高さ場の位相はランダムとして求めます。
import numpy as np
import matplotlib.pyplot as plt
#PSD関数の定義
def Gaussian_PSD(KX, KY, corr=5e-6):
return np.exp(-(KX**2 + KY**2)*corr**2)
N = 1024
L = 100e-6
wavelength = 532e-9 #波長532nm
x = np.linspace(-L/2, L/2, N)
y = x.copy()
dx = x[1] - x[0]
X, Y = np.meshgrid(x, y)
# --- ifft 用の周波数軸 ---
kx = np.fft.fftfreq(N, d=dx) # 長さ N
ky = np.fft.fftfreq(N, d=dx) # 長さ N
KX, KY = np.meshgrid(kx, ky)
# ランダム位相
phase = 2 * np.pi * np.random.rand(N, N)
PSD = Gaussian_PSD(KX,KY)
F = np.sqrt(PSD) * np.exp(1j * phase)
# 高さ場
h = np.fft.ifft2(F, s=(N, N))高さ場をhとして求めることが出来ました。一点、値を確認してみます。
h[512,512]np.complex128(4.195453653622023e-05+1.785581902037705e-05j)
値が複素数になっています。高さ場は明らかに実数なので、これはおかしな結果です。
エルミート対称性を用いた高さ場の実数化
高さ場を求める際に位相はランダムとして乱数を引きましたが、本来は逆フーリエ変換した後に実数となるように制約を加える必要があるのです。
高さ場のフーリエ変換をFとして高さ場hは下記と書けます。
$$h(x, y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(k_x, k_y) e^{i 2\pi (k_x x + k_y y)} dk_x dk_y$$
この時、Fに対して以下の制約を付けると変換後が実数になります。イメージとしては正の周波数と負の周波数がペアにして打ち消す感じです。*は複素共役の意味です。
$$F(-k_x, -k_y) = F^*(k_x, k_y)$$
この特性をエルミート対称性(Hermite symmetry)と呼びます。教科書だといきなりエルミート共役の定義が出てきて???となりますが、自力で計算していると自然と必要な概念だと分かります。
実際に計算する時は、エルミート共役自体をコード化するよりnumpyのrfft2を使う方が安定するようです。
import numpy as np
import matplotlib.pyplot as plt
# 1. パラメータ設定
N = 1024
L = 100e-6
dx = L / N
# 周波数空間のPSDを (N, N//2 + 1) の形状で定義
# kx_r は np.fft.rfftfreq(N, d=dx) を使用
kx_r = np.fft.rfftfreq(N, d=dx)
ky_r = np.fft.fftfreq(N, d=dx)
KXr, KYr = np.meshgrid(kx_r, ky_r) # shape: (N, N//2+1)
# 2. PSDと複素振幅の生成
def Gaussian_PSD(KX, KY, corr=5e-6):
return np.exp(-(KX**2 + KY**2)*corr**2)
PSD_r = Gaussian_PSD(KXr, KYr)
PSD_r[0, 0] = 0 # DC除去
# 3. 位相を振る
phase = 2 * np.pi * np.random.rand(N, N // 2 + 1)
F = np.sqrt(PSD_r) * np.exp(1j * phase)
# 4.irfft2で逆変換
h = np.fft.irfft2(F, s=(N, N))
# 5. 正規化
h = (h - h.mean()) / h.std() * 1e-9最後に高さ場の標準偏差が1nmとなるように最後にスケールを調整しました。
高さ場を表示させると下図です。高さ1nmスケールといいましたが、最大で4nm程度までばらつきがありますね。
ガウシアンPSDの相関長さを短くすれば、高周波成分が許容されるので山谷を多くしたり出来そうです。
スペックル計算方法
前節で求めた高さ場から、付加される位相を求めます。位相スクリーン近似を使います。
位相スクリーン近似
位相スクリーン近似（Phase Screen Approximation）」は、複雑な形状の物体や乱れた媒質を通過する光の挙動を、「薄い膜（スクリーン）を通過した際に、光の位相が局所的に変化する」というモデルに置き換えて計算する手法です。
表面の高さ分布に従って「波面が歪められる」と考えて計算します。付加される位相を$\phi$とすると下式であらわせます。
$$\phi(x, y) = \frac{4\pi}{\lambda} h(x, y) \cos \theta$$
$\theta$は入射角で、今は$\cos \theta$は1とします(垂直照明)。k<1でマスクをつけます。
# 位相スクリーン近似
phi = 4* np.pi / wavelength * h
高さばらつきでどの程度位相が付加されるか”phi”を図示してみます。
カラーバーの値を見ると、最大でも2$\pi$の1.5%程度の位相が付加されてるだけです。位相の大きさとしてはとても小さいと言えます。
続きです。
# 試料面の複素電場(電場強度一定)
field = np.exp(1j * phi)
# フーリエ変換（遠方回折 / フラウンフォーファー回折）
# FFTのシフトを行い、強度（絶対値の2乗）を計算
fft_field = np.fft.fft2(field)
speckle = np.abs(fft_field)**2
speckle_shifted = np.fft.fftshift(speckle)
# スペックルコントラストの計算 (標準偏差 / 平均値)
contrast = np.std(speckle) / np.mean(speckle)
#波数を計算
k = 1 / wavelength#532nmの波数[cycle/m]
kx_s = np.fft.fftshift(kx)
kx_s /= k
ky_s = np.copy(kx_s)
KX_s,KY_s = np.meshgrid(kx_s,ky_s)
#k<1でマスク
mask = KX_s**2 +KY_s**2 <1
fig, ax = plt.subplots(figsize=(6,6))
c = ax.contourf(KX_s,KY_s,speckle_shifted*mask ,levels=50, cmap="gray")
ax.set_aspect("equal")
ax.set_xlabel("KX")
ax.set_ylabel("KY")
fig.colorbar(c, ax=ax, label="Intensity",shrink=0.7)
plt.title(f"Speckle Pattern\nContrast: {contrast:.2f}", fontsize=14)
plt.tight_layout()
plt.show()
真ん中に輝点があるだけになってしまいました。(見えない人はズームしてみてください。)
さらにコントラストが100近くになってしまいました。どういうことでしょうか？
今の条件だと試料を垂直照明しているため、正反射光が入っているのです。正反射光はとても明るいため、そこだけが見えているのです。更にコントラストの値も輝点の影響を受けて大きな値になっているのです。
ちょっと数字を調整して、対数表示してみます。
#定数を足して、対数を取る
c = ax.contourf(KX_s,KY_s,np.log(speckle_shifted*mask+1000) ,levels=50, cmap="gray")
明るさは均一になっておらず、今回のように試料表面の1nm程度の高さばらつきでもぼつぼつとし干渉縞がみえます。
考察
マスクを作って正反射光をブロックした上で線形スケールで画像作成しましたが、その場合境界付近のみが明るくなってしまいます。つまり正反射光はここまでだと切り分けること自体が難しそうでした。
今回はPSDの相関長さを5μmとしましたが、これを変えて様子をみていくのも面白いかもしれません。
前回計算では全く見えなかった正反射光が強く出たのは、付加した位相量が2$\pi$の1.5%程度と小さかったためと思います。位相自体は小さいですがランダムなためスペックルは現れますが、正反射光は元の性質を維持するため輝度が高くなったものと思います。
位相スクリーン近似の式を見ると、高さばらつきが波長の半分程度であれば、付加される位相は2$\pi$に達して、完全発達スペックルになり正反射は気にならなくなると考えられます。
表面がある程度滑らかな物質では、正反射光の影響が強いためスペックルコントラストの計算には測定領域を制限する等の工夫が必要となりそうです。
まとめ
PSDの値から試料表面の高さ場を求めて、位相スクリーン近似でスペックルを計算しました。
・PSDは一般の粗面を表すのに適した物理量である
・高さ場のような実数量をフーリエ変換(逆フーリエ変換)から求めるときは、エルミート共役化する必要がある
・高さ1nm程度のばらつきでもスペックルは発生するが、正反射光が強く見える
例えばSiの格子定数が約0.5nmということを考えると、この結果はスペックルは一般的に起こっている現象だと言えそうです。
次回記事
今回は実際の形状を取り入れましたが、次はこれに加えて、垂直照明でなく斜方照明条件として更に屈折率の影響も取り入れて計算できればと思います。
参考文献
論文の引用元や海外コードの解説で必ず登場するグッドマンの不朽の名著。英語の専門用語や理論的なニュアンスをそのまま正確に理解したい、本格志向の技術者におすすめの一冊です。
散乱光学について体系的に学びたい方には、John C. Stover 著 Optical Scattering: Measurement and Analysis がおすすめです。実際の測定を意識した解説が多く、表面粗さや散乱理論、BRDFなど、研究や開発で役立つ内容が幅広くまとまっています。
PSDについても説明があり数式だけでなく物理的な解釈にも触れられているため、散乱現象をより深く理解したい方に適した一冊です。
私自身も散乱光のシミュレーションや理論を調べる際に、何度も読み返している参考書です。
回折計算（FFT）をコードに落とし込む際の物理的・数学的な根拠を最も丁寧に解説してくれる不朽の名著です。光学シミュレーションで迷ったときの道標として、手元に置きたい一冊。
コメント