粵港澳大灣區(qū)臺(tái)風(fēng)危險(xiǎn)性分析

陳 煜，楊 劍，段忠東，吳甜甜，歐進(jìn)萍

（1.哈爾濱工業(yè)大學(xué)（深圳）土木與環(huán)境工程學(xué)院，廣東 深圳 518055；2.深圳市土木工程智能結(jié)構(gòu)系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，廣東 深圳 518055）

引言

西北太平洋是世界上熱帶氣旋發(fā)生次數(shù)最頻繁、強(qiáng)度最大的區(qū)域，每年生成的臺(tái)風(fēng)總數(shù)一般超過20個(gè)，個(gè)別年份甚至高達(dá)40個(gè)［1］。中國處于西北太平洋的西岸，是世界上受臺(tái)風(fēng)影響最為嚴(yán)重的國家之一?；浉郯拇鬄硡^(qū)位于中國華南地區(qū)，由香港、澳門兩個(gè)特別行政區(qū)和廣東省廣州、深圳、珠海、佛山、惠州、東莞、中山、江門、肇慶等九個(gè)珠三角城市組成，總面積5.6 萬km2。大灣區(qū)作為中國經(jīng)濟(jì)最為發(fā)達(dá)、人口最稠密的沿海地區(qū)，同時(shí)也是受臺(tái)風(fēng)影響最為嚴(yán)重的地區(qū)之一，平均每年大約有2～3 次臺(tái)風(fēng)登陸，臺(tái)風(fēng)引起的強(qiáng)風(fēng)、暴雨、風(fēng)暴潮和洪澇等災(zāi)害會(huì)對(duì)臺(tái)風(fēng)途經(jīng)的地區(qū)造成嚴(yán)重的損失，其中僅2013～2015 年造成的直接經(jīng)濟(jì)損失就超過150 億元［2］。因此，大灣區(qū)的抗風(fēng)設(shè)計(jì)、區(qū)域防災(zāi)規(guī)劃、保險(xiǎn)及再保險(xiǎn)業(yè)都需要對(duì)臺(tái)風(fēng)風(fēng)險(xiǎn)進(jìn)行量化評(píng)估。由于臺(tái)風(fēng)的發(fā)生、移動(dòng)和強(qiáng)度具有隨機(jī)性，大灣區(qū)內(nèi)部的某一局部區(qū)域遭遇臺(tái)風(fēng)卻是小概率事件，而且人類科學(xué)觀測(cè)臺(tái)風(fēng)的歷史也較短，僅通過有限的歷史觀測(cè)資料來進(jìn)行臺(tái)風(fēng)危險(xiǎn)性分析具有較大的局限性，特別是對(duì)長重現(xiàn)期的極值風(fēng)速無法獲得可靠的估計(jì)。

工程抗風(fēng)設(shè)計(jì)和臺(tái)風(fēng)巨災(zāi)風(fēng)險(xiǎn)分析往往需要足夠時(shí)間長度的臺(tái)風(fēng)樣本來估算臺(tái)風(fēng)危險(xiǎn)性。然而復(fù)雜的臺(tái)風(fēng)氣象數(shù)值模擬以其高昂的計(jì)算成本和耗時(shí)限制了我們?cè)诙虝r(shí)間內(nèi)模擬大量臺(tái)風(fēng)樣本的能力，所以工程界利用臺(tái)風(fēng)參數(shù)化數(shù)值模型，結(jié)合歷史臺(tái)風(fēng)觀測(cè)資料的統(tǒng)計(jì)建模，模擬大量反映臺(tái)風(fēng)統(tǒng)計(jì)規(guī)律的臺(tái)風(fēng)事件，并據(jù)此進(jìn)行臺(tái)風(fēng)危險(xiǎn)性分析已成為國際上展開此類研究普遍接受的方法［3-23］。傳統(tǒng)的臺(tái)風(fēng)事件模擬方法主要包括：分區(qū)圓法［5-10］和全路徑法［11-23］。分區(qū)圓法是以目標(biāo)點(diǎn)為中心，向外擴(kuò)展成一個(gè)半徑為200～300 km 的圓形局部區(qū)域，然后從經(jīng)過該圓形區(qū)域的歷史臺(tái)風(fēng)樣本中提取移動(dòng)和強(qiáng)度參數(shù)的概率分布，采用Monte Carlo 方法隨機(jī)模擬該區(qū)域的臺(tái)風(fēng)事件。這種方法適用于臺(tái)風(fēng)樣本充足的局部小區(qū)域，但對(duì)于多城市的大區(qū)域以及長距離分布的基礎(chǔ)設(shè)施，如鐵路、公路和電網(wǎng)系統(tǒng)等并不適用。為了克服上述缺點(diǎn)，全路徑法將臺(tái)風(fēng)活動(dòng)分為：生成、移動(dòng)、強(qiáng)度等部分，然后分段模擬，最后組合成完整的臺(tái)風(fēng)事件。這種方法可以模擬臺(tái)風(fēng)從生成到消亡的路徑和強(qiáng)度，從而對(duì)大范圍區(qū)域進(jìn)行臺(tái)風(fēng)危險(xiǎn)性評(píng)估。

具體來說，首先對(duì)臺(tái)風(fēng)生成位置進(jìn)行模擬，主要方法包括：直接抽樣法和估算概率模擬法。直接抽樣法［11］是僅從歷史生成位置樣本中隨機(jī)抽樣，這會(huì)導(dǎo)致初始模擬位置分布范圍非常有限。估算概率模擬法首先估算臺(tái)風(fēng)生成位置的概率，然后通過Monte Carlo 隨機(jī)模擬生成位置，可覆蓋整個(gè)海洋區(qū)域。Rumpf等［13-14］、Emanuel 等［15］、Hall 和Jewson［16］使用歷史生成點(diǎn)，通過高斯核密度函數(shù)估算海洋中每個(gè)位置的生成概率。Emanuel［24］、Tippet和Camargo［25］、Lee 等［26-27］、Chen 和Duan［28］通過建立臺(tái)風(fēng)生成位置與海氣環(huán)境參數(shù)之間的關(guān)系，將全球氣候模式模擬的海氣環(huán)境參數(shù)輸入到該關(guān)系中以模擬臺(tái)風(fēng)的生成位置，這種方法可應(yīng)用于氣候變化對(duì)臺(tái)風(fēng)活動(dòng)的影響研究中。其次是對(duì)臺(tái)風(fēng)移動(dòng)路徑的模擬，一種方法是通過歷史臺(tái)風(fēng)的移動(dòng)速度和方向來模擬臺(tái)風(fēng)的路徑［11，13，16-17，20］。另一種方法是建立臺(tái)風(fēng)移動(dòng)與環(huán)境氣流和β 漂移項(xiàng)之間的關(guān)系［15，27-28］。由于環(huán)境氣流可以很容易地從再分析數(shù)據(jù)中獲得，所以這種方法更適用于歷史臺(tái)風(fēng)數(shù)據(jù)不足的地區(qū)。對(duì)于臺(tái)風(fēng)強(qiáng)度的模擬也有2 種方法：統(tǒng)計(jì)方法和物理方法。前者是基于歷史臺(tái)風(fēng)強(qiáng)度樣本的統(tǒng)計(jì)強(qiáng)度模型，如自回歸模型［11，17，20］、馬爾科夫模型［14］等；后者是基于臺(tái)風(fēng)物理機(jī)制的簡化數(shù)值強(qiáng)度模型［15，29］，不過計(jì)算耗時(shí)也相對(duì)較大。

將模擬的臺(tái)風(fēng)事件參數(shù)輸入到臺(tái)風(fēng)風(fēng)場模型可以獲得臺(tái)風(fēng)影響區(qū)域各位置的風(fēng)速。工程上為了快速地計(jì)算臺(tái)風(fēng)風(fēng)場，一般采用參數(shù)化的邊界層動(dòng)力學(xué)臺(tái)風(fēng)風(fēng)場模型［7，30-34］。參數(shù)化風(fēng)場模型最早由Chow［35］提出，它是由一組流體動(dòng)量方程構(gòu)成的均勻厚度的水平二維大氣邊界層模型。該模型隨后被進(jìn)一步修正和改進(jìn)并用于評(píng)估臺(tái)風(fēng)危險(xiǎn)性［17，32］。這類二維風(fēng)場模型只能模擬邊界層內(nèi)沿高度平均的風(fēng)速，然后通過經(jīng)驗(yàn)的風(fēng)速折減關(guān)系得到邊界層內(nèi)的風(fēng)速垂直結(jié)構(gòu)，該方法不能完全反映臺(tái)風(fēng)風(fēng)場在空間分布上的特性。為了更加精細(xì)地模擬臺(tái)風(fēng)風(fēng)場結(jié)構(gòu)，Meng等［31］、Kepert［33-34］通過對(duì)邊界層內(nèi)N-S方程的離散簡化，建立了沿高度解析的三維參數(shù)化風(fēng)場模型，不過他們忽略了臺(tái)風(fēng)的垂直平流和垂直擴(kuò)散過程，這將導(dǎo)致邊界層內(nèi)臺(tái)風(fēng)風(fēng)場的超梯度現(xiàn)象比觀測(cè)結(jié)果弱。事實(shí)上，臺(tái)風(fēng)登陸后由于陸地表面呈現(xiàn)出非均勻性，會(huì)導(dǎo)致地表粗糙度和局地地形發(fā)生劇烈變化，這會(huì)對(duì)臺(tái)風(fēng)邊界層的垂直運(yùn)動(dòng)產(chǎn)生影響，進(jìn)而影響風(fēng)場結(jié)構(gòu)和風(fēng)速計(jì)算的準(zhǔn)確性。

總的來說，全路徑合成臺(tái)風(fēng)模型可以模擬臺(tái)風(fēng)從生成到消亡的路徑和強(qiáng)度，與臺(tái)風(fēng)風(fēng)場模型相結(jié)合便可在大區(qū)域上進(jìn)行臺(tái)風(fēng)危險(xiǎn)性評(píng)估，從而量化粵港澳大灣區(qū)不同城市的臺(tái)風(fēng)危險(xiǎn)性。在文中，我們將基于統(tǒng)計(jì)動(dòng)力學(xué)全路徑隨機(jī)合成臺(tái)風(fēng)方法，模擬臺(tái)風(fēng)事件的發(fā)生、路徑、強(qiáng)度。與Chen 和Duan［28］建立的全路徑模型相比，文中的全路徑合成臺(tái)風(fēng)模型在建模數(shù)據(jù)、生成、路徑和強(qiáng)度的模擬等方面均有所改進(jìn)。此外，本文采用Yang 等［36］開發(fā)的三維數(shù)值風(fēng)場模型模擬臺(tái)風(fēng)過境期間大灣區(qū)各地區(qū)的風(fēng)速。與Meng 等［31］、Kepert［33-34］的三維風(fēng)場模型相比，Yang等［36］三維數(shù)值風(fēng)場模型考慮了垂直平流和垂直擴(kuò)散過程對(duì)風(fēng)場的影響，從而可以更好地模擬風(fēng)場的空間結(jié)構(gòu)。另外，本文將從防災(zāi)減災(zāi)的角度對(duì)模擬的臺(tái)風(fēng)事件集進(jìn)行綜合檢驗(yàn)，包括大灣區(qū)登陸臺(tái)風(fēng)的強(qiáng)度、頻率、風(fēng)場，驗(yàn)證模型對(duì)大灣區(qū)歷史臺(tái)風(fēng)統(tǒng)計(jì)特性的模擬能力。最后通過改進(jìn)的全路徑合成臺(tái)風(fēng)模型生成萬年臺(tái)風(fēng)事件集，與三維數(shù)值風(fēng)場模型相結(jié)合模擬大灣區(qū)的年最大風(fēng)速，并分析大灣區(qū)內(nèi)各個(gè)區(qū)域的臺(tái)風(fēng)危險(xiǎn)性。除非特別說明，文中的臺(tái)風(fēng)和熱帶氣旋不進(jìn)行區(qū)別。

1 建模數(shù)據(jù)

臺(tái)風(fēng)危險(xiǎn)性建模和分析所需的數(shù)據(jù)包括臺(tái)風(fēng)最佳路徑數(shù)據(jù)集、全球環(huán)境再分析數(shù)據(jù)、全球陸地高程數(shù)據(jù)、全球土地覆蓋數(shù)據(jù)等?；浉郯拇鬄硡^(qū)位于西北太平洋的西岸，大灣區(qū)遭受到的臺(tái)風(fēng)全部都來自西北太平洋（包含南海）。常用的西北太平洋臺(tái)風(fēng)最佳路徑數(shù)據(jù)集分別來自中國氣象局上海臺(tái)風(fēng)研究所（CMA）［37］、日本氣象廳（JMA）和美國聯(lián)合臺(tái)風(fēng)預(yù)警中心（JTWC）。由于判別臺(tái)風(fēng)事件的方法和標(biāo)準(zhǔn)不同，這3個(gè)機(jī)構(gòu)的最佳路徑數(shù)據(jù)集在臺(tái)風(fēng)數(shù)量、強(qiáng)度等參數(shù)的記錄上均存在差異。對(duì)于登陸中國的臺(tái)風(fēng)，CMA 的記錄相較于其他2個(gè)最佳路徑數(shù)據(jù)集更可靠。CMA 最佳路徑數(shù)據(jù)記錄了1949-2019年所有臺(tái)風(fēng)每6 h的中心位置和強(qiáng)度。由于1980 年以后衛(wèi)星觀測(cè)成為臺(tái)風(fēng)的主流觀測(cè)手段，所以1980 年之后的數(shù)據(jù)質(zhì)量相較于之前的數(shù)據(jù)更可靠，所以本文使用CMA最佳路徑數(shù)據(jù)集（1980-2019）作為建模初始數(shù)據(jù)，最佳路徑數(shù)據(jù)集的詳細(xì)描述如表1所示。

表1 CMA最佳路徑數(shù)據(jù)集說明Table 1 The best track data of CMA

全球環(huán)境參數(shù)數(shù)據(jù)來自美國國家環(huán)境預(yù)測(cè)中心（NCEP）/國家大氣研究中心（NCAR）全球再分析數(shù)據(jù)［38］和COBE-SST 數(shù)據(jù)［39］，包含1980-2019 年每6 h 的大氣環(huán)境風(fēng)速度、月平均大氣環(huán)境溫度、月平均大氣比濕和月平均海洋表面溫度等，如表2所示。NCEP/NCAR全球再分析數(shù)據(jù)的大氣環(huán)境風(fēng)速和溫度包含了從1 000 hPa至10 hPa，共17個(gè)氣壓層的數(shù)據(jù)，比濕只有底部的8個(gè)氣壓層數(shù)據(jù)。由于臺(tái)風(fēng)中心的位置不一定會(huì)落在再分析數(shù)據(jù)的網(wǎng)格點(diǎn)上，所以在模擬中我們將再分析數(shù)據(jù)的環(huán)境參數(shù)線性插值到臺(tái)風(fēng)中心的位置和時(shí)間上。

表2 全球環(huán)境再分析數(shù)據(jù)Table 2 The global environmental reanalysis data

為了考慮地形對(duì)臺(tái)風(fēng)風(fēng)場的抬升和沉降作用，本研究使用全球大洋地勢(shì)圖GEBCO（General Bathymetric Chart of the Oceans，水平空間分辨率為15 km×15 km）和土地覆蓋數(shù)據(jù)ESA-CCI（European Space Agency-Climate Change Initiative，水平空間分辨率為300 m×300 m），建立粵港澳大灣區(qū)的地形和地表粗糙長度空間分布。由于地面覆蓋物（植被、森林等）在不同季節(jié)的變化會(huì)導(dǎo)致地表粗糙長度發(fā)生變化，考慮到影響大灣區(qū)的臺(tái)風(fēng)多發(fā)生在夏季，文中利用表3中ESA-CCI 提供的土地覆蓋類型與地表粗糙長度對(duì)應(yīng)關(guān)系確定大灣區(qū)土地粗糙長度的空間分布［36］。此外，考慮到ESA-CCI 土地覆蓋數(shù)據(jù)分辨率較高且風(fēng)速易受周圍場地的影響，為了與文中所用的三維數(shù)值臺(tái)風(fēng)模型匹配（風(fēng)場模型分辨率為km量級(jí)），本研究將ESA-CCI地表粗糙長度按照1.2 km×1.2 km的窗寬進(jìn)行平滑。

表3 土地覆蓋類型與地表粗糙長度關(guān)系Table 3 Relationship between land cover type and surface roughness length

2 臺(tái)風(fēng)危險(xiǎn)性分析模型

臺(tái)風(fēng)危險(xiǎn)性模型包括全路徑隨機(jī)合成臺(tái)風(fēng)模型和臺(tái)風(fēng)風(fēng)場模型，其中全路徑隨機(jī)合成臺(tái)風(fēng)模型可細(xì)分為生成模型、移動(dòng)模型、強(qiáng)度模型3 個(gè)部分。首先，我們對(duì)西北太平洋的統(tǒng)計(jì)動(dòng)力學(xué)全路徑合成臺(tái)風(fēng)模型［28］進(jìn)行改進(jìn)和驗(yàn)證；然后利用該模型生成西北太平洋區(qū)域的萬年事件集，臺(tái)風(fēng)事件集參數(shù)包括模擬生成的臺(tái)風(fēng)每6 h 的時(shí)間、路徑、強(qiáng)度等參數(shù)。然后篩選出萬年事件集中影響粵港澳大灣區(qū)的臺(tái)風(fēng)事件，將模擬的臺(tái)風(fēng)事件參數(shù)輸入臺(tái)風(fēng)三維數(shù)值風(fēng)場模型，估算大灣區(qū)內(nèi)各地區(qū)的臺(tái)風(fēng)風(fēng)速。

2.1 生成模型

本研究通過建立生成模型模擬臺(tái)風(fēng)的年生成數(shù)量、生成的位置以及時(shí)間。通過對(duì)歷史臺(tái)風(fēng)的年發(fā)生數(shù)量進(jìn)行K-S 檢驗(yàn)，我們發(fā)現(xiàn)在5%的顯著性水平下，臺(tái)風(fēng)的年生成個(gè)數(shù)服從泊松分布。于是，不同于Chen 和Duan 模型［28］中所使用的核密度方法，這里我們直接簡單地假設(shè)西北太平洋的臺(tái)風(fēng)年發(fā)生數(shù)分布為泊松分布，如式（1）所示：

式中，k是年生成個(gè)數(shù)；λ是年平均生成個(gè)數(shù)，通過歷史最佳路徑數(shù)據(jù)集確定。本研究根據(jù)上述估算的年生成數(shù)概率分布，通過Monte Carlo法隨機(jī)抽樣臺(tái)風(fēng)的年生成個(gè)數(shù)。

對(duì)于臺(tái)風(fēng)初始生成的位置和時(shí)間，我們首先將臺(tái)風(fēng)首次達(dá)到15 m/s時(shí)的位置作為臺(tái)風(fēng)的初始位置，采用高斯核密度法估算臺(tái)風(fēng)初始生成位置的三維（臺(tái)風(fēng)中心位置的經(jīng)度、緯度、時(shí)間）時(shí)空概率密度分布，其概率密度估算如下式所示：

式中，x是生成位置向量；xi是生成位置的歷史樣本向量；n是樣本大??；S是生成位置的標(biāo)準(zhǔn)差矩陣；σxx、σyy、σzz分別為經(jīng)度、緯度、時(shí)間的方差；γ1、γ2、γ3是標(biāo)準(zhǔn)化后的3個(gè)變量之間相關(guān)系數(shù)的特征向量；wi是生成概率密度的權(quán)重；λ1、λ2、λ3分別為特征值；hotp1、hotp2、hotp3是最佳帶寬，通過有偏交叉驗(yàn)證方法確定［41］。同樣地，本研究根據(jù)估算的生成位置的概率密度分布，通過Monte Carlo法隨機(jī)抽樣西北太平洋海域的臺(tái)風(fēng)生成位置。

2.2 移動(dòng)模型

臺(tái)風(fēng)的移動(dòng)路徑主要受大尺度環(huán)境氣流和β漂流控制。于是，本研究將臺(tái)風(fēng)移動(dòng)速度分解為大尺度環(huán)境氣流引起的引導(dǎo)速度分量和臺(tái)風(fēng)與環(huán)境大氣相互作用的β漂流分量，建立臺(tái)風(fēng)移動(dòng)速度與引導(dǎo)速度和β漂流之間的關(guān)系。引導(dǎo)氣流速度通過NCEP/NCAR 再分析數(shù)據(jù)的每6小時(shí)大氣環(huán)境風(fēng)速進(jìn)行估算。另外，為了進(jìn)一步在模擬移動(dòng)速度中考慮其隨機(jī)變化，這里β漂流采用帶有隨機(jī)噪聲的自回歸模型進(jìn)行建模，而不是像Chen和Duan模型［28］那樣采用區(qū)域平均值，后者會(huì)導(dǎo)致更平滑的模擬結(jié)果。具體的移動(dòng)速度方程如式（3）所示：

式中，Ui和Vi分別是當(dāng)前時(shí)刻的緯向和經(jīng)向移動(dòng)速度；i和i-1 分別表示當(dāng)前和上一個(gè)6 h 的時(shí)刻；Usteer，i和Vsteer，i是當(dāng)前時(shí)刻引導(dǎo)氣流的緯向和經(jīng)向速度，定義為250 至850 hPa 氣壓層距臺(tái)風(fēng)中心5°環(huán)帶上的平均風(fēng)速；γu和γv分別是移動(dòng)速度的回歸系數(shù)；βx和βy是緯向和經(jīng)向的β漂流速度；a0，a1，b0，b1是β漂流的自回歸系數(shù)；εx和εy分別為緯向和經(jīng)向隨機(jī)項(xiàng)。此外，為了消除再分析環(huán)境風(fēng)場數(shù)據(jù)中的歷史氣旋對(duì)引導(dǎo)氣流計(jì)算的干擾，本研究采用Kurihara等［42］提出的濾波方法去除渦旋。

2.3 強(qiáng)度模型

強(qiáng)度模型包括海上強(qiáng)度模型和陸地衰減模型兩部分。對(duì)于海上強(qiáng)度模型，本研究以臺(tái)風(fēng)中心10 m高度處的最大風(fēng)速作為臺(tái)風(fēng)強(qiáng)度的度量，建立臺(tái)風(fēng)海上自回歸強(qiáng)度模型。首先，我們以相對(duì)強(qiáng)度對(duì)數(shù)的6 h變化量作為因變量，自變量包括氣候和持久性變量，以及從環(huán)境再分析數(shù)據(jù)中提取的大氣和海洋環(huán)境參數(shù)，如海溫、相對(duì)濕度、風(fēng)切變等。與Chen和Duan模型［28］相比，本研究考慮了相對(duì)濕度、風(fēng)切變等對(duì)臺(tái)風(fēng)強(qiáng)度的影響因素，這些之前模型沒有考慮的變量可能會(huì)對(duì)強(qiáng)度的模擬結(jié)果產(chǎn)生較大地影響。此外，這些變量在6 h內(nèi)的變化也作為待選自變量。通過逐步回歸方法，從候選的自變量中選出最佳的自變量以保證回歸方程的顯著性。海上強(qiáng)度模型如式（4）所示，其中式（4）的第2 個(gè)方程估算第2 時(shí)刻的強(qiáng)度，式（4）的第1 個(gè)方程估算第3及隨后時(shí)刻的強(qiáng)度。

式中，i+1、i-1、i分別為下一個(gè)6 h 時(shí)刻、前一個(gè)6 h時(shí)刻和當(dāng)前時(shí)刻；LR是線性回歸算子；RI為相對(duì)強(qiáng)度，RI=Vmax/PI，Vmax是最大風(fēng)速，PI是潛在強(qiáng)度［27］；SST是月平均海洋表面溫度；VS是月平均垂直風(fēng)切變；RH是月平均相對(duì)濕度；U，V分別是緯向和經(jīng)向的移動(dòng)速度；ε是海洋強(qiáng)度的隨機(jī)項(xiàng)。

本研究同時(shí)使用中心氣壓差來定義臺(tái)風(fēng)的強(qiáng)度，中心氣壓差將用在登陸臺(tái)風(fēng)強(qiáng)度模型和臺(tái)風(fēng)風(fēng)場模型中。我們分別建立海上和陸地的中心壓差與最大風(fēng)速之間的關(guān)系，如式（5）所示，其可決系數(shù)R2分別為0.97和0.91。

其中，ΔPocean和ΔPland分別是海洋和陸地的中心壓差。臺(tái)風(fēng)一旦登陸，其水汽熱源會(huì)被切斷，加上地面的摩擦，臺(tái)風(fēng)的強(qiáng)度會(huì)逐漸減弱。最佳路徑數(shù)據(jù)集中的登陸強(qiáng)度記錄表明，登陸強(qiáng)度隨時(shí)間的衰減變化呈現(xiàn)自然指數(shù)遞減關(guān)系。由于臺(tái)風(fēng)登陸后因地形不同導(dǎo)致其衰減速度也會(huì)不同，所以我們用大灣區(qū)和地形相似的廣東省其他區(qū)域臺(tái)風(fēng)登陸數(shù)據(jù)來估計(jì)登陸衰減系數(shù)。于是，本研究基于自然指數(shù)函數(shù)建立登陸廣東省臺(tái)風(fēng)的強(qiáng)度衰減模型以估算登陸之后的強(qiáng)度，具體方程如式（6）所示，

式中，ΔP（t）為臺(tái)風(fēng)登陸后t時(shí)刻的中心壓差；ΔP0為臺(tái)風(fēng)登陸時(shí)的中心壓差；a是衰減系數(shù)；a0和a1是擬合常數(shù)，使用登陸廣東省的臺(tái)風(fēng)強(qiáng)度進(jìn)行擬合；ε是正態(tài)隨機(jī)殘差項(xiàng)。其他區(qū)域使用Chen 和Duan 模型［28］中的k近鄰法估算登陸強(qiáng)度。

2.4 臺(tái)風(fēng)風(fēng)場模型

為了估算臺(tái)風(fēng)產(chǎn)生的風(fēng)速，本研究采用Yang等［36］開發(fā)的三維數(shù)值風(fēng)場模型來計(jì)算臺(tái)風(fēng)產(chǎn)生的風(fēng)速。臺(tái)風(fēng)風(fēng)場中的垂直平流和垂直擴(kuò)散過程是合理模擬三維風(fēng)場的基本要素，Yang等［36］開發(fā)的三維數(shù)值風(fēng)場模型合理地考慮了垂直平流和垂直擴(kuò)散過程對(duì)風(fēng)場的影響，從而可以更好地模擬風(fēng)場的空間結(jié)構(gòu)。臺(tái)風(fēng)風(fēng)場模型的具體介紹參考Yang 等［36］，文中不再贅述。此外，傳統(tǒng)的臺(tái)風(fēng)風(fēng)場［30-32］一般是基于下墊面為海洋的臺(tái)風(fēng)風(fēng)場進(jìn)行開發(fā)的，對(duì)于我們關(guān)心的粵港澳大灣區(qū)陸上城市，需要考慮陸地下墊面對(duì)登陸臺(tái)風(fēng)風(fēng)場的影響。為了在臺(tái)風(fēng)風(fēng)場模型中引入陸地地形效應(yīng)，本研究考慮了地形對(duì)氣流的抬升和沉降作用，地形導(dǎo)致的垂直風(fēng)速可近似表示為：

式中，V為水平風(fēng)速矢量；H為地形高度，使用全球數(shù)值高程GEBCO 獲取。我們將地形產(chǎn)生的垂直風(fēng)速疊加到風(fēng)場模型計(jì)算的垂直風(fēng)速上，代入風(fēng)場模型進(jìn)行迭代計(jì)算，直至收斂。另外，臺(tái)風(fēng)風(fēng)場模型中的摩阻力系數(shù)通過考慮地形和土地覆蓋類型的粗糙長度進(jìn)行計(jì)算［40］，如下式所示，

式中，z1是距離地面的高度；z0為表面粗糙度長度；κ是卡曼常數(shù)，取0.4。我們將全路徑合成臺(tái)風(fēng)模型生成的臺(tái)風(fēng)事件關(guān)鍵參數(shù)，包括臺(tái)風(fēng)的中心氣壓差、中心位置、移動(dòng)速度、移動(dòng)方向等參數(shù)輸入到臺(tái)風(fēng)三維數(shù)值風(fēng)場模型中，便可估算出臺(tái)風(fēng)產(chǎn)生的風(fēng)速。

3 模型驗(yàn)證

本研究通過改進(jìn)的全路徑合成臺(tái)風(fēng)模型生成了萬年臺(tái)風(fēng)事件集，包含大約300 000 個(gè)模擬臺(tái)風(fēng)事件。接下來，我們對(duì)模擬的萬年臺(tái)風(fēng)事件集進(jìn)行綜合檢驗(yàn)，包括臺(tái)風(fēng)的頻率、強(qiáng)度、風(fēng)場結(jié)構(gòu)等，以驗(yàn)證本模型模擬西北太平洋臺(tái)風(fēng)事件的能力。本研究分別從西北太平洋海域和登陸大灣區(qū)臺(tái)風(fēng)2個(gè)角度對(duì)模擬臺(tái)風(fēng)事件集進(jìn)行驗(yàn)證。首先，我們對(duì)比西北太平洋模擬臺(tái)風(fēng)與CMA 最佳路徑數(shù)據(jù)集的年發(fā)生頻率。圖1 顯示模擬的臺(tái)風(fēng)年平均生成數(shù)量與歷史最佳路徑數(shù)據(jù)集記錄均集中在每年24～35個(gè)之間。通過對(duì)模擬和歷史臺(tái)風(fēng)的年發(fā)生數(shù)量進(jìn)行K-S 檢驗(yàn)，結(jié)果表明在5%顯著性水平下模擬和歷史記錄的年發(fā)生數(shù)分布一致。此外，圖2還比較了西北太平洋的模擬和歷史各月平均生成數(shù)量，可以看出模擬結(jié)果與歷史記錄相當(dāng)吻合，最大偏差出現(xiàn)在8月份，但也不超過2.5%。

圖1 西北太平洋的模擬和歷史臺(tái)風(fēng)年生成個(gè)數(shù)分布Fig.1 Distribution of simulated and historical annual numbers in the western North Pacific

圖2 西北太平洋的模擬和歷史臺(tái)風(fēng)各月平均生成個(gè)數(shù)Fig.2 The monthly average number of simulated and historical typhoon in the western North Pacific

本研究進(jìn)一步將西北太平洋海域劃分為2.5°×2.5°的網(wǎng)格，分別統(tǒng)計(jì)各個(gè)網(wǎng)格內(nèi)的臺(tái)風(fēng)參數(shù)，包括年平均通過個(gè)數(shù)、平均移動(dòng)速度、平均移動(dòng)方向（方向以正北為起始，按順時(shí)針為正方向）以及平均最大強(qiáng)度。對(duì)比圖3（a）與（b）可以發(fā)現(xiàn)，西北太平洋的模擬與歷史年平均通過個(gè)數(shù)整體分布一致，二者之間的相關(guān)系數(shù)為0.97。值得注意的是，經(jīng)過日本東部海域的模擬數(shù)量比歷史記錄偏高，這可能是由于我們的模型沒有考慮日本東部副熱帶高壓對(duì)引導(dǎo)氣流的干擾，但這并不影響大灣區(qū)臺(tái)風(fēng)危險(xiǎn)性的評(píng)估。圖3（c）～（f）分別展示了歷史和模擬平均移動(dòng)速度和平均移動(dòng)方向，可以看出歷史與模擬臺(tái)風(fēng)平均移動(dòng)速度和方向整體分布一致，二者相關(guān)系數(shù)分別為0.96和0.97。最后，圖（g）和（h）對(duì)比了每個(gè)網(wǎng)格內(nèi)歷史和模擬臺(tái)風(fēng)的平均最大強(qiáng)度，雖然在某些區(qū)域內(nèi)兩者強(qiáng)度略有差異，但是二者整體分布仍較為吻合，生命最大強(qiáng)度都主要出現(xiàn)在120°E～152.5°E、12.5°N～32.5°N海域內(nèi)，二者的相關(guān)系數(shù)為0.90。

圖3 西北太平洋海域2.5°×2.5°網(wǎng)格內(nèi)的模擬和歷史臺(tái)風(fēng)參數(shù)對(duì)比Fig.3 Comparison of simulated and historical typhoon parameters in 2.5°×2.5°grids of the western North Pacific

臺(tái)風(fēng)造成的財(cái)產(chǎn)和經(jīng)濟(jì)損失大部分都是臺(tái)風(fēng)在近岸或登陸之后造成的，所以本研究更關(guān)心本模型對(duì)登陸粵港澳大灣區(qū)的模擬結(jié)果?；浉郯拇鬄硡^(qū)的地理位置如圖4 所示，本研究對(duì)比了模擬與歷史登陸大灣區(qū)臺(tái)風(fēng)事件的統(tǒng)計(jì)結(jié)果。如圖5所示，從模擬和歷史臺(tái)風(fēng)的年登陸個(gè)數(shù)頻率分布來看，模擬臺(tái)風(fēng)的年登陸個(gè)數(shù)比歷史最佳路徑數(shù)據(jù)集的統(tǒng)計(jì)結(jié)果略微偏高。不過從臺(tái)風(fēng)的年平均登陸個(gè)數(shù)來看，模擬和歷史平均值相當(dāng)，分別為1.27和1.21。

圖4 粵港澳大灣區(qū)的地理位置（粉色線部分）Fig.4 Geographical location of the Guangdong-HongKong-Macao Greater Bay Area（the pink line）

圖5 大灣區(qū)模擬和歷史臺(tái)風(fēng)的年登陸個(gè)數(shù)分布Fig.5 Distribution of annual landfalls of simulated and historical typhoon in the Greater Bay Area

圖6（a）～（b）對(duì)比了模擬和歷史臺(tái)風(fēng)登陸粵港澳大灣區(qū)的移動(dòng)速度和中心氣壓差。從圖6（a）登陸移動(dòng)速度來看，登陸大灣區(qū)之后的臺(tái)風(fēng)移動(dòng)速度較為吻合，模擬和歷史移動(dòng)速度的大部分都落在2～6 m/s 區(qū)間，其中超過三分之一登陸臺(tái)風(fēng)的移動(dòng)速度小于4 m/s，這些低移動(dòng)速度的臺(tái)風(fēng)可能會(huì)帶來長時(shí)間的集中降雨，從而造成嚴(yán)重的洪澇災(zāi)害。圖6（b）表明登陸后的模擬和歷史臺(tái)風(fēng)中心氣壓差大部分都落在10～30 hPa 區(qū)間，相當(dāng)于熱帶低壓到熱帶風(fēng)暴，這說明臺(tái)風(fēng)登陸大灣區(qū)之后強(qiáng)度都會(huì)迅速地衰減。

圖6 大灣區(qū)模擬和歷史臺(tái)風(fēng)登陸大灣區(qū)的參數(shù)分布Fig.6 Distribution of parameters of simulation and historical typhoon landing in the Greater Bay Area

為了進(jìn)一步驗(yàn)證登陸臺(tái)風(fēng)的登陸特性，本研究將大灣區(qū)沿海岸線以每100 km 分作8 段，海岸線分段如圖7 所示。為了更好地考慮臺(tái)風(fēng)對(duì)粵港澳大灣區(qū)全域的影響，我們將大灣區(qū)海岸線向北和向南分別延伸到汕頭和湛江。這里，我們將模擬的萬年臺(tái)風(fēng)事件集以40 年為一組分作250 組，分別統(tǒng)計(jì)了各海岸線段模擬登陸臺(tái)風(fēng)參數(shù)的整體平均值、中位值、90%區(qū)間。圖8（a）～（b）分別對(duì)比了模擬和歷史臺(tái)風(fēng)登陸大灣區(qū)海岸線的4 個(gè)登陸特性參數(shù)，包括年登陸臺(tái)風(fēng)個(gè)數(shù)、登陸臺(tái)風(fēng)的移動(dòng)速度、移動(dòng)方向以及登陸強(qiáng)度。從圖8（a）中的年登陸臺(tái)風(fēng)個(gè)數(shù)來看，大灣區(qū)內(nèi)的海岸線第4、5段的模擬平均結(jié)果比歷史平均值偏大，而第6段則偏小。不過從模擬結(jié)果的90%區(qū)間來看，所有海岸線段的歷史年平均登陸個(gè)數(shù)均能落在模擬結(jié)果的90%區(qū)間內(nèi)。同樣地，從圖8（b）～（d）中可以看出，所有海岸線段內(nèi)的歷史移動(dòng)速度、移動(dòng)方向以及登陸強(qiáng)度的平均值均能落在模擬結(jié)果的90%區(qū)間內(nèi)。

圖7 粵港澳大灣區(qū)海岸線分段Fig.7 Coastline segmentation of the Guangdong-Hong Kong-Macao Greater Bay Area

圖8 模擬和歷史臺(tái)風(fēng)登陸大灣區(qū)海岸線的參數(shù)分布Fig.8 Parameter distribution of simulated and historical typhoon landing on the coastline of the Greater Bay Area

最后，為了驗(yàn)證采用的三維數(shù)值風(fēng)場模型對(duì)登陸大灣區(qū)臺(tái)風(fēng)風(fēng)場的模擬能力，本研究使用ESA-CCI 土地覆蓋數(shù)據(jù)獲取大灣區(qū)的土地覆蓋類型，通過表3 中土地覆蓋類型與地表粗糙長度關(guān)系確定大灣區(qū)的粗糙長度空間分布。由于ESA-CCI土地覆蓋數(shù)據(jù)的分辨率與風(fēng)場模型的分辨率不同，為了與三維數(shù)值臺(tái)風(fēng)模型相匹配，本研究將ESA-CCI 地表粗糙長度按照1.2 km×1.2 km 的窗寬進(jìn)行平滑，如圖9 所示。文中模擬了歷史臺(tái)風(fēng)“黑格比”（Hagupit）于2008 年9 月23 日1 時(shí)～9 月24 日12 時(shí)（UTC）的逐小時(shí)風(fēng)場，并與大灣區(qū)內(nèi)觀測(cè)站記錄該時(shí)段的風(fēng)速和風(fēng)向進(jìn)行對(duì)比。圖10（a）和（b）分別對(duì)比了臺(tái)風(fēng)“黑格比”（Hagupit）在大灣區(qū)的江門市上川島觀測(cè)站和珠海市斗門觀測(cè)站的模擬和實(shí)測(cè)結(jié)果，包括逐小時(shí)風(fēng)速（10 m高度處的10 min平均風(fēng)速）和風(fēng)向（角度以正北為起始，順時(shí)針為正方向）。

圖9 大灣區(qū)地表粗糙長度空間分布Fig.9 Spatial distribution of surface roughness length in the Greater Bay Area

上川島站位于廣東省江門市臺(tái)山上川鎮(zhèn)，該站點(diǎn)西側(cè)臨海，東側(cè)有大塘山、白灰山等環(huán)繞，站點(diǎn)附近土地覆蓋類型為：農(nóng)作物、水體和樹木等。斗門站位于珠海市斗門區(qū)，該站點(diǎn)位于郊區(qū)，四周建筑物和高大植被都比較少，站點(diǎn)附近土地覆蓋類型為：城市區(qū)、農(nóng)作物和水體等。上川島站和斗門站采用平滑后的地表粗糙長度，分別為：0.21 m和0.29 m。肖玉鳳等［3］使用沿高度平均風(fēng)場模型對(duì)上川島站風(fēng)速進(jìn)行模擬時(shí)，認(rèn)為上川島站為B類地貌，故使用地表粗糙長度為0.02 m，與文本方法取值差距較大。這主要是因?yàn)楸狙芯繉?duì)ESACCI 地表粗糙長度的空間分布進(jìn)行了平滑處理，所以上川島站的地表粗糙長度受到周圍場地的影響。不過從圖10 的對(duì)比結(jié)果來看，本研究模擬的2 個(gè)觀測(cè)站點(diǎn)的風(fēng)速和風(fēng)向與實(shí)測(cè)結(jié)果較為一致，表明文中使用的三維數(shù)值風(fēng)場模型在配合大灣區(qū)的地表粗糙長度能較好地模擬臺(tái)風(fēng)風(fēng)場。

圖10 臺(tái)風(fēng)“黑格比”在不同觀測(cè)站實(shí)測(cè)與模擬的風(fēng)速和風(fēng)向?qū)Ρ菷ig.10 Comparison of measured and simulated wind speed and direction of Typhoon Hagupit at different observatories

4 大灣區(qū)臺(tái)風(fēng)危險(xiǎn)性分析

為了研究粵港澳大灣區(qū)的臺(tái)風(fēng)危險(xiǎn)性，本研究以模擬的萬年臺(tái)風(fēng)事件集作為輸入，采用三維臺(tái)風(fēng)數(shù)值風(fēng)場模型模擬經(jīng)過大灣區(qū)范圍內(nèi)的臺(tái)風(fēng)風(fēng)場，以估算大灣區(qū)的臺(tái)風(fēng)危險(xiǎn)性。表4 展示了本研究估算的大灣區(qū)11個(gè)城市的50年和100年重現(xiàn)期臺(tái)風(fēng)風(fēng)速，并與已有的研究以及《建筑結(jié)構(gòu)荷載規(guī)范GB50009-2012》（下面簡稱“規(guī)范”）中的粵港澳大灣區(qū)城市的50年和100年重現(xiàn)期風(fēng)速進(jìn)行對(duì)比［9，28，43］。由于“規(guī)范”中只給出中國各城市的風(fēng)壓，所以本研究依據(jù)“規(guī)范”中的風(fēng)壓與風(fēng)速關(guān)系式（E.2.4-1），將基本風(fēng)壓轉(zhuǎn)換為風(fēng)速。此外，“規(guī)范”中基本風(fēng)壓表（表E.5）缺失部分大灣區(qū)城市，故缺失的風(fēng)速依據(jù)“規(guī)范”中的全國基本風(fēng)壓分布圖（圖E.6.3）進(jìn)行近似估算?！耙?guī)范”只給出50年重現(xiàn)期的風(fēng)速分布圖，并未給出100年重現(xiàn)期的風(fēng)速分布圖。

從表4 可以看出，本研究估算結(jié)果與“規(guī)范”及已有研究文獻(xiàn)的最大偏差出現(xiàn)在深圳的100 年重現(xiàn)期風(fēng)速，估算結(jié)果比“規(guī)范”和Chen 和Duan［28］值分別小3 m/s和1.7 m/s。這可能是由于本文采用的三維數(shù)值風(fēng)場模型考慮了深圳東部和香港北部的高山對(duì)臺(tái)風(fēng)風(fēng)場的快速削弱作用，這一點(diǎn)同樣也反映在粵港澳大灣區(qū)50年和100年重現(xiàn)期臺(tái)風(fēng)等風(fēng)速線上。從圖11（a）和（b）可以看出，海岸線附近風(fēng)速的變化梯度非常大，這說明三維數(shù)值風(fēng)場模型可以較好地捕捉到陸地對(duì)風(fēng)場的快速削弱作用。模擬結(jié)果顯示，珠江入海口兩岸的沿海區(qū)域的極值風(fēng)速都較大，但是深入內(nèi)陸數(shù)公里風(fēng)速就會(huì)迅速減弱。直接暴露在海洋上的海島（江門市上川島和香港橫瀾島）的極值風(fēng)速普遍都很大，50年重現(xiàn)期風(fēng)速均超過40 m/s。另外，市區(qū)靠海的城市的極值風(fēng)速也比較大，香港、澳門、珠海的50年重現(xiàn)期風(fēng)速均超過35 m/s。相反，對(duì)于市區(qū)深入內(nèi)陸的城市，其極值風(fēng)速則明顯偏小，比如肇慶市廣寧縣的50年重現(xiàn)期風(fēng)速僅為23.6 m/s。

表4 粵港澳大灣區(qū)11個(gè)城市的50年和100年重現(xiàn)期臺(tái)風(fēng)風(fēng)速V50和V100Table 4 Typhoon wind speeds for the 50-year and 100-year return periods for 11 cities in the Guangdong-Hong Kong-Macao Greater Bay Area

圖11 粵港澳大灣區(qū)50年和100年重現(xiàn)期臺(tái)風(fēng)風(fēng)速Fig.11 Typhoon wind speeds for the 50-year and 100-year return periods for the Guangdong-Hong Kong-Macao Greater Bay Area

5 結(jié)論

為了量化粵港澳大灣區(qū)的臺(tái)風(fēng)危險(xiǎn)性，本研究通過改進(jìn)和驗(yàn)證的統(tǒng)計(jì)動(dòng)力學(xué)全路徑合成臺(tái)風(fēng)模型生成萬年臺(tái)風(fēng)事件集，并對(duì)逐個(gè)臺(tái)風(fēng)事件進(jìn)行風(fēng)場模擬。首先，本研究使用CMA 最佳路徑數(shù)據(jù)集確定了西北太平洋海域的臺(tái)風(fēng)年平均生成數(shù)服從泊松分布，并采用核密度方法估算臺(tái)風(fēng)生成時(shí)空位置的概率。其次，本研究建立了臺(tái)風(fēng)移動(dòng)速度與引導(dǎo)氣流和β 漂流的回歸關(guān)系，并在移動(dòng)模型中考慮了臺(tái)風(fēng)移動(dòng)速度的隨機(jī)誤差。本研究基于自回歸模型建立了海上強(qiáng)度回歸模型，模型中考慮了垂直風(fēng)切、相對(duì)濕度、移動(dòng)速度對(duì)強(qiáng)度的影響，并基于自然指數(shù)遞減函數(shù)建立登陸大灣區(qū)臺(tái)風(fēng)的強(qiáng)度衰減模型。最后，本研究采用考慮了垂直平流和垂直擴(kuò)散過程的三維數(shù)值風(fēng)場模型模擬臺(tái)風(fēng)風(fēng)場。

本研究分別從西北太平洋海域和登陸大灣區(qū)臺(tái)風(fēng)的兩個(gè)角度對(duì)模擬臺(tái)風(fēng)事件集進(jìn)行驗(yàn)證。首先，對(duì)模擬和歷史臺(tái)風(fēng)年發(fā)生總數(shù)分布的非參數(shù)檢驗(yàn)表明二者分布一致。其次，西北太平洋海域2.5°×2.5°網(wǎng)格內(nèi)模擬與歷史臺(tái)風(fēng)的4 個(gè)特征參數(shù)的整體相關(guān)系數(shù)均能達(dá)到0.90 以上。對(duì)于登陸大灣區(qū)的臺(tái)風(fēng)，模擬和歷史年平均登陸個(gè)數(shù)接近，模擬和歷史登陸臺(tái)風(fēng)的移動(dòng)速度和強(qiáng)度的分布也較為吻合。大灣區(qū)沿海岸線各分段上的歷史年登陸數(shù)、移動(dòng)速度、方向以及登陸強(qiáng)度的平均值均落在模擬結(jié)果的90%區(qū)間內(nèi)?？偟膩碚f，本研究的改進(jìn)統(tǒng)計(jì)動(dòng)力學(xué)全路徑合成臺(tái)風(fēng)模型可以較好地模擬登陸大灣區(qū)的臺(tái)風(fēng)事件。此外，通過與大灣區(qū)兩個(gè)觀測(cè)站的臺(tái)風(fēng)“黑格比”實(shí)測(cè)風(fēng)速記錄對(duì)比，驗(yàn)證了本文采用的臺(tái)風(fēng)近地面三維風(fēng)場模型可以較好地模擬登陸大灣區(qū)臺(tái)風(fēng)的風(fēng)速。

最后本研究以模擬的萬年臺(tái)風(fēng)事件集作為輸入，與三維數(shù)值臺(tái)風(fēng)風(fēng)場模型相結(jié)合對(duì)大灣區(qū)進(jìn)行臺(tái)風(fēng)危險(xiǎn)性評(píng)估，并繪制了大灣區(qū)50年和100年重現(xiàn)期臺(tái)風(fēng)風(fēng)速圖。通過與已有研究以及“規(guī)范”中對(duì)大灣區(qū)城市的50年和100年重現(xiàn)期風(fēng)速記錄進(jìn)行對(duì)比，我們發(fā)現(xiàn)估算結(jié)果與“規(guī)范”及已有研究文獻(xiàn)的偏差不超過3 m/s，其中最大偏差出現(xiàn)在深圳的100年重現(xiàn)期風(fēng)速上。這是由于風(fēng)場模型考慮了深圳東部和香港北部的高山對(duì)臺(tái)風(fēng)風(fēng)場的快速削弱作用。另外，海岸線附近的風(fēng)速變化梯度非常大，珠江入?？趦砂兜难睾^(qū)域的極值風(fēng)速都較大，但是深入內(nèi)陸數(shù)公里風(fēng)速就會(huì)迅速減弱。