本發明屬于地震工程，尤其涉及一種適用于復雜地形的地震波場模擬譜單元建模方法。

背景技術：

1、在結構抗震設計中，合理的地震動參數是保證結構抗震安全性的前提。對于高壩大庫、核電設施等重要結構的抗震，需采用地震動時程作為輸入。但對于特定工程而言，滿足地震地質相似條件的真實記錄通常難以獲得，因此需采用人工模擬的方法加以補充。

2、目前，常用的地震動模擬方法包括隨機方法及基于物理的直接模擬方法。基于物理的直接模擬方法(如有限差分法、有限元法、譜元法等)能夠充分考慮震源、傳播路徑、場址相關的各種影響因素，從波動傳播的基本物理機制出發，具有天然的合理性。其中，譜元法結合了有限元方法的靈活性及偽譜法的精確性和快速收斂性，尤其擅長模擬地震動在不均勻介質及復雜地形條件下的傳播問題。

3、相較于隨機方法，基于物理的直接模擬方法往往計算量很大。在譜元法地震動模擬中，最大模擬頻率與介質波速成正比，與單元節點間的距離成反比，因此通常需要根據地殼波速結構進行淺部網格細化，以最大化計算效率。然而，在高山峽谷等地形條件復雜的地區，陡峭復雜的地形突變將劣化網格質量，導致模擬過程不收斂。

技術實現思路

1、針對現有技術中的上述不足，本發明提供的一種適用于復雜地形的地震波場模擬譜單元建模方法，能夠解決高山峽谷等地形條件復雜地區地殼結構建模時譜單元的質量問題，同時將計算自由度降到最低，實現基于譜元法的地震動穩定、低成本模擬。

2、為了達到以上目的，本發明采用的技術方案為：一種適用于復雜地形的地震波場模擬譜單元建模方法，包括以下步驟：

3、s1、確定研究區地形參數；

4、s2、根據研究區波速結構參數，確定網格細化倍數，并得到波速等值面距底部最小距離；

5、s3、基于確定的地形參數以及波速等值面距底部最小距離，確定地殼譜單元細化轉換層位置；

6、s4、根據地殼譜單元細化轉換層位置，確定網格細化前的三向單元總數；

7、s5、根據確定的研究區地形參數以及網格細化前的三向單元總數，采用六面體單元在全域進行網格剖分；

8、s6、根據轉換層上部的地殼譜單元層數，進行淺部地殼模型細化處理；

9、s7、根據細化處理結果，計算網格質量，統計地殼譜單元數量得到地殼譜單元預期值，并與估計值進行對比，完成對地殼譜單元的優化建模。

10、本發明的有益效果是：在高山峽谷等地形條件復雜的地區，陡峭復雜的地形突變將劣化譜單元網格質量，導致譜元法模擬過程不收斂，本發明提出了一種綜合考慮譜單元質量與譜元法計算效率的適用于復雜地形的地震波場模擬譜單元建模方法，通過計算轉換層上部和下部的尺寸，并采用六面體單元在全域進行網格剖分，以及計算地殼譜單元的數量，實現對地殼譜單元的優化建模處理，本發明能夠同時解決譜單元質量問題，同時將計算自由度降到最低，實現譜元法地震動的穩定、低成本模擬。

11、進一步地，所述步驟s1包括以下步驟：

12、s101、根據最大斷層深度確定最小模擬海拔，以及根據場址及斷層相對位置確定模擬長度l和模擬寬度w；

13、s102、對研究區地形結構進行建模，計算地形海拔最高點至最小模擬海拔的模擬高度h；

14、s103、計算研究區海拔最低點和海拔最高點之間的高程差δh，完成參數的確定。

15、再進一步地，所述步驟s2包括以下步驟：

16、s201、對地殼波速結構進行三維建模；

17、s202、基于建模結果，通過獲取場址地表剪切波速vs,0和區域最大波速vs,m，計算網格細化倍數；

18、s203、根據網格細化倍數，計算波速等值面距底部最小距離，其中，該最小距離為等波速面高度hr。

19、上述進一步方案的有益效果是：依據場址地表剪切波速vs,0和區域最大波速vs,m，計算網格細化倍數能夠確保全部網格滿足最大模擬頻率的要求。

20、再進一步地，所述波速等值面的表達式如下：

21、vsr＝nr·vs,0

22、

23、其中，vsr表示波速等值面，nr表示網格細化倍數，ceil()表示向上取整函數。

24、上述進一步方案的有益效果是：根據波速確定最小的網格細化倍數，能夠避免因網格過小而降低模擬效率。

25、再進一步地，所述步驟s3中地殼譜單元細化轉換層位置的表達式如下：

26、

27、其中，h2表示地殼譜單元細化轉換層深度，δh表示研究區海拔最低點和海拔最高點之間的高程差，cr和ch均表示常數，nr表示網格細化倍數，h表示地形海拔最高點至最小模擬海拔的模擬高度。

28、上述進一步方案的有益效果是：本發明通過上式，能夠確定使單元數量最少時的轉換層深度。

29、再進一步地，所述步驟s4包括以下步驟：

30、s401、利用下式計算得到地殼譜單元轉換層上部豎向尺寸dnear：

31、

32、其中，dnear(fm)表示擬模擬的最大頻率為fm時地殼譜單元轉換層上部豎向尺寸，fm表示擬模擬的最大頻率，ngll表示地殼譜單元的積分點數量，d0表示地殼譜單元的積分點數量為5、最大頻率為1hz以及地表剪切波速為1000m/s時八節點六面體單元的參考尺寸，vs,0表示場址地表剪切波速，nmsh表示地殼譜單元邊上的節點分段數量，八節點六面體單元時nmsh＝1，二十七節點六面體曲面單元時nmsh＝2；

33、s402、根據dnear(fm)，利用下式計算得到轉換層上部的地殼譜單元層數nz1：

34、

35、其中，h2表示地殼譜單元細化轉換層深度；

36、s403、根據地殼譜單元層數nz1，利用下式計算得到地殼譜單元轉換層上部的水平向尺寸a1；

37、

38、其中，h1表示區域地形最低點至轉換層的深度；

39、s404、根據地殼譜單元轉換層上部的水平向尺寸a1，利用下式分別計算得到轉換層上部的水平向單元數量和

40、

41、其中，l表示模擬長度，w表示模擬寬度；

42、s405、根據地殼譜單元轉換層上部的水平向尺寸a1，利用下式計算得到轉換層下部的三方向單元尺寸a2；

43、a2＝a1nr

44、其中，nr表示網格細化倍數；

45、s406、根據三方向單元尺寸a2，利用下式計算得到轉換層下部的豎向單元數量nz2：

46、

47、其中，h表示地形海拔最高點至最小模擬海拔的模擬高度；

48、s407、根據水平向單元數量nx1、ny1和豎向單元數量nz2，利用下式分別計算得到網格細化前的三向單元總數nx、ny和nz，完成對地殼譜單元尺寸的確定：

49、

50、其中，ceil()表示向上取整函數。

51、上述進一步方案的有益效果是：本發明通過確認網格細化前的三向單元總數，能夠方便地計算單元總數并估計模型總自由度。

52、再進一步地，所述步驟s7中地殼譜單元預期值的表達式如下：

53、

54、其中，n表示地殼譜單元預期值，c表示常數，h2表示地殼譜單元細化轉換層深度，δh表示研究區海拔最低點和海拔最高點之間的高程差，nr表示網格細化倍數，h表示地形海拔最高點至最小模擬海拔的模擬高度，l表示模擬長度，w表示模擬寬度，dnear(fm)表示表示擬模擬的最大頻率為fm時地殼譜單元轉換層上部豎向尺寸。

55、上述進一步方案的有益效果是：計算總單元個數能夠方便地同其他工程進行對比，同時能夠估計模型總自由度。基于總自由度數量及同其他工程的對比結果可以為地震動模擬分配合理的計算資源。