應用特高頻陣列雷達三維成像技術於重建大氣亂流之空間結構

應用特高頻陣列雷達三維成像技術於重建大氣亂流之空間結構

陳振雄

摘要 — 本研究發展多載頻和多接收機之雷達干涉法， 並應用於重建大氣亂流層狀結構之三維空間分布，分析 的資料來自日本 MU 特高頻陣列雷達之觀測。觀測時發 射 5 個頻率並使用 19 個接收通道之資料，再以 Capon 非線性反演方法重建雷達照射體積內之回波強度空間分 布，藉此呈現主要亂流結構之小尺度變化，有助於進行 大氣微觀尺度之動力分析。1 一、簡介 結合多接收機和多載頻之雷達觀測技術可以進行目標 物之空間定位，本技術已應用於脈波式的特高頻(VHF) 陣列雷達，藉以執行大氣折射指數不規則體結構(通稱 亂流)之小尺度三維成像和定位分析。 本文呈現實際觀測的分析結果，資料來源為日本 MU 特高頻陣列雷達(由京都大學所建立與管理)，資料時間 為 2006 年 2 月 9 日 20:16-20:25UT，每一點資料取樣時 間約為 0.512 秒。雷達陣列組態如圖一所示，觀測時以 整 個 天 線 場 送 出 雷 達 波 ， 垂 直 發 射 5 個 頻 率 (46.0、46.25、46.5、46.75、47.0MHz)，並啟用 21 個接 收通道，包含整個天線陣列場接收和中間部份的 20 個 接 收 天 線 組 ( 編 號 A2、A3、A4、B2、B3、B4、C2、C3、C4、D2、D3、 D4、E2、E3、E4、F1、F2、F3、F4、F5)，其中 F1 通 道之資料未使用。脈波長 1 s，因此取樣資料之徑向解 析度約相當於 150m，雷達波束垂直向上，波束寬約為 3.7o _{(-3dB 的寬度)，因此雷達脈波照射的空間尺度為數} 百公尺。成像分析則使用 Capon 演算法[1]，此演算法 可有效處理多維度訊號且較為省時，但是隨著訊雜比值 的降低，解析效果會減弱。 圖一：日本 MU VHF 雷達組態，A1-F5 為用於接收的天線群組。 二、雷達三維成像技術 設定 M 個載頻 f1, f2,…, fM，依序隨雷達脈波輪流發射， 地面上則啟用 N 個接收機，座標定為 D1, D2,…,Dn。各 個 接 收 機 所 接 收 到 的 各 載 頻 回 波 寫 成 一 個 欄 向 量 (column vector)： T MN M M N N t S t S t S t S t S t S t S t S t S (t) )] ( ) ( ) ( ) ( ... ) ( ) ( ) ( ) ( ) ( [ 2 1 2 22 21 1 12 11     S , (1) 其中上註標 T 為矩陣之轉置(transpose)，Sij(t)為接收機 j 1 _{本研究由國科會贊助，計畫編號 NSC 101-2111-M-270-001 及 NSC} 102-2111-M-270-001。 收到的載頻 i 之回波訊號。理想上，M 個載波須同時發 射，但實際操作上只能於每個雷達脈波間切換。由於亂 流之變化時間約為 1 秒鐘，而雷達脈波發射時間差僅為 數百個微秒(s)，在此極短時間內亂流並未產生明顯變 化，因此各載頻之回波可視為由相同之亂流所散射或反 射回來。 定義 S(t)之互相關函數矩陣如下：





























MM

N

M

M

M

H

V

V

V

V

V

V

t

t

V

V

V

)

(

)(

2

1

2

22

21

1

12

11











S

S

V

, (2) 其中上註標H為Hermitian 運算子，V是 MN×MN之矩陣 (又稱為visibility matrix)，Vpq為N×N 矩陣，為各接收機 訊號之互相關函數矩陣，p和q代表載頻編號，例如第1 個和第2個載頻之互相關函數矩陣為































*

2

1

*

22

1

*

21

1

*

2

12

*

22

12

*

21

12

*

2

11

*

22

11

*

21

11

12

N

N

N

N

N

N

S

S

S

S

S

S

S

S

S

S

S

S

S

S

S

S

S

S

V















, (3) 符號*為複數之共軛運算(conjugate)。使用Capon演算法 之三維成像運算式為 e V e a a ₂ 1 ₁ ) ( 1 ) (  _{H } , R W , R B , (4) 其中

B(a, R): 亮度分布(brightness distribution)

a=[sinθsinφ, sinθcosφ, cosθ]

θ 和 φ: 偵測方向之天頂角和方位角 R: 相對於雷達體積中心點之距離變數

W(a, R): 三維之權重函數(weighting function)，由雷

達波束和距離權重函數所形成。

變數 e 為一個權重向量(weighting vector)，用以反演亮 度分布值，其形式為

T k R k j k R k j k R k j k R k j k R k j k R k j k R k j k R k j k R k j N M M M M M M N N e e e e e e e e e ] [ ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 2 1 2 2 2 2 2 1 2 2 1 1 2 1 1 1 1 1 D a D a D a D a D a D a D a D a D a e                                    , (5) 其中 ki是載波 i 之波數(wavenumber)。 方程式(4)顯示反演之亮度分布是受到權重函數 W(a, R)影響後之結果，因此 W(a, R)之效應需適當移除，才 能取得較為合理之亮度分布，用以呈現亂流之三維結構 和位置[2]。本研究參考[3]針對日本 MU 雷達所提方法 與方程式，進行 W(a, R)效應之移除，其詳細過程不再 贅述，本文僅呈現修正後之成像結果。 三、觀測分析結果 3.1 多載率觀測分析之大氣亂流層狀結構 在三維成像技術中，若只利用一個接收機的多載頻資 料進行分析，則只可得到徑向距離的亂流結構隨時間之 分布和變化，稱為距離成像法(range imaging/RIM) [4]， 其簡介與方程式參考附錄。本研究中 RIM 之資料來自 於以整個天線陣列場為接收通道所接收之回波訊號，此 時雷達的波束最窄，可獲得較佳的徑向結構。 圖二是本研究所分析的一個案例，縱軸為取樣層高度， 橫軸為時間軸，時間由右邊開始。圖中有四個取樣層， 每一計算值之時間解析約為 16 秒。上圖為原始的回波 功率圖，距離解析度為 150 m，所以無法呈現各層內部 的亂流或層狀結構。下圖是距離成像分析後的結果，每 2 m 距離即計算一個回波強度值(亮度)，大大提升了距 離解析度。其中，在 5.25-5.40km 之取樣層中可以看到 一個受到大氣波動所調制的亂流層狀結構，波峰並略微 向右傾斜，此傾斜現象是由於水平風速在垂直方向上的 逐漸增大所致，即風切效應，可參考圖四之水平風場示 意圖。實際分析之水平風場即類似於此(未呈現於本文 中)。 3.2 大氣亂流層之三維結構 分析圖二中受到大氣波動所調制的亂流層狀結構，圖 三是三個連續時刻 T4、T3、T2的三維結構解析圖，即圖 二中編號第 4、3、2 的位置點。觀察圖二之波動性層狀 結構，可判斷第 3 位置點剛好是波峰通過雷達上空之時 間。 為清楚顯示亂流結構在雷達波束照射體積內之變化， 圖三只顯示九層橫切面之輪廓線圖，呈現各高度層回波 功率之二維水平角分布。輪廓線之中心代表較強的回波 有時可見到二個回波中心，例如中間 T3的時刻。參照 圖二，此時刻之亂流層在 5.325 km 下方，即位於取樣 層之中心高度下方。第 3 位置點是波峰處，在其中心高 度下方之橫切面輪廓線可看到有二個明顯的回波中心。 本現象可由圖四之波動層狀結構模式說明之。在第 T3 時刻，雷達波束照射到波峰處，但由於層狀結構在波峰 二邊反方向傾斜，雷達接收到的訊號主要來自波峰二邊 的散射或反射回波。 相較於 T3時刻，T4時刻的傾斜層狀結構使得雷達回 波主要來自雷達正上方之左邊，T2時刻的傾斜層狀結構 則使得雷達回波主要來自雷達正上方之右邊。圖三的橫 切面輪廓線圖清楚顯示這兩種回波位置的特性。各個時 刻的分析結果說明多載頻和多接收機三維成像之可行性。 另外，圖三中的 range-meridional 平面上顯示二條剖 線，黑色線是橫切面成像回波各自加起來後的高度剖線 紅色線則是只取天頂方向的成像回波之高度變化，皆可 顯示層狀結構在徑向上的分布。但前者包含偏離天頂角 的回波，所以若偏離天頂角較遠的地方有強回波，通常 會加大層狀結構在徑向上的厚度。不過，在本案例中二 者差異不大。但若只取天頂方向的成像回波，可以有效 降低非亂流回波之污染，例如飛機反射波或其他訊號源 干擾，這些非亂流回波通常來自較大的天頂角方向[5] [6]。 回 波 功 率 (2006/02/09) (dB) 距 離 成 像 (2006/02/09) (dB) 圖二：多載頻技術反演之層狀結構，圖中 30、25 等數字表示不同時刻 之編號，時間從右邊開始。 30 25 20 15 10 5 T4 T3 T2

圖三：三維成像之結果，呈現九個橫切面之輪廓線，由左至右為圖二中第 4、3、2 時刻。粗虛線顯示各橫切面回波中心的連接。在 Range-meridional 平面上的兩條剖線是亮度垂直分布，紅色線為天頂角方向的亮度垂直分布，黑色線是橫切面各自加總後的亮度垂直分布。 水平風場 T4 T3 T2 圖 四 ： 層 狀 結 構 模 式 與 雷 達 回 波 示 意 圖 ， 雷 達 波 束 垂 直 向 上 。 T2、T3、T4 表示雷達波束照射的時間點。雷達波束內向下箭號表示回 波來向。 結論 多載頻和多接收機技術已使用於特高頻陣列雷達，可 以重建大氣亂流層狀結構之三維空間分布。本研究取用 日本 MU 特高頻陣列雷達的觀測資料，針對一個波動性 的亂流層狀結構，不僅利用多載頻成像法呈現出層狀結 構在雷達體積內的高度和時間變化，也解析出層狀結構 在雷達體積內的回波來向角變化，說明三維成像技術的 適用性。由於應用於大氣研究的特高頻陣列雷達建構經 費較高，所以無法像傳統的氣象雷達一樣普遍設立。不 過，隨著陣列雷達觀測技術不斷更新和精進，例如脈波 相位加碼、相控快速切換雷達波束方向、取樣和訊號處 理速度的提升等，使得特高頻陣列雷達探測大氣的能力 不斷提升，資料分析之解析度也不斷提高，這對大氣微 觀物理現象的觀測和研究會有很大的助益。 附錄 多載頻之距離成像法須設定 M 個載頻 f1, f2,…, fM，依 序隨雷達脈波輪流發射與接收。利用一些非線性的反演 方法，可由這些載頻回波訊號來重建亂流之徑向結構。 適用於目前雷達和觀測資料的非線性訊號反演方法為 Capon 分析法，並以在時間域之計算較為方便，定義 B(r)為反演後在距離 r 的亮度分布，其值代表 r 處的回 波強度。Capon 分析法如下：

e

R

e

1

1

)

(

r



_{H }

B

(6) 其 中 e = [ ej2k1r _ej2k2r _{… e}j2knr _]T_{， 為 為 一 個 權 重 向 量} (weighting vector)，用以反演亮度分布值，ki是載率 i 之 波數，T 表示矩陣轉置動作。(6)式中的註標 H 和-1 分別 代表 Hermitian 運算子(conjugate 加 transpose)和反矩陣 運算。R 是二載頻回波訊號間沒有時間延遲的互相關函 數之矩陣，形式為





































nn

n

n

n

n

V

V

V

V

V

V

V

V

V

...

.

.

.

.

.

.

...

...

)0

(

2

1

2

22

21

1

12

11

R

(7)









₀

_,

₎

*

(

_i

_j

ij

k

V

V

V



(8) Vij為互相關分析(cross-correlation analysis)，並為一對載 頻波數差 Δk 之函數，Vi和 Vj是第 i 個和第 j 個載率所收 之時間序列訊號。 對應於(6)式，另有 Fourier 分析法，其亮度計算為 B(r) = e+_{R e} (9) 相較於 Capon 分析法，Fourier 法的解析效果較差。 圖五是模擬的案例，呈現 Fourier 和 Capon 分析法的 解析度差異。模擬中所使用之雷達回波方程式和相關運 算式可參考[7]之附錄，本文略過其詳細描述。在本次模 擬案例中設定雷達脈波寬為 2μs，即取樣資料之徑向解 析 度 約 相 當 於 300m ， 並 使 用 六 個 載 率 52.75, 53.05, 53.35, 53.65, 53.95, 54.25 MHz。圖五(a)之模型中設定二 個層狀結構，在 300 m 的取樣層中厚度和位置都有變化， 圖五(b)和(c)分別是 Fourier 和 Capon 分析法之結果，相 較之下，Capon 分析法確實有較佳的效果。Fourier 分析 法雖然可以分辨出這二層結構，但無法解析厚度之變化 更多模擬案例顯示，當層狀結構更靠近時，Fourier 分析 法無法分離這二層結構。 圖 五 ： 多 載 頻 距 離 成 像 法 模 擬 分 析 ， 脈 波 寬 2μs ， 六 個 載 率 52.75、53.05、53.35、53.65、53.95、54.25 MHz。(a) 二層狀模型， (b) Fourier 分析法，(c) Capon 分析法。 誌謝 本研究之雷達資料由日本 MU 雷達站所提供，該雷達 站由日本京都大學生存圈研究所(RISH)所經營與管理。 感謝 RISH 之 G. Hassenpflug 博士、T. Tsuda 和 M. Yamamoto 教授之協助觀測。

參考文獻

[1] Palmer, R. D., S. Gopalam, T.-Y. Yu, and S. Fukao, “Coherent radar imaging using Capon’s method,” Radio Sci., vol. 33, 1585-1598, doi:10.1029/98RS02200, 1998.

[2] Yu, T.-Y. and R. D. Palmer, “Atmospheric radar imaging using multiple-receiver and multiple-frequency techniques,” Radio Sci., vol. 36, 1493-1503, doi:10.1029/2000RS002622, 2001

[3] Chen, J.-S. and J. Furumoto, “A novel approach to mitigation of radar beam weighting effect on coherent radar imaging using VHF atmospheric radar,” IEEE Trans. Geosci. Remote Sens., to be published, 2011.

[4] Palmer, R. D., T.-Y. Yu, and P. B. Chilson, “Range imaging using frequency diversity,” Radio Sci., vol. 34, 1485–1496, 1999.

[5] Hassenpflug, G., M. Yamamoto, H. Luce, and S. Fukao, “Description and demonstration of the new Middle and Upper atmosphere Radar imaging system: 1-D, 2-D and 3-D imaging of troposphere and

stratosphere,” Radio Sci., vol. 43, RS2013,

doi:10.1029/2006RS003603, 2008.

[6] Yu, T.-Y., J. Furumoto, and M. Yamamoto, “Clutter suppression for high-resolution atmospheric observations using multiple receivers

and multiple frequencies,” Radio Sci., 45, RS4011,

doi:10.1029/2009RS004330, 2010.

[7] Chen, J.-S., J. Furumoto, and T. Nakamura, “Effects of radar beam width and scatterer anisotropy on multiple-frequency range imaging using VHF atmospheric radar,” Radio Sci., 45, RS5001, doi:10.1029/2009RS004267, 2010.

(b) Fourier