伶仃洋河口水域納潮特性分析

方神光, 陳文龍, 崔麗琴

(珠江水利委員會(huì) 珠江水利科學(xué)研究院, 廣東 廣州 510611)

伶仃洋河口水域納潮特性分析

方神光, 陳文龍, 崔麗琴

(珠江水利委員會(huì) 珠江水利科學(xué)研究院, 廣東 廣州 510611)

為闡明伶仃洋河口水域納潮現(xiàn)狀和特性, 分別采用經(jīng)驗(yàn)公式和建立貼體曲線坐標(biāo)系下的珠江河口二維數(shù)學(xué)模型的方法, 對(duì)河口伶仃洋水域的納潮特性進(jìn)行了分析和探討。數(shù)學(xué)模型的離散和求解采用純隱格式的混合有限分析法, 根據(jù)實(shí)測(cè)資料對(duì)其進(jìn)行了驗(yàn)證。比較顯示, 采用實(shí)測(cè)資料和經(jīng)驗(yàn)公式計(jì)算到的伶仃洋水域納潮量與數(shù)學(xué)模型的計(jì)算結(jié)果相差在5%以內(nèi), 從而證實(shí)了數(shù)學(xué)模型的可靠性。結(jié)果表明, 伶仃洋海域的主要納潮口門為伶仃洋口門, 約占總納潮量的 87.7%, 通過(guò)香港暗士頓水道進(jìn)出的潮量?jī)H占 12.3%左右, 顯示伶仃洋口門為該海域的主要徑潮通道, 對(duì)該河口的防洪和生態(tài)具有重要的作用。

伶仃洋; 納潮量; 經(jīng)驗(yàn)公式; 數(shù)學(xué)模型

隨著近 30 a來(lái)填海造陸等工程建設(shè), 伶仃洋海域河口岸線及地形條件變化顯著, 對(duì)該水域的納潮量影響巨大。納潮量是一個(gè)海灣可以接納潮水的體積, 是海灣環(huán)境評(píng)價(jià)的重要指標(biāo), 對(duì)伶仃洋海域的納潮演變進(jìn)行分析和研究, 有助于理解伶仃洋海域近些年來(lái)的生態(tài)環(huán)境變化, 并對(duì)開展的河口治理提供科學(xué)依據(jù)。有關(guān)河口納潮量的研究和分析成果較多, 如喬貫宇等[1]采用 ADCP(聲學(xué)多普勒海流剖面儀)方法對(duì)膠州灣口進(jìn)行走航觀測(cè), 得到了膠州灣的納潮量; 陳紅霞等[2]同時(shí)采用海圖和 ADCP方法計(jì)算膠州灣的納潮, 并進(jìn)行了比較和分析; 葉海桃等[3]將灘涂和水域分開, 計(jì)算了三沙灣的納潮量和水體交換時(shí)間; 吳隆業(yè)等[4]采用遙感影像資料得到了?？诟鄣乃蛎娣e進(jìn)而計(jì)算了納潮量的變化; 楊世倫等[5]考慮到河口潮灘被圍墾的情況, 改進(jìn)了傳統(tǒng)的納潮計(jì)算公式, 對(duì)膠州灣的納潮量進(jìn)行了重現(xiàn)計(jì)算和分析。

有關(guān)伶仃洋河口的納潮研究成果相對(duì)較少, 但與納潮變化密切相關(guān)的河口水動(dòng)力、水環(huán)境及水文地形變化的研究成果很多。如徐峰俊等[6]和陳水森等[7]根據(jù)以往研究資料分析了伶仃洋的水動(dòng)力環(huán)境、水域面積、灘槽面積和沖淤變化情況等, 并提出了相應(yīng)的治理策略。本文將分別通過(guò)建立伶仃洋河口二維數(shù)學(xué)模型及采用經(jīng)驗(yàn)公式的方法來(lái)計(jì)算伶仃洋河口水域的納潮量, 將兩種方法相互進(jìn)行對(duì)比驗(yàn)證。水文驗(yàn)證資料采用2007年8月份的大潮和中潮資料以及2001年2月小潮實(shí)測(cè)資料。在此基礎(chǔ)上, 對(duì)伶仃洋水域的納潮特性進(jìn)行分析和探討。

1 納潮量計(jì)算方法

通常情況下, 納潮量是指平均潮差條件下某一海灣可能接納的海水量, 其計(jì)算公式為:

式中P為平均潮差條件下的納潮量, HΔ為平均潮差, S為平均水域面積(即平均高潮位與平均低潮位面積之均值)。平均高潮位和平均低潮位之差值, 可根據(jù)驗(yàn)潮站的多年實(shí)測(cè)資料求得, 該計(jì)算方法簡(jiǎn)單,計(jì)算結(jié)果趨于一致。然而, S的計(jì)算較為復(fù)雜。地形圖上一般只標(biāo)出0 m線(理論最低潮位)和岸線(相當(dāng)于最大高潮線), 位于兩者之間的平均高潮線和平均低潮線需經(jīng)推算方能得到。

對(duì)伶仃洋的納潮量計(jì)算, 此處一方面采用公式(1)進(jìn)行計(jì)算; 另一方面基于建立的珠江河口二維數(shù)學(xué)模型, 通過(guò)在河口口門位置設(shè)置特征斷面來(lái)統(tǒng)計(jì)伶仃洋的納潮量。

2 數(shù)學(xué)模型及驗(yàn)證

2.1 數(shù)學(xué)模型

采用的正交曲線坐標(biāo)系下的二維水深平均水動(dòng)力數(shù)學(xué)模型如下:

連續(xù)方程:

動(dòng)量方程:

式中: ξ, η 為曲線坐標(biāo)系(無(wú)量綱); U, V 分別為 ξ, η 方向的垂向平均流速(m/s); ξ為水位(基準(zhǔn)面到自由水面的距離, m); Q為源匯項(xiàng)(m3/s); H為總水深(m);,為曲線坐標(biāo)到直角坐標(biāo)的變換參數(shù)(m);ρ為水體密度(kg/m3); f為柯氏系數(shù)(s-1); g為重力加速度(m/s2); vt為紊動(dòng)黏性系數(shù)(m2/s); t為時(shí)間(s)。

2.2 數(shù)值計(jì)算方法及計(jì)算水域

二維潮流數(shù)學(xué)模型中的偏微分方程可以寫為如下統(tǒng)一形式:應(yīng)用本文提出的純隱格式的混合有限分析法對(duì)該方程進(jìn)行離散, 離散詳細(xì)步驟可以參見相關(guān)文獻(xiàn)[8]。同時(shí)此處采用C型網(wǎng)格結(jié)合SIMPLER算法用于速度和水位的耦合求解。由于采用曲線規(guī)則網(wǎng)格, 計(jì)算是沿行或列進(jìn)行, 因此計(jì)算前首先需要對(duì)計(jì)算區(qū)域的復(fù)雜邊界進(jìn)行識(shí)別, 其相關(guān)處理方法以及計(jì)算區(qū)域存在的淺灘的處理方式可以參見相關(guān)文獻(xiàn)[9]。

珠江口水域數(shù)模計(jì)算范圍如圖 1所示, 東西距離約63 km, 南北距離約145 km。整個(gè)計(jì)算域包括獅子洋、伶仃洋西東口門、香港水道、伶仃洋外萬(wàn)山群島等。采用貼體正交曲線網(wǎng)格, 共布網(wǎng)格436個(gè)×665個(gè), 最大網(wǎng)格尺寸250 m×140 m, 最小網(wǎng)格尺寸30 m×8 m。

2.3 初始條件和邊界條件

流場(chǎng)數(shù)值模擬的初始條件為靜止?fàn)顟B(tài); 上邊界四個(gè)河口(虎門、蕉門、洪奇門和橫門)給定為實(shí)測(cè)流量邊界, 外海開邊界不同時(shí)刻的潮位由實(shí)測(cè)值通過(guò)插值方法給定, 流速采用二類邊界條件, 即各時(shí)刻邊界上的流速值賦值為相鄰內(nèi)部網(wǎng)格點(diǎn)上的計(jì)算值;岸邊界給定為無(wú)滑移邊界條件。

2.4 露灘處理

計(jì)算區(qū)域存在淺灘, 由于淺灘在漲潮時(shí)淹沒,退潮時(shí)露出, 對(duì)淺灘的處理方式影響到計(jì)算的穩(wěn)定、收斂以及計(jì)算結(jié)果的精度。此處將 Leendertse等[10]在推出ADI數(shù)值模式的同時(shí)給出的干濕網(wǎng)格判別準(zhǔn)則與河道水流數(shù)值模擬中常用的凍結(jié)法結(jié)合起來(lái),準(zhǔn)確有效地處理計(jì)算區(qū)域的淺灘。其基本原理是:(1)某瞬時(shí)任一行(或列)計(jì)算完畢, 發(fā)現(xiàn)某網(wǎng)格點(diǎn)處總水深H小于設(shè)定的最小水深Hmin, 則認(rèn)為干出, 將該處糙率設(shè)為極大值(如 108), 即將該干節(jié)點(diǎn)處的流速凍結(jié), 同時(shí)該點(diǎn)水位為該點(diǎn)底部高程加上最小水深Hmin; (2)在任一時(shí)層計(jì)算前, 根據(jù)干網(wǎng)格節(jié)點(diǎn)周圍網(wǎng)格節(jié)點(diǎn)水位與該干網(wǎng)格節(jié)點(diǎn)水位比較來(lái)確定該干網(wǎng)格節(jié)點(diǎn)是否被淹沒, 若被淹沒, 則將其糙率恢復(fù)為正常值。

2.5 模型驗(yàn)證

分別選取大、中、小三個(gè)潮型對(duì)模型進(jìn)行了驗(yàn)證, 用于潮位驗(yàn)證的有 12個(gè)站點(diǎn), 用于流速和流向驗(yàn)證的有 15個(gè)站點(diǎn), 由于篇幅所限, 此處只給出大潮(2007年8月13日21時(shí)至8月14日22時(shí))的部分站點(diǎn)(4個(gè)潮位站和4個(gè)測(cè)流站, 如圖1所示)驗(yàn)證曲線, 如圖2、圖3和圖4所示, 其他可詳見文獻(xiàn)[11]。從模型驗(yàn)證結(jié)果來(lái)看, 無(wú)論是計(jì)算的量值還是位相,均與實(shí)測(cè)值基本吻合, 可真實(shí)復(fù)演伶仃洋的潮流運(yùn)動(dòng), 證明了本文所用模型的實(shí)用性。

3 納潮量的實(shí)測(cè)資料分析

圖1 計(jì)算區(qū)域及網(wǎng)格剖分圖Fig. 1 Calculation area and grids

圖2 潮位驗(yàn)證Fig. 2 Verifications for tidal levels

圖3 流速驗(yàn)證Fig. 3 Verifications for tidal velocities

圖4 流向驗(yàn)證Fig. 4 Verifications for flowing directions

采用公式(1)來(lái)計(jì)算伶仃洋納潮量時(shí), 選取處于伶仃洋正中間的內(nèi)伶仃島站、赤灣站和金星門站或淇澳站的實(shí)測(cè)潮位資料。統(tǒng)計(jì)顯示包括淺灘在內(nèi)的珠江口門以內(nèi)的伶仃洋和獅子洋的總水域面積約為1 900 km2, 由于存在淺灘, 潮位不同時(shí), 淺灘所占面積不同, 導(dǎo)致納潮水域面積也不一樣。表1列出了典型潮位站大、中、小潮型下實(shí)測(cè)高低潮位和潮差。根據(jù)建立的伶仃洋河口數(shù)學(xué)模型計(jì)算了不同高低潮位對(duì)應(yīng)的水域面積和平均納潮水域面積, 由此得到采用不同站點(diǎn)的潮差計(jì)算到的大潮、中潮和小潮的納潮量。

采用實(shí)測(cè)結(jié)果計(jì)算分析顯示, 大潮時(shí), 采用內(nèi)伶仃島站、赤灣站和金星門站的潮差計(jì)算到伶仃洋和獅子洋的總納潮量分別為 29.70億, 27.97億,29.46億m3, 平均后大潮期的納潮量為29.04億m3;中潮時(shí), 采用內(nèi)伶仃島站、赤灣站和金星門站的潮差計(jì)算到的納潮量分別為25.40億, 24.27億, 26.15億m3, 平均后中潮的納潮量為25.27億m3; 2001年2月小潮時(shí), 采用金星門和赤灣站計(jì)算到的納潮量分別為 18.83億和 20.02億 m3, 平均后小潮的納潮量為19.24億 m3; 顯示大潮納潮量比中潮納潮量多 13%,中潮納潮量比小潮多23.9%。計(jì)算顯示伶仃洋平均納潮水域面積是小潮最大, 中潮次之, 大潮最小, 主要是大潮的低潮位太低, 導(dǎo)致淺灘出露面積大, 相應(yīng)平均納潮水域面積減小幅度也大的緣故。

4 數(shù)值模擬計(jì)算納潮量分析

4.1 納潮斷面布置

根據(jù)建立的數(shù)學(xué)模型, 在珠江河口設(shè)置特征觀測(cè)斷面, 統(tǒng)計(jì)特征斷面上不同時(shí)刻的進(jìn)出潮量, 以此得到該水域在一個(gè)潮周期內(nèi)的納潮量。為此, 在珠江河口海域分別設(shè)置了11個(gè)特征潮量觀測(cè)斷面, 分別為伶仃洋口門斷面、暗士頓水道斷面、金星門斷面、內(nèi)伶仃斷面、赤灣斷面、深圳灣斷面、虎門斷面、青州斷面、九州斷面、西槽斷面和銅鼓斷面, 如圖 1所示。這些特征斷面的設(shè)置主要考慮兩方面的因素, 一方面要保證從數(shù)值模擬試驗(yàn)結(jié)果中準(zhǔn)確無(wú)誤地提取和統(tǒng)計(jì)該斷面潮量, 斷面的設(shè)置都與某一條網(wǎng)格線重合; 另一方面, 根據(jù)特征觀測(cè)斷面通過(guò)的潮量, 能掌握伶仃洋海域當(dāng)前的納潮量及其特性。

表1 納潮量實(shí)測(cè)資料計(jì)算結(jié)果Tab. 1 Calculation results and measured data of tide prisms

4.2 伶仃洋納潮特性分析

伶仃洋海域?qū)儆诘湫偷牟灰?guī)則半日潮, 本文根據(jù)建立的珠江河口數(shù)學(xué)模型對(duì)大、中、小潮型的數(shù)值模擬計(jì)算結(jié)果, 統(tǒng)計(jì)了 7個(gè)主要納潮觀測(cè)斷面在不同時(shí)刻的進(jìn)出潮量及一個(gè)潮周期內(nèi)的總納潮量,結(jié)果如圖 5、表 2、表 3和表 4所示。分析可見:(1)一個(gè)潮周期內(nèi), 伶仃洋海域在本次大、中、小潮通過(guò)珠江河口口門(伶仃洋口門和暗士頓口門的平均進(jìn)出潮量之和)的平均潮量分別是27.88億, 24.71億和 19.53億 m3, 與采用實(shí)測(cè)資料和公式(1)計(jì)算的納潮量(表1)比較來(lái)看, 兩者分別相差–0.77億, –0.56億和0.29億m3, 誤差在5%以內(nèi), 這一方面驗(yàn)證了采用數(shù)學(xué)模型和設(shè)置特征觀測(cè)斷面計(jì)算海域納潮量的準(zhǔn)確性和有效性, 另一方面體現(xiàn)了采用實(shí)測(cè)資料和公式計(jì)算海域納潮量在平均意義上的有效性; (2)在大、中、小潮期, 通過(guò)伶仃洋口門斷面和暗士頓水道斷面的平均納潮量分別占該海域總潮量的87.7%和12.3%;可見, 伶仃洋海域的主要納潮口門為伶仃洋口門,通過(guò)暗士頓水道的潮量?jī)H占伶仃洋海域總納潮量的較小部分; (3)伶仃洋海域?yàn)榈湫筒灰?guī)則半日潮, 大潮期, 通過(guò)伶仃洋口門斷面的潮量中, 第一次漲潮量(19.83億m3)大于第一次落潮量(6.82億m3), 第二次漲潮量(30.91億 m3)則明顯小于第二次落潮量(40.21億m3); 在中潮和小潮期間也呈現(xiàn)同樣的變化規(guī)律; 且潮型越大, 該現(xiàn)象越明顯; (4)淇澳島–內(nèi)伶仃島間的內(nèi)伶仃斷面是內(nèi)伶仃洋海域的主要納潮斷面, 由該斷面進(jìn)出的平均潮量約占內(nèi)伶仃海域總納潮量的 59.1%, 其次是赤灣斷面(內(nèi)伶仃島–赤灣),約占總納潮量的28.8%, 金星門斷面最小, 約占總納潮量的12.1%。大、中、小潮時(shí), 內(nèi)伶仃海域納潮占整個(gè)伶仃洋海域總納潮量比例分別是58.2%, 60.3%,62.4%, 內(nèi)伶仃海域納潮量的比例隨潮差減小呈增加趨勢(shì), 平均約為60.3%, 外伶仃海域納潮量占總納潮比例約 39.7%, 內(nèi)伶仃洋海域納潮量要大于外伶仃海域; (5)本次大、中、小潮型中, 通過(guò)虎門口門的納潮量占整個(gè)伶仃洋納潮量的平均比例為 5.3%, 且納潮比例隨潮差減小略呈增加趨勢(shì), 說(shuō)明此處模擬計(jì)算獅子洋的納潮量占總納潮量的比例為 5.3%; 深圳灣水域的平均納潮量分別為 0.57億, 0.44億和0.28億m3, 僅占總納潮量的1.8%; (6)西槽斷面和青州斷面進(jìn)出潮量占伶仃洋總納潮量的比例平均約為46.6%, 可見伶仃洋口門位置靠近西航道深水區(qū)是伶仃洋的主要納潮通道, 位于淺灘區(qū)的九州斷面僅占總納潮量的16.7%左右; (7)銅鼓斷面的進(jìn)出潮量占伶仃洋總納潮量的比例平均約 12.2%, 顯示伶仃洋西槽和東槽通過(guò)銅鼓斷面水域存在一定程度的潮量交換, 該斷面的漲潮量明顯大于落潮量, 說(shuō)明落潮時(shí)該部分潮量從暗士頓口門退出。

表2 珠江河口口門斷面納潮量計(jì)算成果表Tab. 2 Tide prisms in and out of cross sections at the Pearl River mouth by numerical simulations

表3 伶仃洋內(nèi)各特征口門納潮量及占總納潮量的比例Tab. 3 Proportions of tide prisms in and out of cross sections to the totality of the Lingdingyang Bay

表4 內(nèi)伶仃海域各口門納潮統(tǒng)計(jì)表Tab. 4 Tide prisms in and out of the cross sections in the inner Lingdingyang Sea

圖5 不同潮型下各斷面進(jìn)出潮量速率變化曲線Fig. 5 Tidal prism velocities in and out of cross sections in spring, middle and neap tides

總結(jié)以上分析可見, 采用實(shí)測(cè)資料和數(shù)值模擬計(jì)算結(jié)果統(tǒng)計(jì)到的伶仃洋海域納潮量結(jié)果基本一致,確保了潮量計(jì)算結(jié)果的精度和準(zhǔn)確性; 伶仃洋海域的主要納潮口門為伶仃洋口門(澳門–大濠島斷面),約占總納潮量的 87.7%, 通過(guò)香港暗士頓水道進(jìn)出的潮量?jī)H占 12.3%左右; 內(nèi)伶仃海域的納潮量大于外伶仃海域, 內(nèi)伶仃海域的納潮量約占整個(gè)伶仃洋海域納潮量的 60.0%, 外伶仃海域納潮約占總納潮量的40.0%。

5 結(jié)論

分別采用經(jīng)驗(yàn)公式和數(shù)學(xué)模型方法, 對(duì)珠江河口伶仃洋水域的納潮特性進(jìn)行了計(jì)算和分析, 主要結(jié)論如下: (1)建立了貼體曲線坐標(biāo)系下的二維珠江河口水流數(shù)學(xué)模型, 并采用純隱格式的混合有限分析法對(duì)數(shù)學(xué)模型進(jìn)行了離散和求解; 通過(guò)將數(shù)值模擬結(jié)果與實(shí)測(cè)站點(diǎn)的潮位和流速及流向比較證實(shí)了本數(shù)學(xué)模型和技術(shù)方法的有效性和可靠性; (2)采用實(shí)測(cè)資料和經(jīng)驗(yàn)公式計(jì)算和分析了伶仃洋海域的納潮量, 顯示本次所選大、中、小潮型下, 伶仃洋海域的納潮量分別為29.04億, 25.27億和19.24億m3; 與數(shù)值模擬計(jì)算結(jié)果相差–0.77億, –0.56億和0.29億m3, 誤差在5%以內(nèi); (3)伶仃洋海域的主要納潮口門為伶仃洋口門(澳門–大濠島斷面), 約占總納潮量的 87.7%, 通過(guò)香港暗士頓水道進(jìn)出的潮量?jī)H占12.3%左右; (4)內(nèi)伶仃海域的納潮量大于外伶仃海域, 內(nèi)伶仃海域的納潮量約占整個(gè)伶仃洋海域納潮量的 60.0%, 外伶仃海域納潮約占總納潮量的40.0%。

[1] 喬貫宇, 華鋒, 范斌, 等. 基于 ADCP灣口測(cè)流的納潮量計(jì)算[J]. 海洋科學(xué)進(jìn)展, 2008, 26(3): 285-291.

[2] 陳紅霞, 華鋒, 劉娜, 等. 不同方式的納潮量計(jì)算比較——以膠州灣 2006年秋季小潮為例[J]. 海洋科學(xué)進(jìn)展, 2009, 27(1): 11-15.

[3] 葉海桃, 王義剛, 曹兵. 三沙灣納潮量及灣內(nèi)外的水交換[J]. 河海大學(xué)學(xué)報(bào), 2007, 35(1): 96-98.

[4] 吳隆業(yè), 孫玉星, 王振先. 海口港納潮量及其變化遙感研究[J]. 海岸工程, 1997, 16(2): 1-5.

[5] 楊世倫, 陳啟明, 朱駿, 等.半封閉海灣潮間帶部分圍墾后納潮量計(jì)算的商榷——以膠州灣為例[J]. 海洋科學(xué), 2003, 27(8):43-47.

[6] 徐峰俊, 朱士康, 王華. 伶仃洋水動(dòng)力環(huán)境分析及治理策略探討[J]. 人民珠江, 2004, (1): 11-14.

[7] 陳水森, 黎夏, 鄒春洋, 等. 20年來(lái)珠江口伶仃洋灘槽變化及演變分析[J]. 海洋科學(xué), 2001, 25(6): 52-54.

[8] 李煒. 黏性流體的混合有限分析解法[M].北京: 科學(xué)出版社, 2000.

[9] 方神光, 陳純, 劉濤. 一種純隱格式在近海潮流數(shù)值模擬中的應(yīng)用[J]. 水道港口, 2008, 29(5): 305-309.

[10] Leendertse J J. A water quality simulation model for well-mixed estuaries and coastal seas, Vol.1, principle of computation[M]. New York: The Rand Corporation,1970.

[11] 方神光. 水流、水質(zhì)數(shù)學(xué)模型軟件研究和開發(fā)——以港珠澳大橋?qū)χ榻谒虻募{潮和水體交換的影響研究為例[R]. 廣州: 珠江水利科學(xué)研究院, 2010.

Characteristics of tidal prism of the Lingdingyang Bay

FANG Shen-guang, CHEN Wen-long, CUI Li-qin
(Pearl River Hydraulic Research Institute of Pearl River Water Resource Committee, Guangzhou 510611, China)

Sep., 20, 2010

the Lingdingyang Bay; tidal prism; empirical formula; mathematical model

Tidal prism characteristics of the Lingdingyang Bay were analyzed on the basis of both empirical formula and a two-dimensional mathematical model under body-fitted curvilinear coordinates. The pure implicit scheme of hybrid finite analytic method was applied to discrete and solute equations. The mathematical model was validated with the measured data. The error between tidal prism calculated by empirical formula using measured data and tidal prism of the Lingdingyang Bay calculated by the mathematical model was less than 5%. The results indicated that the Lingdingyang Bay mouth was the main path of tidal current which accounted for about 87.7% of total tidal prism and Anshidun waterway accounted for about 12.3%. The Lingdingyang Bay mouth was the main channel of flood and tide for the sea and played important role on flood control and estuarine ecosystem.

TV148+.2

A

1000-3096(2012)06-0105-08

2010-09-20;

2012-02-10

國(guó)家自然科學(xué)基金項(xiàng)目(51109232); 廣東省自然科學(xué)基金項(xiàng)目(10151061101000001)

方神光(1978-), 男, 湖北監(jiān)利人, 副研究員, 博士, 從事環(huán)境水動(dòng)力學(xué)研究, E-mail: fangshenguang@163.com

(本文編輯：劉珊珊)