提取準(zhǔn)雙周和準(zhǔn)一月低頻振蕩的Lanczos濾波器及其應(yīng)用

姚菊香，李麗平，羅璇，楊瑋，王盤興

(南京信息工程大學(xué)大氣科學(xué)學(xué)院，江蘇南京210044)

0 引言

自從Madden and Julian(1971，1972)發(fā)現(xiàn)熱帶大氣的緯向風(fēng)和氣壓場(chǎng)存在40～50 d周期的低頻振蕩(LFO)以來，大氣中的低頻振蕩現(xiàn)象一直受到氣象學(xué)家高度重視。對(duì)于低頻振蕩的重要性，國(guó)內(nèi)外學(xué)者做過較為系統(tǒng)總結(jié)(李崇銀，2004;Zhang，2005)。青藏高原熱力作用對(duì)中高緯度地區(qū)低頻振蕩的產(chǎn)生具有重要作用(章基嘉等，1991;劉式適等，2000)。通過1979年第一次青藏高原氣象科學(xué)試驗(yàn)，章基嘉等(1984)發(fā)現(xiàn)高原熱源分量存在中周期(4～10 d)、短周期(2～3 d)的振蕩;周兵等(2000)指出高原東部地面感熱通量存在30～50 d低頻振蕩和準(zhǔn)雙周振蕩。

徐國(guó)強(qiáng)和朱乾根(2000)在分析1998年高原大氣低頻振蕩結(jié)構(gòu)特征時(shí)指出，高原存在多尺度的大氣低頻振蕩，高原和長(zhǎng)江中下游的低頻降水存在同時(shí)加強(qiáng)和同時(shí)減弱的關(guān)系。鞏遠(yuǎn)發(fā)等(2007)指出，2001和2003年夏季的亞洲季風(fēng)區(qū)，夏季5—8月期間，高原中南側(cè)較強(qiáng)的低頻熱源(熱匯)對(duì)后期江淮流域降水偏多(少)起重要作用。王躍男等(2009)指出1978(1999)年夏季青藏高原東部大氣熱源存在以10～20 d(30～60 d)周期為主的振蕩，并與長(zhǎng)江中下游降水在10～20 d(30～60 d)頻段存在顯著相關(guān)。

上述研究發(fā)現(xiàn)，青藏高原熱源存在多時(shí)間尺度特征，并對(duì)東亞大氣環(huán)流、中國(guó)氣候變化以及災(zāi)害性天氣的形成和發(fā)展有著重要的影響。為了更好地研究高原季內(nèi)振蕩的多尺度特征，一個(gè)優(yōu)良的濾波器是必要的。姚菊香等(2005)以30～60 d低頻振蕩為對(duì)象，定量分析了常用的Butter worth帶通濾波器(Murakami，1984;簡(jiǎn)記為B.f.)和Lanczos帶通濾波器(Duchon，1979;簡(jiǎn)記為L(zhǎng).f.)性能的優(yōu)劣，結(jié)果表明，對(duì)于取樣間隔為1 d的序列，當(dāng)L.f.濾波器參數(shù)l≥69時(shí)，其濾波性能全面優(yōu)于B.f.。本文的目的是根據(jù)青藏高原大氣熱源季內(nèi)振蕩分析的實(shí)際需要，討論兩個(gè)波段(10～20 d、20～40 d)L.f.的濾波特性，確定L.f.濾波器臨界參數(shù)，并以拉薩附近格點(diǎn)1950—2006年逐年整層熱源日序列為例，驗(yàn)證兩個(gè)波段L.f.功能，為高原及鄰近區(qū)域熱源季內(nèi)振蕩及傳播特征研究提供基本方法。

1 Lanczos濾波器

Lanczos濾波器(L.f.)是一種數(shù)字濾波器(Duchon，1979)。對(duì)采樣間隔Δt=1 d的時(shí)間序列

要濾出頻段［f1，f2］的分量，需要先計(jì)算L.f.的權(quán)重函數(shù)Lw(τ):

它由2l+1個(gè)數(shù)構(gòu)成(l為L(zhǎng).f.的參數(shù))。通過

可從序列(1)中濾出頻段［f1，f2］的分量;濾波過程使序列(1)縮短2l(單位:d)。(3)式中“*”為卷積算符。

L.f.的頻率響應(yīng)函數(shù)為

理想濾波器的響應(yīng)函數(shù)記為R(f)(燃料化學(xué)工業(yè)部石油地球物理勘探局計(jì)算中心站，1974)。

圖1和圖2分別給出了提取10～20 d與20～40 d兩個(gè)波段的L.f.的響應(yīng)函數(shù)LR(f)與R(f)的

圖1 提取10～20 d振蕩的L.f.的響應(yīng)函數(shù)LR(f)與理想濾波器響應(yīng)函數(shù)R(f)的比較(曲線、粗實(shí)水平線分別為L(zhǎng)R(f)、R(f))a.l=20;b.l=50;c.l=121Fig.1 Comparison of the response functionLR(f)of L.f.for filtering 10—20 d oscillation and the response function R(f)of ideal filter(The curve and thick solid horizontal line indicateLR(f)and R(f)filters，respectively)a.l=20;b.l=50;c.l=121

圖2 提取20～40 d振蕩的L.f.的響應(yīng)函數(shù)LR(f)與理想濾波器響應(yīng)函數(shù)R(f)的比較(曲線、粗實(shí)水平線分別為L(zhǎng)R(f)、R(f))a.l=20;b.l=50;c.l=121Fig.2 Comparison of the response functionLR(f)of L.f.for filtering 20—40 d oscillation and the response function R(f)of ideal filter(The curve and thick solid horizontal line indicateLR(f)and R(f)filters，respectively)a.l=20;b.l=50;c.l=121

比較?？梢?，隨著l的增大，LR(f)趨向于R(f)。其中:

2 L.f.濾波效果的定量分析

2.1 L.f.、B.f.和H.f.三種濾波器響應(yīng)函數(shù)的比較

判斷L.f.濾波效果的優(yōu)劣，可以通過比較它與Butter worth濾波器(B.f.)、諧波濾波器(H.f.)(馮康，1978)的響應(yīng)函數(shù)進(jìn)行。B.f.與H.f.都是分析季內(nèi)振蕩常用的濾波器，它們的頻率響應(yīng)函數(shù)分別記為BR(f)、HR(f)。

B.f.是頻率連續(xù)型濾波器。從序列(1)中濾出頻段［f1，f2］的B.f.的頻率響應(yīng)函數(shù)為

其中，Z、a、b1、b2取決于Δt、f、f1、f2的值的大小。通過

從序列(1)中濾出頻段［f1，f2］的分量。濾波不改變序列的長(zhǎng)度。

H.f.是頻域(或周期域)上的離散型濾波器(馮康，1978)。對(duì)(1)式的序列，有諧波分解式

式中:ak、bk是波數(shù)k的余弦、正弦諧波波數(shù);Δf=

1/n(d－1)是基頻。從(1)式中濾出的H.f.的頻率響應(yīng)函數(shù)為

濾波結(jié)果為

與B.f.類似，H.f.濾波也不改變序列的長(zhǎng)度。

圖3給出了提取夏季(n=92)中10～20 d、20～40 d波段(分別記為a、b波段)的B.f.、H.f.的頻率響應(yīng)函數(shù)BR(f)、HR(fk)分別與R(f)的比較。

HR(fk)在fk=k/n上與R(f)相等;BR(f)與R(f)、HR(fk)差異明顯。除中心頻率f0=上BR(f)達(dá)到R(f)外，其余頻率上BR(f)均不等于R(f)。

與圖1、圖2的LR(f)比較，對(duì)較小l(圖1a、2a)，兩波段上的BR(f)較LR(f)更接近R(f)，故B.f.優(yōu)于L.f.;而對(duì)較大l(圖1c、2c)，LR(f)較BR(f)更接近R(f)，故在l取值較大時(shí)，L.f.優(yōu)于B.f.。根據(jù)表1可知，a，b波段近似可以由的H.f.獲得。

2.2 L.f.濾波性能優(yōu)于B.f.的臨界參數(shù)l0的確定

圖3 BR(f)、HR(fk)與R(f)的比較(曲線為BR(f)，○為HR(fk)，粗實(shí)水平線為R(f))a.提取周期為10～20 d;b.提取周期為20～40 dFig.3 Comparison ofBR(f)，HR(fk)and R(f)(The curve indicatesBR(f)，○indicatesHR(fk)，and the thick solid horizontal line indicates R(f))a.the period of 10—20 d;b.the period of 20—40 d

表1 序列(1)的k=0，10的周期Tk和頻率fk(n=92)Table 1 Period Tkand frequency fkof time series(1)for k=n=92)

表1 序列(1)的k=0，10的周期Tk和頻率fk(n=92)Table 1 Period Tkand frequency fkof time series(1)for k=n=92)

k 012345678910 Tk/d∞92.046.030.723.018.415.313.111.510.29.2 fk/d－100.0110.0220.0330.0430.0540.0650.0760.0870.0980.109

率區(qū)間［0 d－1，0.5 d－1］分為3段:通過帶［f1，f2］、低頻抑制帶［0，f1)、高頻抑制帶(f2，0.5 d－1］，分別以i=標(biāo)記;則L.f.與B.f.在通過帶、低頻抑制帶、高頻抑制帶的輸出輸入方差比分別為

和

由此，對(duì)提取頻段［f1，f2］的濾波器，可定義L.f.與B.f.在通過帶和兩個(gè)抑制帶上的濾波效率比(以下簡(jiǎn)稱效率比)為

它是確定L.f.濾波性能超過B.f.的l臨界值的統(tǒng)計(jì)量。

表2和表3給出了a、b波段L.f.與B.f.的兩種效率比隨l的變化。由ir(l)定義知，1r(l)＞1表示L.f.濾出帶性能優(yōu)于B.f.;2r(l)＜1、3r(l)＜1表示L.f.抑制帶性能優(yōu)于B.f.。故當(dāng)1r(l)＞1和2r(l)、3r(l)＜1均得到滿足時(shí)，L.f.濾波效果全面優(yōu)于B.f.。

由表2和表3可見，10～20 d和20～40 d波段L.f.的l臨界值l0分別為24、46;當(dāng)l≥l0時(shí)，L.f.的濾波效果必定優(yōu)于B.f.。

3 應(yīng)用舉例

以拉薩夏季整層熱源(100 hPa至地面)準(zhǔn)雙周和準(zhǔn)一月振蕩分析為例，驗(yàn)證L.f.的應(yīng)用。

3.1 資料和濾波

Zhang et al.(2009)用NCEP/NCAR逐日再分析資料(Kalnay et al.，1996)，按Yanai et al.(1992)倒算法求得了全球1948—2006年逐日大氣熱源數(shù)據(jù)集。本文取格點(diǎn)(90°E，30°N)1950—2006年的序列代表拉薩(91.13°E，29.67°N)整層熱源，記為

表2 10～20 d L.f.與B.f.的效率比ir(l)Table 2 The efficiency ratioir(l)for L.f.and B.f.of 10—20 d

表3 20～40 d L.f.與B.f.的效率比ir(l)Table 3 The efficiency ratioir(l)for L.f.and B.f.of 20—40 d

稱它為ty年年序列;t'為年日序，年總?cè)諗?shù)n'=365(平年)或366(閏年)。

從(14)式中提取序列

稱它為ty年夏季序列;td為夏季日序。

使用L.f.(參數(shù)l取121)從(14)式的年序列中求得夏季a、b波段(即10～20 d，20～40 d波段)振蕩分量。a、b波段濾波結(jié)果記為

相應(yīng)地，使用H.f.直接從夏季序列((15)式)中求得a、b波段振蕩分量

以1959年為例，對(duì)兩個(gè)波段L.f.的濾波結(jié)果與諧波分析結(jié)果進(jìn)行比較。圖4給出了1959年夏季的Q及由(17)、(18)式給出的a、b波段LQ'、HQ'。

圖4 拉薩1959年夏季整層熱源及濾波結(jié)果(圖4a中水平線為ˉQ(1959)=79.5 W·m－2;圖4b、c中實(shí)線分別為L(zhǎng)Q'a和LQ'b，虛線分別為HQ'a和HQ'b;單位:W·m－2)a.Q;b.LQ'a和HQ'a;c.LQ'b和HQ'bFig.4 The whole heat source and its filtering results over Lhasa in summer 1959(The horizontal line indicatesˉQ(1959)=79.5 W·m－2in Fig.4a;The solid lines indicateLQ'aandLQ'b，and the dashed lines indicateHQ'aandHQ'bin Fig.4b and Fig.4c，respectively;units:W·m－2)a.Q;b.LQ'aandHQ'a;c.LQ'bandHQ'b

由圖4可知，通過L.f.濾波得到的a、b波段振蕩分量與諧波分析的結(jié)果非常接近。

3.2 a、b頻段振蕩強(qiáng)度分析

用(16)、(17)式的資料定義了ty年a、b頻段的振蕩強(qiáng)度，即它們是ty年夏季a、b頻段振蕩的日均方差。圖5給出了它們的年際變化曲線。可見，LE、HE做準(zhǔn)同步變化;a、b頻段振蕩存在明顯的年際變化。

由圖4和圖5可知，兩個(gè)波段L.f.的濾波結(jié)果在加熱季內(nèi)變化與波段能量年際變化兩方面均與相應(yīng)諧波分析接近，由此可以說明，選用l=121的L.f.可以很好地濾出逐日熱源序列中的10～20 d與20～40 d振蕩。l=121的L.f.濾波結(jié)果可用于青藏高原大氣加熱季內(nèi)振蕩的實(shí)際分析。

表4據(jù)LEa(ty)、LEb(ty)值列出了a、b頻段的6個(gè)最強(qiáng)年。

3.3 a、b頻段的振蕩過程

圖5 拉薩1950—2006年a、b頻段振蕩強(qiáng)度的年際變化(實(shí)線為L(zhǎng)Ea、LEb;虛線為HEa、HEb;單位:104W2·m－4;1970年LEa、LEb缺損)a.LEa和HEa;b.LEb和HEbFig.5 The interannual variations of oscillation intensity in a and b frequency ranges in Lhasa from 1950 to 2006(The solid lines indicateLEaandLEb，and the dashed lines indicate HEaandHEb;units:104W2·m－4;In 1970，theLEaandLEbis deficient)a.LEa andHEa;b.LEbandHEb

圖6 拉薩a、b頻段振蕩強(qiáng)年振蕩過程(單位:102W·m－2)a.LQ'a(1997，1988，1967年);b.LQ'b(1968，1959，1960年)Fig.6 The oscillation processes in the strong oscillation years in a and b frequency ranges in Lhasa(units:102W·m－2)a.LQ'a(in 1997，1988 and 1967);b.LQ'b(in 1968，1959 and 1960)

表4 拉薩1950—2006年夏季a、b頻段振蕩的最強(qiáng)年份Table 4 The strong oscillation years in a and b frequency ranges in Lhasa in summer from 1950 to 2006

利用(16)式的LQ'a、LQ'b，可以分析拉薩夏季a、b頻段的振蕩過程。圖6給出了a、b頻段振蕩3個(gè)最強(qiáng)年振蕩過程?？梢姡?997、1988年a頻段振蕩均在7月上旬前后達(dá)最強(qiáng)，且1988年8月中、下旬還有一次較強(qiáng)過程;而1967年a頻段振蕩在7月中、下旬前后最強(qiáng)。1968、1959年b頻段振蕩分別在7月上旬、6月中旬前后最強(qiáng)，而在1960年強(qiáng)振蕩則從6月上旬一直維持至8月上旬間。因此，拉薩夏季熱源a、b頻段振蕩過程存在明顯的年際差異和多樣性。

4 結(jié)論

根據(jù)青藏高原大氣熱源季內(nèi)振蕩分析的實(shí)際需要，設(shè)計(jì)了適于提取準(zhǔn)雙周(10～20 d)、準(zhǔn)一月(20～40 d)振蕩的Lanczos濾波器(L.f.)。通過與Butterworth濾波器(B.f.)濾波效果的定量分析，確定了準(zhǔn)雙周、準(zhǔn)一月L.f.窗寬參數(shù)l的臨界值l0=24、46;當(dāng)l≥l0時(shí)，L.f.濾波器性能全面優(yōu)于B.f.。

選用拉薩附近格點(diǎn)(90°E，30°N)的1950—2006年整層大氣熱源資料Q作了應(yīng)用試驗(yàn)，結(jié)果表明，由l=121的L.f.得到的拉薩夏季準(zhǔn)雙周、準(zhǔn)一月振蕩分量LQ'a、LQ'b質(zhì)量可靠，可用于青藏高原大氣熱源季內(nèi)振蕩強(qiáng)度的年際差異和季內(nèi)過程分析。

致謝:張博、陳隆勛先生提供了拉薩整層大氣熱源資料，謹(jǐn)致謝忱。

馮康.1978.數(shù)值計(jì)算方法［M］.北京:國(guó)防工業(yè)出版社.

鞏遠(yuǎn)發(fā)，段廷揚(yáng)，張菡.2007.夏季亞洲大氣熱源匯的變化特征及其與江淮流域旱澇的關(guān)系［J］.大氣科學(xué)，31(1):89-98

李崇銀.2004.大氣季節(jié)內(nèi)振蕩研究的新進(jìn)展［J］.自然科學(xué)進(jìn)展，14(7):734-741.

劉式適，柏晶瑜，徐祥德，等.2000.青藏高原大地形的動(dòng)力、熱力作用與低頻振蕩［J］.應(yīng)用氣象學(xué)報(bào)，11(3):312-321.

彭永清，王盤興，吳洪寶.1997.大氣低頻變化的分析與應(yīng)用［M］.北京:氣象出版社.

燃料化學(xué)工業(yè)部石油地球物理勘探局計(jì)算中心站.1974.地震勘探數(shù)字技術(shù):第二冊(cè)［M］.北京:科學(xué)出版社.

王躍男，陳隆勛，何金海，等.2009.夏季青藏高原熱源低頻振蕩對(duì)我國(guó)東部降水的影響［J］.應(yīng)用氣象學(xué)報(bào)，20(4):419-427.

徐國(guó)強(qiáng)，朱乾根.2000.1998年青藏高原大氣低頻振蕩的結(jié)構(gòu)特征分析［J］.南京氣象學(xué)院學(xué)報(bào)，23(4):505-512.

姚菊香，王盤興，李麗平.2005.季節(jié)內(nèi)振蕩研究中兩種數(shù)字濾波器的性能對(duì)比［J］.南京氣象學(xué)院學(xué)報(bào)，28(2):248-253.

章基嘉，彭永清，王鼎良.1984.夏季青藏高原各熱源分量的時(shí)頻特征及高度場(chǎng)對(duì)它們的響應(yīng)［C］//青藏高原氣象科學(xué)實(shí)驗(yàn)文集(一).北京:科學(xué)出版社:182-192.

章基嘉，孫國(guó)武，陳葆德.1991.青藏高原地區(qū)的大氣低頻振蕩［M］.北京:氣象出版社.

周兵，何金海，徐海明.2000.青藏高原氣象要素場(chǎng)低頻特征及其與夏季區(qū)域降水的聯(lián)系［J］.南京氣象學(xué)院學(xué)報(bào)，23(1):93-100.

Duchon C E.1979.Lanczos filtering in one and two dimensions［J］.J Appl Meteor，18(8):1016-1022.

Kalnay E，Kanamitsu M，Kistler R，et al.1996.The NCEP/NCAR 40-year reanalysis project［J］.Bull Amer Meteor Soc，77(3):437-471.

Madden R A，Julian P R.1971.Detection of a 40-50 day oscillation in the zonal wind in the tropical Pacific［J］.J Atmos Sci，28(5):702-708.

Madden R A，Julian P R.1972.Discription of global-scale circulation cells in the tropics with a 40-50 day period［J］.J Atmos Sci，29(6):1109-1123.

Murakami M.1984.Analysis of deep convective activity over the western Pacific and southeast Asia.Part II:Seasonal and in traseasonal variation during the northern summer［J］.J Meteor Soc Japan，62(1):88-108.

Yanai M，Li C，Song Z S.1992.Seasonal heating of the Tibetan Plateau and its effects on the evolution of the Asian summer monsoon［J］.J Meteor Soc Japan，70(1):319-350.

Zhang Bo，Chen Longxun，He Jinhai，et al.2009.Characteristics of atmospheric heat sources over Asia in summer:Comparison of results calculated using multiple reanalysis datasets［J］.Acta Meteor Sinica，23(5):585-597.

Zhang C.2005.Madden-Julian Oscillation［J］.Reviews of Geophysics，43，RG2003.doi:10.1029/2004RG000158.