本發(fā)明涉及地震波數值模擬研究技術領域，特別涉及一種全局優(yōu)化的交錯網格有限差分正演模擬方法和裝置。

背景技術：

地震波正演模擬即模擬地震波在地球介質中的傳播過程，并研究地震波的傳播特征與地球介質參數的關系，達到與實際觀測地震記錄的最優(yōu)逼近。由于實際地球介質的異常復雜，很難得到地震波波動方程的解析解。故地震波數值模擬方法是人們正確認識地震波的傳播規(guī)律，驗證各種假設地球模型正確性的主要方法技術；是地震勘探和地震學的重要基礎，對實際地震資料的地質解釋和儲層預測等工作具有重要的理論指導作用和實際利用價值。由于廣泛地被應用在地震勘探和天然地震領域中，地震數值模擬方法也得到迅速發(fā)展和有效應用。目前，地震波正演模擬的數值方法主要有：射線追蹤法、積分方程法和波動方程法三大類。波動方程法是建立在地震波傳播方程的理論基礎上，對復雜介質中地震波傳播具有廣泛適應性，在地震波數值模擬中應用最廣泛。

有限差分法是偏微分方程的主要數值解法之一，是最早用于地震波數值模擬中的數值模擬方法。波動方程有限差分法相比于其他數值模擬方法如偽譜法、有限元法，因為其存儲的有效利用、計算量低且簡單易于實現而更受人偏愛廣泛應用在地震數值模擬、偏移和反演中。1968年Alterman和Karal最先嘗試將有限差分法應用于層狀介質彈性波傳播的數值模擬中，產生了理論地震圖。1972年，Boore將有限差分法用于非均勻介質地震波傳播的模擬。隨后，有限差分法被進一步用粘彈性、各向異性等復雜介質中和三維模型中地震波傳播的模擬。

常規(guī)的有限差分法正演時，會產生嚴重的數值假頻，稱為數值頻散，無法使波傳播較大距離。而且時間域算法的誤差會逐漸積累，導致計算精度不足，數值模擬結果分辨率降低。這是因為對模型進行離散成網格，并用差分算子逼近微分算子，從而使波動方程的系數發(fā)生變化即使相速度變成離散網格間隔的函數。因此當每個波長內采樣點太少時，就會產生網格數值頻散；如果增加單位波長內采樣點，又會增加計算量和存儲量，使正演成本增加。

如何消除這種數值頻散，提高模擬精度是應用有限差分法正演的關鍵問題之一。為了達到這一目標，前人從各方面作了很多努力：在時間-空間域有限差分方面，從Alford等在1974年的二階有限差分法到Dablain于1986年提出了高階有限差分法算子，隨著精度階數增加，單位波長需要的采樣點越來越少。一般情況下，時間方向采用二階差分格式，空間方向采用高階差分方式。從規(guī)則網格發(fā)展到交錯網格(Virieux在1986首先提出；Ozdenvar&Mcmechan，1997；董良國，2000；KetilHokstad，2003)，交錯網格方法提高了地震模擬的精確程度及計算穩(wěn)定性，并消除部分假像。此外，Boris和Book1973年、1975提出用求解守恒方程的通量校正傳輸方法(FCT)，有效的壓制了在粗網格上用差分算子引起的數值頻散。

由于用高階有限差分算子計算每個格點值時，參與計算的點數增加，使高階有限差分的計算量大大增加。而且有人提出高階有限差分算子存在飽和效應，且計算不穩(wěn)定。極限有限差分法即偽譜法，不存在數值頻散問題且每個波長只需要兩個采樣點，但計算量太大。

Zhou和Greenhalgh(1992)使用了廣義的漢寧(Hanning)窗截斷偽譜法算子得到了一種優(yōu)化的有限差分算子；Igel et al.(1995)使用了高斯窗截斷得到了一種優(yōu)化的有限差分算子。優(yōu)化的有限差分算子相對于常規(guī)有限差分算子其譜覆蓋范圍增大，即利用低階的有限差分算子可以達到常規(guī)高階算子的精度。Chu and stoffa2012利用二項式窗函數統(tǒng)一了有限差分法和偽譜法：用二項式窗函數對偽譜法算子進行截斷，得到常規(guī)有限差分系數；并改進了二項式窗函數得到了優(yōu)化的有限差分系數，有效提高了差分精度。wang等2014、2015年利用組合的窗函數直觀可視化地控制其主瓣和旁瓣的形狀，進而調節(jié)有限差分算子逼近微分算子的精度。Sun在2015年通過對一階導數的有限寬頻譜進行傅里葉(Fourier)反變換，得出了與chu在2012文中一致的交錯網格差分算子形式，并用正規(guī)hanning窗函數截斷后得到優(yōu)化的差分系數。

Holberg于1987用最優(yōu)化方法最小化群速度的相對峰值誤差，而Etgen于2007提出用相速度構建優(yōu)化問題的目標函數，由此均得到了優(yōu)化的差分系數。Liu等2009提出同時考慮時間導數和空間導數，基于時空域頻散關系，用泰勒展開法求解地震波方程的有限差分算子；并在2013、2014用最小二乘理論確定了的時空域優(yōu)化的空間有限差分算子。Zhang等于2013年提出了在給定目前最小的誤差限下(萬分之一)，用模擬退火算法直接對近似空間二階導數的有限差分算子進行優(yōu)化，使其覆蓋更大的波數范圍。Yang、Yan等于2014、2015、2016年分別用最小二乘法、樣點逼近法及雷米茲(Remez)算法確定了優(yōu)化的空間交錯網格有限差分算子。Liang等在2013、2014用線性化方法確定了時空域優(yōu)化的空間有限差分算子。辛維等2015分別用線性化算子和模擬退火法確定了給定波數上限下優(yōu)化的空間交錯網格有限差分算子，并進行了對比。以上這些方法，都一定程度上減小了有限差分模擬的數值頻散，并提高了數值模擬的計算效率。

綜上所述，對有限差分算子進行進一步的優(yōu)化，可分為窗函數優(yōu)化算法和直接優(yōu)化算法。窗函數優(yōu)化方法是選擇合適的窗函數去截斷偽普法算子，從而得到優(yōu)化的有限差分算子。窗函數優(yōu)化法的一個關鍵問題就是窗函數的選擇，其直接影響有限差分算子精度。該方法雖然優(yōu)化過程直觀，但是無法直接控制有限差分算子的精度，且受窗函數限制也無法達到很高的精度。直接優(yōu)化算法是則是利用最小二乘法、Remez算法等最優(yōu)化方法直接求使差分算子在波數域的誤差曲線盡可能滿足需求的差分系數。該類方法使有限差分算子的設計簡化成一優(yōu)化問題，給定譜覆蓋范圍，用優(yōu)化方法搜索使誤差盡可能小的有限差分系數。

針對現有技術的實現方案，Yang在2014年提出的基于最小二乘法理論的交錯網格有限差分算子優(yōu)化方案中，是用最小二乘法對給定波數范圍內的誤差的積分求極小，從而確定一種優(yōu)化的差分系數。Yang在2015年提出的線性化方法求解優(yōu)化交錯網格有限差分算子方案中，也需要給定波數覆蓋范圍，并進行等間隔采樣；帶入誤差函數后形成一線性方程組，然后用高斯消元法求解得唯一解。辛維于2015年應用最大最小原則作為模擬退火法優(yōu)化算法求解交錯網格有限差分算子的數值頻散誤差判定標準，得到了特定波數覆蓋范圍下優(yōu)化的差分系數。

然而，上述兩種方案以及Yang在2016年提出的Remez算法均存在以下缺陷：1)均需要預先設定優(yōu)化差分算子所能覆蓋的最大波數范圍；2)均不能直接控制頻散誤差，需要先畫出相對誤差曲線才能進行誤差控制；3)僅對相對誤差進行了分析，以相對誤差大小判斷差分算子精度不夠精確。因此，靈活性差、應用局限。

技術實現要素：

本發(fā)明實施例提供了一種全局優(yōu)化的交錯網格有限差分正演模擬方法，以解決現有技術中優(yōu)化交錯網格有限差分算子的方案中存在的靈活性差、應用局限的技術問題。該方法包括：獲得交錯網格有限差分算子；基于最大化范數建立目標函數，求解所述目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化所述交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，其中，所述交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的；利用所述優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。

在一個實施例中，所述目標函數為：

其中，E(c)為所述目標函數的極值；是優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍；c_n是第n個所述交錯網格有限差分權系數；k_x是波數覆蓋范圍；T是優(yōu)化的交錯網格有限差分算子的最大譜誤差容限；_N是得到所述交錯網格有限差分算子所用的網格點數的長度；Δx是采樣間隔。

在一個實施例中，求解所述目標函數，包括：采用模擬退火算法求解所述目標函數。

在一個實施例中，采用模擬退火算法求解所述目標函數，包括：設定所述優(yōu)化的交錯網格有限差分算子的最大譜誤差容限為預設譜誤差容限數值；通過循環(huán)以下步驟求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的相應交錯網格有限差分權系數：將波數在0到尼奎斯特波數的范圍內離散采樣，按照波數采樣點由小到大的順序，先將最小的一個波數采樣點作為預設波數范圍的上限，該預設波數范圍的下限為0，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在該得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在得到的交錯網格有限差分權系數使得所述目標函數的極值大于所述預設譜誤差容限數值時，終止計算，將上一個波數采樣點采用模擬退火算法計算得到的交錯網格有限差分權系數輸出，作為所述優(yōu)化的交錯網格有限差分算子的交錯網格有限差分權系數；在計算得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，繼續(xù)再將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數。

在一個實施例中，在求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的每個交錯網格有限差分權系數之前，還包括：設定模擬退火算法的自變量維度和各維度自變量的取值范圍，其中，自變量維度表示所述交錯網格有限差分權系數的個數，各維度自變量的取值范圍表示各個所述交錯網格有限差分權系數的振幅絕對值范圍，將自變量維度設定為得到所述交錯網格有限差分算子所用的網格點數的長度的二分之一；將各維度自變量的取值范圍設定為區(qū)間【0，2】；設定各個所述交錯網格有限差分權系數的振幅圍繞中心位置呈衰減的震蕩。

本發(fā)明實施例還提供了一種全局優(yōu)化的交錯網格有限差分正演模擬裝置，以解決現有技術中優(yōu)化交錯網格有限差分算子的方案中存在的靈活性差、應用局限的技術問題。該裝置包括：有限差分算子獲取模塊，用于獲得交錯網格有限差分算子；目標函數求解模塊，用于基于最大化范數建立目標函數，求解所述目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化所述交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，其中，所述交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的；地震波正演模塊，用于利用所述優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。

在一個實施例中，所述目標函數為：

其中，E(c)為所述目標函數的極值；是優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍；c_n是第n個所述交錯網格有限差分權系數；k_x是波數覆蓋范圍；T是優(yōu)化的交錯網格有限差分算子的最大譜誤差容限；_N是得到所述交錯網格有限差分算子所用的網格點數的長度；Δx是采樣間隔。

在一個實施例中，所述目標函數求解模塊具體用于采用模擬退火算法求解所述目標函數。

在一個實施例中，所述目標函數求解模塊，包括：誤差設定單元，用于設定所述優(yōu)化的交錯網格有限差分算子的最大譜誤差容限為預設譜誤差容限數值；目標函數求解單元，用于通過循環(huán)以下步驟求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的相應交錯網格有限差分權系數：將波數在0到尼奎斯特波數的范圍內離散采樣，按照波數采樣點由小到大的順序，先將最小的一個波數采樣點作為預設波數范圍的上限，該預設波數范圍的下限為0，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在該得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在得到的交錯網格有限差分權系數使得所述目標函數的極值大于所述預設譜誤差容限數值時，終止計算，將上一個波數采樣點采用模擬退火算法計算得到的交錯網格有限差分權系數輸出，作為所述優(yōu)化的交錯網格有限差分算子的交錯網格有限差分權系數；在得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，繼續(xù)再將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數。

在一個實施例中，所述目標函數求解模塊，還包括：參數設定單元，用于在求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的每個間交錯網格有限差分權系數之前，設定模擬退火算法的自變量維度和各維度自變量的取值范圍，其中，自變量維度表示所述交錯網格有限差分權系數的個數，各維度自變量的取值范圍表示各個所述交錯網格有限差分權系數的振幅絕對值范圍，將自變量維度設定為得到所述交錯網格有限差分算子所用的網格點數的長度的二分之一；將各維度自變量的取值范圍設定為區(qū)間【0，2】；設定各個所述交錯網格有限差分權系數的振幅圍繞中心位置呈衰減的震蕩。

在本發(fā)明實施例中，通過基于最大化范數建立目標函數，并求解目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，進而利用優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。由于目標函數是基于最大化范數建立的，相對于現有技術中最小平方法建立的目標函數，使得優(yōu)化交錯網格有限差分算子的過程更加靈活；同時，由于通過求解目標函數得到的交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的，有利于更加靈活地控制優(yōu)化的交錯網格有限差分算子的絕對譜誤差，使得優(yōu)化的交錯網格有限差分算子在較緊的譜誤差容限內還能夠覆蓋盡可能大的波數范圍，從而更有效地壓制數值頻散，提高地震模擬精度，有利于為逆時偏移、全波形反演提供高精度、高效的正演模擬方案，降低應用局限的問題。

附圖說明

此處所說明的附圖用來提供對本發(fā)明的進一步理解，構成本申請的一部分，并不構成對本發(fā)明的限定。在附圖中：

圖1是本發(fā)明實施例提供的一種全局優(yōu)化的交錯網格有限差分正演模擬方法的流程圖；

圖2是本發(fā)明實施例提供的一種不同譜誤差容限T下的8階交錯網格有限差分權系數絕對誤差對比示意圖；

圖3是本發(fā)明實施例提供的一種圖2的放大示意圖；

圖4是本發(fā)明實施例提供的一種采用模擬退火算法求解目標函數的流程圖；

圖5(a)是本發(fā)明實施例提供的一種譜誤差容限T＝0.0001時的高階優(yōu)化的交錯網格差分算子絕對誤差對比示意圖；

圖5(b)是本發(fā)明實施例提供的一種圖5(a)的放大示意圖；

圖6(a)是本發(fā)明實施例提供的一種700ms時刻8階常規(guī)交錯網格有限差分算子模擬X分量的示意圖；

圖6(b)是本發(fā)明實施例提供的一種700ms時刻8階優(yōu)化的交錯網格有限差分算子模擬X分量的示意圖；

圖6(c)是本發(fā)明實施例提供的一種700ms時刻8階常規(guī)交錯網格有限差分算子模擬Z分量的示意圖；

圖6(d)是本發(fā)明實施例提供的一種700ms時刻8階優(yōu)化的交錯網格有限差分算子模擬Z分量的示意圖；

圖7是本發(fā)明實施例提供的一種縱波速度模型示意圖；

圖8(a)是本發(fā)明實施例提供的一種marmousi2模型彈性波正演的垂直分量波場記錄的示意圖；

圖8(b)是本發(fā)明實施例提供的一種圖8(a)中區(qū)域1的放大示意圖；

圖8(c)是本發(fā)明實施例提供的一種圖8(a)中區(qū)域2的放大示意圖；

圖9是本發(fā)明實施例提供的一種全局優(yōu)化的交錯網格有限差分正演模擬裝置的結構框圖。

具體實施方式

為使本發(fā)明的目的、技術方案和優(yōu)點更加清楚明白，下面結合實施方式和附圖，對本發(fā)明做進一步詳細說明。在此，本發(fā)明的示意性實施方式及其說明用于解釋本發(fā)明，但并不作為對本發(fā)明的限定。

從現有技術中，我們發(fā)現用交錯網格有限差分數值進行數值模擬時，由于網格離散導致的數值頻散仍不可避免。關鍵是在保證精度的情況下，如何增大譜覆蓋范圍，以減少單位波長內采樣點數，從而減少存儲量和計算量；在單位波長采樣點固定時，則能更好地壓制高波數數值頻散以適應高頻地震模擬?；谔├照归_的常規(guī)交錯有限差分算子所覆蓋的譜范圍最小，但是波數較小時能最精確的逼近微分算子?；谧钚《朔▋?yōu)化的差分算子考慮奈奎斯特(Nyquist)波數內所有波數下的誤差積分極小，結果所覆蓋的譜范圍最大；但是在波數較小時與微分算子的譜誤差較大，無法有效壓制數值頻散。因此，壓制交錯網格有限差分算子的關鍵是在誤差范圍和譜覆蓋范圍之間找到一種平衡，本申請發(fā)明人提出上述全局優(yōu)化的交錯網格有限差分正演模擬方法，利用最大化范數建立的目標函數，則可以靈活地控制譜誤差容限，使給定波數范圍內所有樣點的譜誤差均在該譜誤差容限內。因此，首先找到合適的優(yōu)化算法求解該目標函數，然后選擇誤差容限，使得到的優(yōu)化的交錯網格有限差分算子在較緊的誤差限內還能夠覆蓋盡可能大的波數范圍。從而更有效地壓制數值頻散，為逆時偏移、全波形反演提供高精度、高效的正演模擬方案。

在本發(fā)明實施例中，提供了一種全局優(yōu)化的交錯網格有限差分正演模擬方法，如圖1所示，該方法包括：

步驟101：獲得交錯網格有限差分算子；

步驟102：基于最大化范數建立目標函數，求解所述目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化所述交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，其中，所述交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的；

步驟103：利用所述優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。

由圖1所示的流程可知，在本發(fā)明實施例中，通過基于最大化范數建立目標函數，并求解目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，進而利用優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。由于目標函數是基于最大化范數建立的，相對于現有技術中最小平方法建立的目標函數，使得優(yōu)化交錯網格有限差分算子的過程更加靈活；同時，由于通過求解目標函數得到的交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的，有利于更加靈活地控制優(yōu)化的交錯網格有限差分算子的絕對譜誤差，使得優(yōu)化的交錯網格有限差分算子在較緊的譜誤差容限內還能夠覆蓋盡可能大的波數范圍，從而更有效地壓制數值頻散，提高地震模擬精度，有利于為逆時偏移、全波形反演提供高精度、高效的正演模擬方案，降低應用局限的問題。

具體實施時，在進行地震波正演模擬的過程中，利用優(yōu)化的交錯網格有限差分算子近似計算空間導數，時間導數則可以用常規(guī)二階交錯網格有限差分算子近似計算。

具體實施時，本申請對獲得上述交錯網格有限差分算子(該交錯網格有限差分算子可以為高階交錯網格有限差分算子)的方式不做限定，用sinc函數插值理論可以推導出高階交錯網格有限差分算子，當然用連續(xù)傅里葉反變換也可以推導出的一致的結果。例如，我們用sinc函數插值理論通過以下步驟獲得高階交錯網格有限差分算子：

根據離散信號的采樣理論，一個帶限的連續(xù)信號f(x)可以被以一個均勻采樣的信號f_n通過sinc函數插值重建：

其中，Δx為采樣間隔，為截止波數。

我們?yōu)榈玫胶瘮礷(x)的一階導數，可以對公式(1)左右兩邊分別微分，即：

對于一階導數，交錯網格有限差分算子相對于常規(guī)有限差分算子具有更高的精度和更佳的穩(wěn)定性，通常在波場模擬中獲得更好的結果。因此，為推導出一階導數的高階交錯網格有限差分算子，我們選擇用公式(2)在處來估計一階導數如下：

因此，用一個長度為N點的窗函數(其中N為偶數)，去截斷公式(3)，可得到N點空間交錯網格有限差分算子：

對于常規(guī)交錯網格有限差分算子，為二項式窗函數(Chu和Stoffa，2012)，表示如下：

式中，l≤m；

是二項式系數公式，l≤m。也就是說，常規(guī)交錯網格有限差分算子權系數對應于微分的sinc插值系數乘以加權的的二項式窗函數。對于一階導數的交錯網格有限差分權系數c_n可以用加權的二項式窗函數找到，如下：

由公式(8)注意到，權系數因此公式(7)可以寫成：

至此，我們基于sinc函數插值理論推導出了N點常規(guī)交錯網格有限差分算子。

對公式(9)進行Fourier變換可得：

式中，k_x為波數覆蓋范圍。此外，將公式(10)右邊的正弦函數用泰勒級數展開，可建立N階線性方程組，通過求解該方程組也能推導出交錯網格常規(guī)有限差分算子。正因為常規(guī)交錯網格有限差分算子無論是高階還是低階算子均是圍繞零波數進行泰勒展開逼近，使得其在小波數范圍內模擬精度很高，而在大波數范圍內由于誤差較大而導致寬頻地震模擬中嚴重的空間頻散。

具體實施時，為適應寬頻帶地震模擬需要，我們必須拓寬交錯網格有限差分算子的波數覆蓋范圍，提高高頻地震模擬精度，因此，在本實施例中，利用最大化范數控制一定波數范圍內的最大絕對譜誤差，建立上述目標函數如下：

其中，E(c)為所述目標函數的極值；是優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍；c_n是第n個所述交錯網格有限差分權系數；k_x是波數覆蓋范圍；T是優(yōu)化的交錯網格有限差分算子的最大譜誤差容限；_N是得到所述交錯網格有限差分算子所用的網格點數的長度；Δx是采樣間隔。使得我們可以得到譜覆蓋范圍更大的優(yōu)化的交錯網格有限差分權系數。很明顯，基于最大化范數建立的目標函數比最小平方法建立的目標函數更加靈活。

具體實施時，在建立目標函數之后，可知譜誤差容限T決定了交錯網格有限差分算子能否尋優(yōu)成功以及其最終逼近精度。圖2為譜誤差容限T分別設為0.00005、0.0001、0.0002、0.0004和0.0005時，相應得到的8階交錯網格有限差分權系數的絕對誤差譜曲線，并與常規(guī)交錯網格差分系數(圖2中實線曲線表示常規(guī)交錯網格差分系數，實線曲線上的4、6、8、10等數字表示常規(guī)交錯網格差分系數的階數)進行對比。從圖2中可以發(fā)現隨著譜誤差容限T逐漸變松，交錯網格有限差分權系數的譜覆蓋范圍逐漸增大，T＝0.0005時，交錯網格有限差分權系數最大且超過了16階常規(guī)交錯網格有限差分系數；從圖3中則發(fā)現，隨著譜誤差容限T逐漸變松，交錯網格有限差分權系數的譜誤差在小波數時波動逐漸變大。因此，我們可以得出結論：設定譜誤差容限T若太小，則優(yōu)化的交錯網格有限差分算子覆蓋的譜范圍過小使寬頻帶彈性波模擬精度偏低；若設定譜誤差容限T太大，則優(yōu)化的交錯網格有限差分算子覆蓋譜范圍較大但是誤差積累影響比較嚴重，無法適應長時間步、深層彈性波模擬。

模擬退火法(SA)于1953年Metropolis提出屬于一種隨機算法，并由Kirkpatrick在1983年將其應用到組合優(yōu)化問題。模擬退火算法采用隨機采樣方式尋找新解，并以隨溫度下降而逐漸變小的概率接受較差的解。因此，模擬退火算法是局部搜索算法的擴展，具有能夠跳出局部最優(yōu)陷阱的優(yōu)點，較快的收斂到全局最優(yōu)解。用模擬退火法尋優(yōu)時初始溫度和終止溫度的選擇非常重要，直接決定搜索多長，以及能否成功搜索到滿足誤差上限的解。因此，針對如此復雜的目標函數，線性搜索方法和最小二乘法無法靈活求解。本申請?zhí)岢霾捎媚M退火法對上述目標函數進行求解，直接搜索最優(yōu)的交錯網格有限差分權系數。

具體實施時，采用模擬退火算法求解所述目標函數，包括：設定所述優(yōu)化的交錯網格有限差分算子的最大譜誤差容限為預設譜誤差容限數值；通過循環(huán)以下步驟求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的相應交錯網格有限差分權系數：將波數在0到尼奎斯特波數的范圍內離散采樣，按照波數采樣點由小到大的順序，先將最小的一個波數采樣點作為預設波數范圍的上限，該預設波數范圍的下限為0，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在該得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在得到的交錯網格有限差分權系數使得所述目標函數的極值大于所述預設譜誤差容限數值時，終止計算，將上一個波數采樣點采用模擬退火算法計算得到的交錯網格有限差分權系數輸出，作為所述優(yōu)化的交錯網格有限差分算子的交錯網格有限差分權系數；在計算得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，繼續(xù)再將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數。

具體的，本申請采用模擬退火算法求解上述目標函數，并不直接設定優(yōu)化的交錯網格有限差分算子所能覆蓋的最大波數范圍，而是在0到尼奎斯特波數的范圍內離散采樣后，在波數采樣點為最大波數上限的預設波數范圍中，搜索優(yōu)化的交錯網格有限差分算子在滿足譜誤差容限數值上限的情況下所能覆蓋的最大波數范圍。

具體實施時，為了減少優(yōu)化代價，在本實施例中，在求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的每個交錯網格有限差分權系數之前，還包括：設定模擬退火算法的自變量維度和各維度自變量的取值范圍，其中，自變量維度表示所述交錯網格有限差分權系數的個數，各維度自變量的取值范圍表示各個所述交錯網格有限差分權系數的振幅絕對值范圍，將自變量維度設定為得到所述交錯網格有限差分算子所用的網格點數(即窗函數)的長度的二分之一；將各維度自變量的取值范圍設定為區(qū)間【0，2】；設定各個所述交錯網格有限差分權系數的振幅圍繞中心位置呈衰減的震蕩。

具體的，本申請發(fā)明人根據sinc插值理論發(fā)現：(1)優(yōu)化的N點交錯網格有限差分算子中有N/2個交錯網格有限差分權系數，比如N＝16，則有8個系數要優(yōu)化；(2)各個交錯網格有限差分權系數的振幅絕對值均在[0，2]區(qū)間內；(3)交錯網格有限差分權系數的振幅應該是圍繞中心位置0呈衰減的震蕩，即|c_n|＞|c_n+1|，且|c_n||c_n+1|＜0，for n＝1,2,...,N/2。這樣我們可以通過只確定c₁到c_n而優(yōu)化整個交錯網格有限差分算子。在設定模擬退火算法的關鍵參數時，可以根據上述(1)和(2)模擬退火算法的自變量維度M和各維度自變量取值范圍，再用(3)原則對解進行限制大大減少了搜索范圍，使模擬退火法更高效。事實上通過觀察傳統(tǒng)差分算子，發(fā)現原始系數也滿足這三點約束。這就表明這三個規(guī)則對通常的差分算子都是合理、適用的。

具體的，采用模擬退火算法求解上述目標函數并進行地震波傳播模擬的過程如圖4所示：

(1)、優(yōu)化的交錯網格有限差分算子的最大譜誤差容限T為預設譜誤差容限；

(2)、將波數在范圍內均勻離散采樣，按照波數采樣點由小到大的順序，_K表示波數采樣點的序數，是以一個波數采樣點作為上限的預設波數范圍，從K＝1開始對每個波數采樣點執(zhí)行第3-4步，搜索滿足T的最優(yōu)交錯網格有限差分權系數。

(3)、用模擬退火算法進行全局搜索，其中要設定關鍵參數有：馬爾可夫(markov)鏈長度、起始溫度、終止溫度和步長因子；此外，還有自變量維度M和各維度自變量取值范圍，即2M階交錯網格有限差分算子要優(yōu)化的交錯網格有限差分系數個數及各交錯網格有限差分權系數的取值范圍。

(4)、若搜索到有交錯網格有限差分權系數使目標函數的極值小于預設譜誤差容限，則回到第二步，且k＝k+1，即將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解目標函數，進行下一個波數采樣點的全局尋優(yōu)；若使目標函數的極值大于預設譜誤差容限，則終止，取上一個波數采樣點的交錯網格有限差分系數作為最優(yōu)解輸出。

(5)、利用2M階優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波傳播的高精度模擬。

優(yōu)化的交錯網格有限差分算子更佳地逼近一階空間微分算子，1)可以用于各種介質中的高精度聲波模擬和逆時偏移成像，包括各向同性介質和各向異性介質；2)可以用于各種介質中的高精度彈性波模擬和逆時偏移成像；3)對于基于CPU/GPU協(xié)同并行的地震波模擬裝置，僅需要用優(yōu)化的交錯網格差分權系數代替原來的差分系數即可得到數值頻散更小的模擬結果。其中以彈性波的模擬為應用實例，給出二維情況下、各向同性介質中的一階交錯網格速度-應力方程如下：

式中，λ、μ是兩個拉梅常數，

具體實施時，在應用上述全局優(yōu)化的交錯網格有限差分正演模擬方法時，我們設定較緊的譜誤差容限T為0.0001，即以保證波在傳播很長距離后積累誤差仍然很小可以接受。然后采用模擬退火算法，求解交錯網格有限差分權系數，較快地收斂到全局最優(yōu)解，由此得到優(yōu)化的交錯網格有限差分算子：(1)使譜絕對誤差控制在萬分之一內，這是目前所達到的最小誤差；(2)盡可能地覆蓋更大的譜范圍，這也是目前此譜誤差容限下所達到的最大譜覆蓋范圍。從而在利用優(yōu)化的交錯網格有限差分算子進行波場模擬過程中，可以有效地降低空間離散引入的頻散誤差，提高模擬結果的精度。因此，我們稱本專利優(yōu)化的交錯網格有限差分算子相比于相同階數下的常規(guī)網格有限差分算子，具有更高的精度。從圖5(a)(圖5(a)中實線曲線表示常規(guī)交錯網格差分系數，實線曲線上的4、6、8、10等數字表示常規(guī)交錯網格差分系數的階數)中可以看出，優(yōu)化的8階交錯網格有限差分算子達到常規(guī)12階交錯網格差分算子的精度；優(yōu)化的12階交錯網格有限差分算子則超過了常規(guī)24階交錯網格差分算子的精度。從圖5(b)中可以發(fā)現，優(yōu)化的8階交錯網格有限差分算子在精度誤差0.0001下的譜覆蓋范圍遠超過16階常規(guī)交錯網格差分算子；而優(yōu)化的10、12階交錯網格有限差分算子在精度誤差0.0001下的譜覆蓋范圍遠超過常規(guī)28階交錯網格差分算子。因此，用低階優(yōu)化的交錯網格有限差分算子代替高階常規(guī)交錯網格有限差分算子，由于算子長度變短而減少了計算量和計算時間；另外，基于GPU并行計算平臺，差分算子長度越短則對應全局內存訪問量越少，減少延遲時間。即本專利優(yōu)化的交錯網格有限差分算子在保證了模擬精度的條件下，還可以有效地提高地震波模擬的效率。

以下結合具體示例來說明上述全局優(yōu)化的交錯網格有限差分正演模擬方法的有效性。

實例一、將優(yōu)化的8階交錯網格有限差分算子近似計算空間一階導數，而時間一階導數仍用常規(guī)二階交錯網格差分算子計算，求解公式(12)，模擬彈性波在均勻各向同性介質中的傳播，以驗證其有效性。

速度模型網格大小為500×500，橫向和縱向網格間隔均為10m，縱波速度為2500m/s，橫波速度為1700m/s，密度為常數1g/cm3；采用主頻為33Hz的雷克(Ricker)子波點源加載在垂向應力上，在模型中心激發(fā)，時間采樣間隔為1ms。如圖6(a)、6(b)、6(c)以及6(d)所示，可以發(fā)現在算子長度相同時(相同階數)，優(yōu)化的交錯網格有限差分算子模擬縱波和橫波的頻散明顯要比常規(guī)交錯網格差分算子模擬的少。從而，證明了本專利方案壓制數值頻散、提高地震波模擬精度方面的有效性。

實例二、用復雜的marmousi2模型驗證優(yōu)化的8階交錯網格有限差分算子的壓制頻散性能。

Marmousi2模型大小為961*561，其縱波速度模型如圖7所示，橫波速度用給出，密度為常數1g/cm3。采用主頻為38Hz的Ricker子波，在靠近地表處(2.9976，0.01249)Km激發(fā)；橫向和縱向網格間隔均為6.245m，時間采樣間隔為0.5ms。圖8(a)為優(yōu)化的8階交錯網格有限差分方法和常規(guī)方法的時間總長度為6s的垂直分量波場記錄，其中，左側是8階常規(guī)交錯網格有限差分模擬；右側為8階優(yōu)化交錯網格有限差分模擬；圖8(b)、8(c)分別是圖8(a)中區(qū)域1和2的放大圖，可以看出優(yōu)化的交錯網格有限差分方法有效地壓制了marmousi2中由橫波低速引起的數值頻散，提高了寬頻彈性波模擬的精度。

基于同一發(fā)明構思，本發(fā)明實施例中還提供了一種全局優(yōu)化的交錯網格有限差分正演模擬裝置，如下面的實施例所述。由于全局優(yōu)化的交錯網格有限差分正演模擬裝置解決問題的原理與全局優(yōu)化的交錯網格有限差分正演模擬方法相似，因此全局優(yōu)化的交錯網格有限差分正演模擬裝置的實施可以參見全局優(yōu)化的交錯網格有限差分正演模擬方法的實施，重復之處不再贅述。以下所使用的，術語“單元”或者“模塊”可以實現預定功能的軟件和/或硬件的組合。盡管以下實施例所描述的裝置較佳地以軟件來實現，但是硬件，或者軟件和硬件的組合的實現也是可能并被構想的。

圖9是本發(fā)明實施例的全局優(yōu)化的交錯網格有限差分正演模擬裝置的一種結構框圖，如圖9所示，該裝置包括：

有限差分算子獲取模塊901，用于獲得交錯網格有限差分算子；

目標函數求解模塊902，用于基于最大化范數建立目標函數，求解所述目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化所述交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，其中，所述交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的；

地震波正演模塊903，用于利用所述優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。

在一個實施例中，所述目標函數為：

其中，E(c)為所述目標函數的極值；是優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍；c_n是第n個所述交錯網格有限差分權系數；k_x是波數覆蓋范圍；T是優(yōu)化的交錯網格有限差分算子的最大譜誤差容限；_N是得到所述交錯網格有限差分算子所用的網格點數的長度；Δx是采樣間隔。

在一個實施例中，所述目標函數求解模塊具體用于采用模擬退火算法求解所述目標函數。

在一個實施例中，所述目標函數求解模塊，包括：誤差設定單元，用于設定所述優(yōu)化的交錯網格有限差分算子的最大譜誤差容限為預設譜誤差容限數值；目標函數求解單元，用于通過循環(huán)以下步驟求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的相應交錯網格有限差分權系數：將波數在0到尼奎斯特波數的范圍內離散采樣，按照波數采樣點由小到大的順序，先將最小的一個波數采樣點作為預設波數范圍的上限，該預設波數范圍的下限為0，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在該得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數，在得到的交錯網格有限差分權系數使得所述目標函數的極值大于所述預設譜誤差容限數值時，終止計算，將上一個波數采樣點采用模擬退火算法計算得到的交錯網格有限差分權系數輸出，作為所述優(yōu)化的交錯網格有限差分算子的交錯網格有限差分權系數；在得到的交錯網格有限差分權系數使得所述目標函數的極值小于所述預設譜誤差容限數值時，繼續(xù)再將下一個波數采樣點作為預設波數范圍的上限，在該預設波數范圍內采用模擬退火算法求解所述目標函數得到交錯網格有限差分權系數。

在一個實施例中，所述目標函數求解模塊，還包括：參數設定單元，用于在求解所述目標函數得到優(yōu)化所述交錯網格有限差分算子的每個交錯網格有限差分權系數之前，設定模擬退火算法的自變量維度和各維度自變量的取值范圍，其中，自變量維度表示所述交錯網格有限差分權系數的個數，各維度自變量的取值范圍表示各個所述交錯網格有限差分權系數的振幅絕對值范圍，將自變量維度設定為得到所述交錯網格有限差分算子所用的網格點數的長度的二分之一；將各維度自變量的取值范圍設定為區(qū)間【0，2】；設定各個所述交錯網格有限差分權系數的振幅圍繞中心位置呈衰減的震蕩。

在本發(fā)明實施例中，通過基于最大化范數建立目標函數，并求解目標函數得到用于優(yōu)化所述交錯網格有限差分算子的交錯網格有限差分權系數，采用求解得到的交錯網格有限差分權系數優(yōu)化交錯網格有限差分算子，得到優(yōu)化的交錯網格有限差分算子，進而利用優(yōu)化的交錯網格有限差分算子計算空間導數來進行地震波正演模擬。由于目標函數是基于最大化范數建立的，相對于現有技術中最小平方法建立的目標函數，使得優(yōu)化交錯網格有限差分算子的過程更加靈活；同時，由于通過求解目標函數得到的交錯網格有限差分權系數是在優(yōu)化的交錯網格有限差分算子的最大波數覆蓋范圍內且在優(yōu)化的交錯網格有限差分算子的最大譜誤差容限下計算得到的，有利于更加靈活地控制優(yōu)化的交錯網格有限差分算子的絕對譜誤差，使得優(yōu)化的交錯網格有限差分算子在較緊的譜誤差容限內還能夠覆蓋盡可能大的波數范圍，從而更有效地壓制數值頻散，提高地震模擬精度，有利于為逆時偏移、全波形反演提供高精度、高效的正演模擬方案，降低應用局限的問題。

顯然，本領域的技術人員應該明白，上述的本發(fā)明實施例的各模塊或各步驟可以用通用的計算裝置來實現，它們可以集中在單個的計算裝置上，或者分布在多個計算裝置所組成的網絡上，可選地，它們可以用計算裝置可執(zhí)行的程序代碼來實現，從而，可以將它們存儲在存儲裝置中由計算裝置來執(zhí)行，并且在某些情況下，可以以不同于此處的順序執(zhí)行所示出或描述的步驟，或者將它們分別制作成各個集成電路模塊，或者將它們中的多個模塊或步驟制作成單個集成電路模塊來實現。這樣，本發(fā)明實施例不限制于任何特定的硬件和軟件結合。

以上所述僅為本發(fā)明的優(yōu)選實施例而已，并不用于限制本發(fā)明，對于本領域的技術人員來說，本發(fā)明實施例可以有各種更改和變化。凡在本發(fā)明的精神和原則之內，所作的任何修改、等同替換、改進等，均應包含在本發(fā)明的保護范圍之內。