基於矩陣降秩的地震數據重建方法與流程

2023-05-18 18:05:47 1

本發明涉及油氣勘探開發技術領域，尤其涉及一種基於矩陣降秩的地震數據重建方法。

背景技術：

隨著我國油氣勘探開發程度的不斷提高，地震勘探對象逐漸從構造油氣藏轉向巖性油氣藏，地震勘探目標也變得日益複雜，這對地震數據的處理質量提出了更高要求。然而，在地震數據的採集和處理過程中，經常會遇到不規則缺失道地震數據。它的存在會對地震數據後續多道處理技術的正確運行產生不良影響，降低地震數據的處理品質。

不規則地震數據產生的原因主要來源於兩個方面：一是在地震數據的採集階段，由於存在障礙物、禁採區、海上地震拖纜羽狀漂移以及採集的經濟成本考慮等因素，使得地震數據在空間方向通常呈現稀疏或不規則分布；二是在地震數據預處理階段，由於剔除廢炮和廢道等因素也會引起地震數據在空間方向的不規則分布。不規則地震數據不但會對地震後續處理產生噪聲幹擾，而且會對地震多道處理技術的正確運行產生不良影響，其中以波動方程疊前深度偏移，自由表面多次波消除，譜估計和時移地震可重複性處理最為嚴重。對於空間稀疏採樣數據，在偏移處理中會產生空間假頻，這不但影響地震偏移速度場的精確建立，而且會降低偏移成像質量，使構造模糊，斷點不清晰，最終達不到油藏精確解釋和描述的目的。對於不規則採樣地震數據，會引起地下面元覆蓋次數不均勻，在疊加成像時出現採集腳印現象。此外，在時移地震油藏監測中由不規則採樣引起的前後兩次監測差異也會反映在處理剖面上，但這種差異實質並不代表地下流體的真實變化。總之，以上面臨的問題都需要藉助地震數據插值重建技術來解決。

此外，在我國3D地震勘探施工方式已經普及，3D地震勘探必然涉及疊前5D不規則缺失道地震數據的重建問題(5D分別指時間t、炮點坐標sx和sy、檢波器rx和ry，也可以指時間t、共中心點cmpx和cmpy、偏移距hx和hy，還可以指時間t、共中心點cmpx和cmpy、偏移距h和方位角θ)。

與2D和3D地震數據重建相比，5D數據重建可以利用更多維數中的已知道信息，並將這些信息作為約束條件來對缺失道進行重建，獲得更加精確的波場重建結果。因此，疊前5D地震數據重建是不規則缺失道重建領域的重點研究對象。

地震數據重建方法在經歷了基於函數變換重建、預測濾波重建和波動方程重建之後，目前已進入第四代基於降秩法的重建，基於矩陣降秩的地震數據重建是降秩法重建中最具代表性的一類方法。矩陣降秩法重建的基本原理為假設原始完整採樣的理想無噪聲地震數據可以被一個低秩矩陣表示，而隨機噪聲和不規則缺失採樣會增加矩陣的秩。因此，地震數據重建問題可以被看做一個高秩矩陣的降秩問題。目前矩陣降秩的方法主要有Cadzow濾波法(Cadzow filtering，簡稱CF)(Rank-reduction-based trace interpolation,SEG Technical Program Expanded Abstracts 2010,3829-3833)、基於隨機SVD分解的多道奇異譜分析方法(Multichannel singular spectrum analysis，簡稱MSSA)(Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis,Geophysics,2011,76(3):V25-V32)以及基於Lanczos雙對角分解的多道奇異譜分析方法(A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions,Geophysics,2013,78(1):V21–V30)。Cadzow濾波方法與多道奇異譜分析方法在方法原理上是一致的，區別在於方法來源不同。Cadzow濾波方法來源於語音信號增強領域，多道奇異譜分析方法來源於古氣候時間序列分析領域。這兩種方法在進行數據重建時需要在頻率域對每個頻率切片數據構建Hankel矩陣或者Toeplitz矩陣，然後運用截斷SVD分解算法，保留前幾個較大奇異值，捨棄較小的奇異值來實現矩陣降秩。由於截斷SVD分解計算速度慢，計算成本高，Oropeza和Sacchi提出運用隨機SVD分解算法取代截斷SVD分解來抬高降秩分解計算效率，對於2D和3D等低維數據重建而言這種方法有一定效果(Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis,Geophysics,2011,76(3):V25-V32；基於多道奇異譜分析的三維地震數據規則化方法,石油地球物理勘探，2014，49(5)：846-851)。但是，對於高維地震數據重建尤其是疊前5D地震數據重建而言，需要構建大型四重塊Hankel矩陣或者四重塊Toeplitz矩陣，隨機SVD分解算法因計算量大，計算成本高無法用於該大型四重塊Hankel矩陣或者四重塊Toeplitz矩陣的降秩分解。Gao等人(A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions,Geophysics,2013,78(1):V21–V30)提出運用Lanczos雙對角分解取代隨機SVD分解實現對四重塊Hankel或者四重塊Toeplitz矩陣的降秩，雖然運用Lanczos雙對角分解方法可以提高降秩計算效率，但是在對降秩後的四重Toeplitz矩陣沿對角線求平均時仍需先生成該大型低秩四重塊Toeplitz矩陣，導致數據的存儲量大，重建效率仍很低。

因此，5D地震數據重建除面臨計算量大的問題之外，還面臨著數據存儲量大的挑戰。如何尋找一種計算效率高，降秩計算速度快的降秩重建方法且能減小或者壓縮數據存儲量也是5D重建面臨的迫切問題。

技術實現要素：

本發明提供一種基於矩陣降秩的地震數據重建方法，以提高地震數據重建速度且減小數據存儲量。

本發明提供一種基於矩陣降秩的地震數據重建方法，包括：獲取地震頻率切片數據的四重塊Toeplitz矩陣，並存儲所述四重塊Toeplitz矩陣的部分行列元素以表示所述四重塊Toeplitz矩陣，所述部分行列元素包括所述四重塊Toeplitz矩陣的各不同元素；利用隨機QR降秩分解算法對所述部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理；利用無展開求平均算法對降秩處理後的所述四重塊Toeplitz矩陣的對角線元素求平均，得到所述地震頻率切片數據的重建數據；對所述地震頻率切片數據的重建數據做傅立葉反變換，得到時間域地震重建數據，用於油氣勘探。

一個實施例中，獲取地震頻率切片數據的四重塊Toeplitz矩陣，包括：對時間域地震數據做關於時間的傅立葉變換，得到頻率域地震數據，並從所述頻率域地震數據中提取得到所述地震頻率切片數據；利用所述地震頻率切片數據構建所述四重塊Toeplitz矩陣，或利用所述地震頻率切片數據構建四重塊Hankel矩陣，並將所述四重塊Hankel矩陣變換為所述四重塊Toeplitz矩陣。

一個實施例中，對時間域地震數據做關於時間的傅立葉變換，得到頻率域地震數據之前，還包括：通過坐標變換將地震勘探數據由炮-檢域轉換到共中心點-偏移距域，得到所述時間域地震數據。

一個實施例中，所述時間域地震數據為5D地震數據。

一個實施例中，所述部分行列元素為所述四重塊Toeplitz矩陣的第一行元素及第一列元素或最後一行元素及最後一列元素。

一個實施例中，利用隨機QR降秩分解算法對所述部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理，包括：通過將表示的四重塊Toeplitz矩陣乘以一隨機矩陣對表示的四重塊Toeplitz矩陣進行矩陣壓縮，得到壓縮矩陣，所述隨機矩陣的行數和列數分別為所述四重塊Toeplitz矩陣的列數和秩；利用隨機QR降秩分解算法對所述壓縮矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理。

一個實施例中，通過將表示的四重塊Toeplitz矩陣乘以一隨機矩陣對表示的四重塊Toeplitz矩陣進行矩陣壓縮，包括：用向量形式表示所述隨機矩陣，並利用基於四維傅立葉快速變換的快速乘法將表示的四重塊Toeplitz矩陣和向量形式表示的隨機矩陣相乘，以對表示的四重塊Toeplitz矩陣進行矩陣維數壓縮。

一個實施例中，利用隨機QR降秩分解算法對所述壓縮矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理，包括：對所述壓縮矩陣實施QR分解，得到由正交矩陣和上三角矩陣的乘積表示的壓縮矩陣，並存儲所述正交矩陣；將所述正交矩陣的共軛轉置矩陣和表示的四重塊Toeplitz矩陣相乘得到子矩陣，存儲所述子矩陣，並利用所述正交矩陣和所述子矩陣的乘積表示降秩後的所述四重塊Toeplitz矩陣。

一個實施例中，利用無展開求平均算法對降秩處理後的所述四重塊Toeplitz矩陣的對角線元素求平均，包括：利用所述正交矩陣和所述子矩陣計算降秩後的所述四重塊Toeplitz矩陣的對角線元素的平均。

一個實施例中，對所述地震頻率切片數據的重建數據做傅立葉反變換，得到時間域地震重建數據之前，還包括：計算對角線元素求平均的結果和所述地震頻率切片數據的差值並判斷該差值是否在設定誤差範圍內；若否，利用對角線元素求平均的結果重新獲取所述地震頻率切片數據的四重塊Toeplitz矩陣，並存儲重新獲取的四重塊Toeplitz矩陣的部分行列元素以表示重新獲取的四重塊Toeplitz矩陣，重新獲取的部分行列元素包括重新獲取的四重塊Toeplitz矩陣的各不同元素；利用隨機QR降秩分解算法對重新獲取的部分行列元素表示的重新獲取的四重塊Toeplitz矩陣進行降秩處理，以對重新獲取的四重塊Toeplitz矩陣進行降秩處理；利用無展開求平均算法對降秩處理後的重新獲取的四重塊Toeplitz矩陣的對角線元素求平均，得到重新獲取的對角線元素平均；計算重新獲取的對角線元素平均和所述對角線元素求平均的結果的差值並判斷該差值是否在所述設定誤差範圍內，若是，將重新獲取的對角線元素平均作為所述地震頻率切片數據的重建數據，若否，利用重新獲取的對角線元素平均再次重新獲取所述地震頻率切片數據的四重塊Toeplitz矩陣，並存儲再次重新獲取的四重塊Toeplitz矩陣的部分行列元素以表示再次重新獲取的四重塊Toeplitz矩陣，再次重新獲取的部分行列元素包括再次重新獲取的四重塊Toeplitz矩陣的各不同元素，利用隨機QR降秩分解算法對再次重新獲取的部分行列元素表示的再次重新獲取的四重塊Toeplitz矩陣進行降秩處理，以對再次重新獲取的四重塊Toeplitz矩陣進行降秩處理，利用無展開求平均算法對降秩處理後的再次重新獲取的四重塊Toeplitz矩陣的對角線元素求平均，計算再次重新獲取的對角線元素平均和所述重新獲取的對角線元素平均的差值，依次迭代進行，直到所得差值在所述設定誤差範圍內。

本發明實施例的基於矩陣降秩的地震數據重建方法，通過僅存儲四重塊Toeplitz矩陣的部分行列元素來表示整個四重塊Toeplitz矩陣能夠降低地震數據重建過程中的數據存儲量；利用隨機QR降秩分解算法可以實現對部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以此可以減少降秩計算量，提高地震數據重建速度；利用無展開求平均算法對降秩處理後的四重塊Toeplitz矩陣的對角線元素求平均值，無需計算並存儲降秩處理後的所述四重塊Toeplitz矩陣的具體形式，無需存儲大型的低秩四重塊Toeplitz矩陣，可以減少數據存儲量，提高地震數據重建效率。本發明的方法能夠減少數據儲存量和計算量且不降低重建質量。本發明實施例的方法不僅能夠用於五維地震數據的重建，還可以很好地用於N維地震數據的重建，N≥2。

附圖說明

為了更清楚地說明本發明實施例或現有技術中的技術方案，下面將對實施例或現有技術描述中所需要使用的附圖作簡單地介紹，顯而易見地，下面描述中的附圖僅僅是本發明的一些實施例，對於本領域普通技術人員來講，在不付出創造性勞動的前提下，還可以根據這些附圖獲得其他的附圖。在附圖中：

圖1是本發明實施例的基於矩陣降秩的地震數據重建方法的流程示意圖；

圖2是本發明一實施例中獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法流程示意圖；

圖3是本發明另一實施例中獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法流程示意圖；

圖4是本發明一實施例中利用隨機QR降秩分解算法對部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理的方法流程示意圖；

圖5是本發明一實施例中利用隨機QR降秩分解算法對壓縮矩陣進行降秩處理的方法流程示意圖；

圖6是本發明另一實施例的基於矩陣降秩的地震數據重建方法的流程示意圖；

圖7是本發明另一個實施例的基於矩陣降秩的地震數據重建方法的流程示意圖；

圖8是分別利用現有方法和本發明一實施例的基於矩陣降秩的地震數據方法進行重建所耗費計算時間的對比示意圖；

圖9是分別利用現有方法和本發明一實施例的基於矩陣降秩的地震數據方法進行重建的重建質量因子對比示意圖；

圖10是本發明一實施例中原始完整數據的cmpx面元切片示意圖；

圖11是本發明一實施例中50％不規則缺失道數據的cmpx面元切片示意圖；

圖12是本發明一實施例中重建數據的cmpx面元切片示意圖；

圖13是圖10所示原始完整數據與圖12所示重建數據的差剖面；

圖14是本發明一實施例中原始實際採集的觀測數據的cmpy面元切片示意圖；

圖15是利用現有基於Lanczos雙對角分解的多道奇異譜分析方法的重建結果的cmpy面元切片示意圖；

圖16是利用本發明實施例的方法的重建結果的cmpy面元切片示意圖。

具體實施方式

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

在現有矩陣降秩類重建方法中最具代表性的是基於隨機SVD分解的多道奇異譜分析方法，該方法可以歸結為：(1)對(t,x1,x2)域3D地震數據做關於時間t的傅立葉變換到頻率-空間域，得到(f,x1,x2)域數據，這裡x1和x2表示空間維變量；(2)對每一個頻率切片數據構建二重塊Hankel矩陣；(3)對二重塊Hankel矩陣利用隨機SVD分解算法進行降秩，得到降秩後的二重塊Hankel矩陣；(4)對降秩後的低秩二重塊Hankel矩陣沿反對角線元素求平均得到重建數據；(5)對該頻率切片重複步驟(2)，(3)和(4)，直到該切片滿足重建誤差要求；(6)對所有頻率成分按步驟(2)-(5)進行重建，最後對頻率f做傅立葉反變換到時間域完成重建。

發明人發現，基於隨機SVD分解的多道奇異譜分析方法存在如下缺陷：

1、該技術方案的第二歩需要生成並存儲二重塊Hankel矩陣，對於實際數據在重建時需要構建大型塊Hankel矩陣，導致存儲量大幅增加，佔用大量內存。然而，Hankel矩陣具有沿反對角線元素相等的特殊結構，只需存儲矩陣的第1行和最後1列就可以表示整個矩陣，不需要將整個矩陣生成並且存儲它的全部元素。例如，對於一個大小為10×10的Hankel矩陣，若直接生成和存儲該矩陣需要存儲100個元素，而利用本發明的方案只需存儲19個元素。

2、在步驟(3)矩陣降秩階段，該技術採用隨機SVD分解算法，雖然該算法與標準的截斷SVD分解算法相比計算量減少，但是仍需依賴SVD分解，對大型塊Hankel矩陣而言，降秩計算量依然很大。例如，對一個秩為k,大小為m×n的二重塊Hankel矩陣(m≤n)利用該技術進行降秩，計算量為8mn+8nk2+5mk2,而本發明提出的隨機QR分解技術計算量為km2+k(m+n-1)log2(m+n-1)。若取m＝n＝10，k＝3，則該技術的運算量為3570次，而本發明的技術只需585次運算。

3、在步驟(4)沿反對角線求平均的過程中，該技術是先生成和存儲降秩後的低秩二重塊Hankel矩陣，這會佔用大量內存。其實該二重塊Hankel矩陣可以不用存儲，沿反對角線進行求平均的元素可以直接從降秩分解後的子矩陣中提取。

4、該技術只能對疊前3D地震數據進行插值，無法對高維的5D數據進行插值重建處理。

基於Lanczos雙對角分解的多道奇異譜分析重建方法(A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions,Geophysics,2013,78(1):V21–V30)，採用Lanczos雙對角分解算法來取代截斷SVD分解算法，提高降秩計算效率，並且採用四維快速傅立葉變換來實現四重塊Toeplitz矩陣與向量的快速乘法，避免矩陣與向量直接相乘。但是，該技術在對降秩後的低秩四重塊Toeplitz矩陣沿對角線元素求平均時需要先生成並存儲該大型四重塊Toeplitz矩陣，然後再沿對角線求平均，得到重建數據。生成並存儲該大型四重塊Toeplitz矩陣必然會佔用大量內存，而這一過程完全可以省略，直接從降秩分解得到的子矩陣中提取出對角線元素並求平均，從而節約數據存儲量。該技術的實施方案可以歸結為：(1)將(t,x1,x2,x3,x4)域5D地震數據做關於時間t的傅立葉變換到頻率域得(f,x1,x2,x3,x4)域地震數據，這裡x1,x2,x3,x4表示空間維變量，既可以指炮點x坐標Sx，炮點y坐標Sy，檢波器x坐標rx和檢波器y坐標ry，也可以指共中心點x坐標CMPx，共中心點y坐標CMPy，偏移距x坐標hx和偏移距y坐標hy)；(2)對每一個頻率切片數據構建四重塊Toeplitz矩陣；(3)對該四重塊Toeplitz矩陣利用Lanczos雙對角分解進行降秩，得到降秩後的低秩四重塊Toeplitz矩陣；(4)對低秩四重塊Toeplitz矩陣沿對角線元素求平均，得到該地震頻率切片對應的重建數據；(5)對該地震頻率切片重複步驟(2)，(3)和(4)，直到該切片滿足重建誤差要求，完成對該頻率切片的數據重建；(6)對所有的頻率成分按步驟(2)-(5)進行重建，最後對頻率f做傅立葉反變換到時間域，完成重建。

發明人發現，基於Lanczos雙對角分解的多道奇異譜分析重建方法存在如下缺陷：

1、該技術直接對由頻率切片數據構建的四重塊Toeplitz矩陣實施降秩，而四重塊Toeplitz矩陣元素存在冗餘，完全可以採用類似於現有技術1的做法，先對該大型四重Toeplitz矩陣乘以一個隨機矩陣進行矩陣壓縮，然後對壓縮後的矩陣進行降秩分解，進而可以減小計算量，但是該技術未採用該做法。

2、該技術步驟(3)採用Lanczos雙對角分解算法進行降秩，與傳統的截斷SVD分解算法相比，計算效率有提高，但還可以進一步被提高。對於一個秩為k的m×n型四重塊Toeplitz矩陣，Lanczos分解的降秩計算量為O(km2)階，若採用本發明的隨機QR分解算法，計算量僅為O(mk2)，這裡k＜＜m。

3、在步驟(4)對降秩後的低秩四重塊Toeplitz矩陣沿對角線求平均時該技術仍需先生成一個大型四重塊Toeplitz矩陣，然後再沿對角線求平均，此處生成和存儲該四重塊Toeplitz矩陣需要佔用大量存儲空間。本發明採用無展開求平均技術，在求平均步驟可以不用生成和存儲大型該四重塊Toeplitz矩陣，而是直接從降秩分解的子矩陣中提取出對角線元素並對它們求平均。

基於上述發現，針對疊前5D重建面臨的降秩分解計算量大，計算成本高以及數據存儲量大的問題，本發明提出了一種基於矩陣降秩的地震數據重建方法。該方法是一種非SVD分解方法，可以有效解決矩陣降秩因嚴重依賴SVD分解算法導致的計算量大、計算成本高、無法應用於大規模工業實際生產數據插值重建等問題。

圖1是本發明實施例的基於矩陣降秩的地震數據重建方法的流程示意圖。如圖1所示，本發明實施例的基於矩陣降秩的地震數據重建方法，可包括步驟：

S110：獲取地震頻率切片數據的四重塊Toeplitz矩陣，並存儲所述四重塊Toeplitz矩陣的部分行列元素以表示所述四重塊Toeplitz矩陣，所述部分行列元素包括所述四重塊Toeplitz矩陣的各不同元素；

S120：利用隨機QR降秩分解算法對所述部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理；

S130：利用無展開求平均算法對降秩處理後的所述四重塊Toeplitz矩陣的對角線元素求平均，得到所述地震頻率切片數據的重建數據；

S140：對所述地震頻率切片數據的重建數據做傅立葉反變換，得到時間域地震重建數據，用於油氣勘探。

在上述步驟S110中，該地震頻率切片數據可以是頻率域某頻率下的地震數據。該四重塊Toeplitz矩陣可以由地震頻率切片數據直接構建，或由地震頻率切片數據構建的其他類型矩陣，例如四重塊Hankel矩陣，通過變換得到。可以根據四重塊Toeplitz矩陣的結構特點，僅存儲四重塊Toeplitz矩陣中的部分行列元素，該部分行列元素中的元素可以覆蓋四重塊Toeplitz矩陣中的各不同元素。利用該部分行列元素可以表示整個四重塊Toeplitz矩陣。

在上述步驟S120中，利用隨機QR降秩分解算法對所述部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理時，可以僅涉及上述部分行列元素，例如僅涉及四重塊Toeplitz矩陣的第一行及第一列，以此可以實現對整個四重塊Toeplitz矩陣進行降秩處理。

在上述步驟S130中，對降秩處理後的所述部分行列元素表示的四重塊Toeplitz矩陣的對角線元素求平均，既可以實現對降秩處理後的所述四重塊Toeplitz矩陣的對角線元素求平均。降秩處理後的低秩四重塊Toeplitz矩陣可以用兩個或多個矩陣相乘的形式表示，利用無展開求平均算法可以直接計算得到降秩處理後的低秩四重塊Toeplitz矩陣對角線元素的平均值，而無需先計算得到降秩處理後的低秩四重塊Toeplitz矩陣的具體形式後再對具體形式中的對角線元素求平均值。若降秩處理後的低秩四重塊Toeplitz矩陣的對角線元素的平均值滿足重建條件，則可得到所述地震頻率切片數據的重建數據。

在上述步驟S140中，地震勘探數據可以對應多個上述地震頻率切片數據。每個地震頻率切片數據可以利用上述步驟S110～步驟S130得到相應的重建數據。對各地震頻率切片數據的重建數據做傅立葉反變換，可得到時間域地震重建數據。該時間域地震重建數據可以填補最初的地震缺失道數據，並可以壓制包含在最初地震勘探數據中的隨機噪聲，從而使該時間域地震重建數據相比於原來含有缺失道地震數據具有較高的品質，將重建後的數據應用於油氣勘探，勘探結果會更準確。

本發明實施例中，通過僅存儲四重塊Toeplitz矩陣的部分行列元素來表示整個四重塊Toeplitz矩陣能夠降低地震數據重建過程中的數據存儲量；利用隨機QR降秩分解算法可以實現對部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以此可以減小降秩計算量，提高地震數據重建速度；利用無展開求平均算法對降秩處理後的四重塊Toeplitz矩陣的對角線元素求平均，無需計算並存儲降秩處理後的低秩四重塊Toeplitz矩陣的具體形式，無需存儲大型的低秩四重塊Toeplitz矩陣，可以減少數據存儲量，提高地震數據重建效率。

一些實施例中，所述部分行列元素可為所述四重塊Toeplitz矩陣的第一行元素及第一列元素或最後一行元素及最後一列元素。利用四重塊Toeplitz矩陣的第一行元素及第一列元素或四重塊Toeplitz矩陣的最後一行元素及最後一列元素即可表示整個四重塊Toeplitz矩陣。

本實施例中，僅存儲四重塊Toeplitz矩陣的第一行元素及第一列元素或最後一行元素及最後一列元素，可以顯著降低數據存儲量。

其他一些實施例中，在上述步驟S110中，可獲取由地震頻率切片數據構建的四重塊Hankel矩陣，可以只需存儲四重塊Hankel矩陣的第1行元素和最後1列元素。本實施例中，針對四重塊Hankel矩陣沿反對角線元素相等的結構特徵，只存儲矩陣的第一行和最後一列，可以減少數據存儲量。

圖2是本發明一實施例中獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法流程示意圖。如圖2所示，在上述步驟S110中，獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法，可包括步驟：

S111：對時間域地震數據做關於時間的傅立葉變換，得到頻率域地震數據，並從所述頻率域地震數據中提取得到所述地震頻率切片數據；

S112：利用所述地震頻率切片數據構建所述四重塊Toeplitz矩陣，或利用所述地震頻率切片數據構建四重塊Hankel矩陣，並將所述四重塊Hankel矩陣變換為所述四重塊Toeplitz矩陣。

在上述步驟S111中，可以通過各種不同方式從所述頻率域地震數據中提取得到所述地震頻率切片數據，例如可以根據傅立葉變換的長度參數、時間域地震數據沿時間維的採樣點數和採樣間隔等確定頻率切片的提取間隔，可以按頻率從小到大的順序提取地震頻率切片數據。

在上述步驟S112中，通過多種不同方法得到的所述地震頻率切片數據的所述四重塊Toeplitz矩陣，可以充分利用四重塊Toeplitz矩陣的結構特點，可僅存儲四重塊Toeplitz矩陣的第一行元素及第一列元素或最後一行元素及最後一列元素，從而減少四重塊Toeplitz矩陣的元素存儲量。

圖3是本發明另一實施例中獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法流程示意圖。如圖3所示，圖2所示的獲取地震頻率切片數據的四重塊Toeplitz矩陣的方法，在步驟S111之前，即對時間域地震數據做關於時間的傅立葉變換，得到頻率域地震數據之前，還可包括步驟：

S113：通過坐標變換將地震勘探數據由炮-檢域轉換到共中心點-偏移距域，得到所述時間域地震數據。

本實施例中，通過將地震勘探數據由炮-檢域轉換到共中心點-偏移距域，並基於共中心點-偏移距域的時間域地震數據進行後續降秩處理，可以充分利用共中心點-偏移距域的時間域地震數據的低秩特徵，進一步減少四重塊Toeplitz矩陣降秩處理的計算量。

一些實施例中，所述時間域地震數據為5D地震數據。本實施例中，對5D(五維)地震數據進行重建可以利用更多維的已知地震道信息，將這些信息作為約束來進行重建，可以獲得更加精確的地震波場重建結果。在其他實施例中，本發明的方法可以用於其他各種維數的地震數據的重建，例如3D(三維)地震數據。

圖4是本發明一實施例中利用隨機QR降秩分解算法對部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理的方法流程示意圖。如圖4所示，在上述步驟S120中，利用隨機QR降秩分解算法對所述部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理的方法，可包括步驟：

S121：通過將表示的四重塊Toeplitz矩陣乘以一隨機矩陣對表示的四重塊Toeplitz矩陣進行矩陣壓縮，得到壓縮矩陣，所述隨機矩陣的行數和列數分別為所述四重塊Toeplitz矩陣的列數和秩；

S122：利用隨機QR降秩分解算法對所述壓縮矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理。

在上述步驟S121中，該隨機矩陣可需滿足兩個條件，其一是該矩陣的各列線性獨立，互不相干；其二是該矩陣的行數和列數分別為所述四重塊Toeplitz矩陣的列數和秩。該隨機矩陣具有較小的維數，通過將上述部分行列元素表示的四重塊Toeplitz矩陣乘以該隨機矩陣，可以對上述部分行列元素表示的四重塊Toeplitz矩陣進行矩陣壓縮，得到壓縮矩陣。在上述步驟S122中，與直接對四重塊Toeplitz矩陣進行降秩處理相比，利用隨機QR降秩分解算法對所述壓縮矩陣進行降秩處理，能夠減少對四重塊Toeplitz矩陣降秩的冗餘，從而可以進一步減少將秩處理的計算量。

一些實施例中，在上述步驟S122中，通過將表示的四重塊Toeplitz矩陣乘以一隨機矩陣對表示的四重塊Toeplitz矩陣進行矩陣壓縮的方法，具體實施方式可為：用向量形式表示所述隨機矩陣，並利用基於四維傅立葉快速變換的快速乘法將表示的四重塊Toeplitz矩陣和向量形式表示的隨機矩陣相乘，以對表示的四重塊Toeplitz矩陣進行矩陣維數壓縮。

本實施例中，利用基於四維傅立葉快速變換的快速乘法計算上述部分行列元素表示的四重塊Toeplitz矩陣和向量形式表示的隨機矩陣的乘積，能夠提高對四重塊Toeplitz矩陣進行矩陣壓縮的計算速度。

一個實施例中，上述實施例的基於四維傅立葉快速變換的快速乘法可以參考「A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions」(Geophysics,2013,78(1)：V21–V30)中基於4D FFT的四重塊Toeplitz矩陣和向量的快速乘法原理實施。

圖5是本發明一實施例中利用隨機QR降秩分解算法對壓縮矩陣進行降秩處理的方法流程示意圖。如圖5所示，在上述步驟S122中，利用隨機QR降秩分解算法對所述壓縮矩陣進行降秩處理，以對所述四重塊Toeplitz矩陣進行降秩處理的方法，可包括步驟：

S1221：對所述壓縮矩陣實施QR分解，得到由正交矩陣和上三角矩陣的乘積表示的壓縮矩陣，並存儲所述正交矩陣；

S1222：將所述正交矩陣的共軛轉置矩陣和表示的四重塊Toeplitz矩陣相乘得到子矩陣，存儲所述子矩陣，並利用所述正交矩陣和所述子矩陣的乘積表示降秩後的所述四重塊Toeplitz矩陣。

在上述步驟S1221中，對上述壓縮矩陣實施QR分解，可以分解得到正交矩陣，記為Q和上三角矩陣，記為R，利用該正交矩陣Q和該上三角矩陣R的乘積形式可以表示該壓縮矩陣。在上述步驟S1222中，利用該正交矩陣R、該正交矩陣的共軛轉置矩陣，記為Q*，及上述部分行列元素表示的四重塊Toeplitz矩陣三者的乘積可以表示降秩後的低秩四重塊Toeplitz矩陣。而該正交矩陣和上述子矩陣維數較小，通過存儲該正交矩陣和該子矩陣，並利用該正交矩陣和該子矩陣表示降秩後的所述四重塊Toeplitz矩陣，可以避免直接計算和存儲該大型低秩四重塊Toeplitz矩陣，以此可以進一步減少數據存儲量和計算量。

一些實施例中，可以利用基於四維傅立葉快速變換的矩陣和向量快速乘法將所述正交矩陣的共軛轉置矩陣Q*和表示的四重塊Toeplitz矩陣相乘得到子矩陣，以此可以進一步提高地震數據重建的速度。

一些實施例中，在上述步驟S130中，利用無展開求平均算法對降秩處理後的所述四重塊Toeplitz矩陣的對角線元素求平均的方法，具體實施方式可以表示為：利用所述正交矩陣和所述子矩陣計算降秩後的所述四重塊Toeplitz矩陣的對角線元素的平均。

本實施例中，在不生成和不存儲降秩後的大型低秩四重塊Toeplitz矩陣的情況下，利用事先存儲的上述正交矩陣和上述子矩陣直接計算降秩後的所述低秩四重塊Toeplitz矩陣的對角線元素平均值，無需計算及存儲大型的降秩後的低秩四重塊Toeplitz矩陣，以此可以進一步減少數據存儲量和計算量。

一個實施例中，對降秩後的所述四重塊Toeplitz矩陣沿對角線元素求平均值可表示為：

其中，

在上式中，D(fi,i1,i2,i3,i4)表示降秩後沿所述四重塊Toeplitz矩陣的對角線元素求取的平均值，fi表示第i個地震頻率切片數據的頻率，i1、i2、i3及i4表示第i個地震頻率切片數據分別在共中心點x坐標cmpx、共中心點y坐標cmpy、偏移距x坐標hx及偏移距y坐標hy的數值；為根據正交矩陣Q構建的四重塊Toeplitz矩陣，為正交矩陣Q的共軛轉置矩陣Q*和四重塊Toeplitz矩陣可記為A(4)相乘得到的子矩陣M的列向量形式，M＝Q*A(4)，r為所述四重塊Toeplitz矩陣的秩，l為正整數，1≤l≤r；四重塊Toeplitz矩陣大小為N4×K4，其中的元素為根據正交矩陣構建的N3×K3的三重塊Toeplitz矩陣，p4＝1,2,...,L4；該三重塊Toeplitz矩陣的每個元素為根據正交矩陣構建的N2×K2的二重塊Toeplitz矩陣，p3＝1,2,...,L3；該二重塊Toeplitz矩陣的每個元素為根據正交矩陣構建的N1×K1的Toeplitz矩陣，p2＝1,2,...,L2；該Toeplitz矩陣的每個元素是正交矩陣Q中第p1行且第l列的元素，p1＝1,2,...,L1；列向量大小為K4×1，其中的每個元素是根據子矩陣M構建的大小為K3×1的列向量，ν4＝1,2,...,K4；向量的每個元素是根據子矩陣M構建的大小為K2×1的列向量，ν3＝1,2,...,K3；向量中的每個元素是根據子矩陣M構建的大小為K1×1的列向量，ν2＝1,2,...,K2；向量中的每個元素是子矩陣M的第l行且第ν1列元素，ν1＝1,2,...,K1；Lj+Kj-1＝Nj，符號表示取整數部分，Nj為地震數據D(fi,i1,i2,i3,i4)在第j維空間的方向長度，j＝1,2,3,4。

其他一些實施例中，由地震頻率切片數據構建的四重塊Hankel矩陣，可以不存儲該四重塊Hankel矩陣，只需存儲該四重塊Hankel矩陣的第一行和最後一列。而沿四重塊Hankel矩陣反對角線進行求平均的元素可以直接從降秩分解後的低秩四重塊Hankel矩陣的子矩陣中提取，然後求平均。本實施例中，直接從降秩分解得到的正交矩陣Q和子矩陣M中提取出對角線元素並求平均，可以減少數據存儲量和計算量。

圖6是本發明另一實施例的基於矩陣降秩的地震數據重建方法的流程示意圖。如圖6所示，圖1所示的基於矩陣降秩的地震數據重建方法，在步驟S140之前，即對所述地震頻率切片數據的重建數據做傅立葉反變換，得到時間域地震重建數據之前，還可包括步驟：

S150：計算對角線元素求平均的結果和所述地震頻率切片數據的差值並判斷該差值是否在設定誤差範圍內；

S160：若否，利用對角線元素求平均的結果重新獲取所述地震頻率切片數據的四重塊Toeplitz矩陣，並存儲重新獲取的四重塊Toeplitz矩陣的部分行列元素以表示重新獲取的四重塊Toeplitz矩陣，重新獲取的部分行列元素包括重新獲取的四重塊Toeplitz矩陣的各不同元素；

S170：利用隨機QR降秩分解算法對重新獲取的部分行列元素表示的重新獲取的四重塊Toeplitz矩陣進行降秩處理，以對重新獲取的四重塊Toeplitz矩陣進行降秩處理；

S180：利用無展開求平均算法對降秩處理後的重新獲取的四重塊Toeplitz矩陣的對角線元素求平均，得到重新獲取的對角線元素平均；

S190：計算重新獲取的對角線元素平均和所述對角線元素求平均的結果的差值並判斷該差值是否在所述設定誤差範圍內，若是，將重新獲取的對角線元素平均作為所述地震頻率切片數據的重建數據，若否，利用重新獲取的對角線元素平均再次重新獲取所述地震頻率切片數據的四重塊Toeplitz矩陣，並存儲再次重新獲取的四重塊Toeplitz矩陣的部分行列元素以表示再次重新獲取的四重塊Toeplitz矩陣，再次重新獲取的部分行列元素包括再次重新獲取的四重塊Toeplitz矩陣的各不同元素，利用隨機QR降秩分解算法對再次重新獲取的部分行列元素表示的再次重新獲取的四重塊Toeplitz矩陣進行降秩處理，以對再次重新獲取的四重塊Toeplitz矩陣進行降秩處理，利用無展開求平均算法對降秩處理後的再次重新獲取的四重塊Toeplitz矩陣的對角線元素求平均，計算再次重新獲取的對角線元素平均和所述重新獲取的對角線元素平均的差值，依次迭代進行，直到所得差值在所述設定誤差範圍內。

本實施例中，當後一次計算得到的對角線元素平均值與前一次計算得到的對角線元素平均值(或地震頻率切片數據)的差值不滿足重建條件，即不小於或等於設定誤差範圍，例如10-4，則利用迭代求解方法，繼續對重建的結果進行重建直到滿足重建條件，以此可以獲得更佳的地震數據重建結果。

圖7是本發明另一個實施例的基於矩陣降秩的地震數據重建方法的流程示意圖。如圖7所示，本發明實施例的基於矩陣降秩的地震數據重建方法，可包括步驟：

(1)對時間域5D地震數據做關於時間t的傅立葉變換到頻率域；

(2)對每一個頻率切片數據構建四重塊Toeplitz矩陣並只存儲該四重塊Toeplitz矩陣的第一行和第一列元素來表示整個矩陣，以此可以達到節約計算機內存，減小數據存儲量的目的；

(3)採用隨機QR快速降秩分解算法對四重塊Toeplitz矩陣進行降秩；

(4)採用無展開求平均算法對降秩後的四重塊Toeplitz矩陣沿對角線元素求平均，得到該地震頻率切片的重建數據；

(5)對該地震頻率切片按照步驟(2)、(3)和(4)進行迭代重建，直至滿足重建誤差條件；

(6)對所有頻率成分按步驟(2)-(5)進行重建，最後做關於頻率f的傅立葉反變換到時間域得到時間域重建數據，完成重建。

一個實施例中，利用隨機QR分解進行矩陣降秩的地震數據重建方法，可包括步驟：

(1)在共中心點-偏移距域進行地震數據重建。首先可通過坐標變換將3D地震勘探中採集的疊前炮-檢域5D數據D(t,sx,sy,rx,ry)轉換到共中心點-偏移距域的5D數據D(t,cmpx,cmpy,hx,hy)。例如，坐標變換公式可表示為hx＝sx-rx和hy＝sy-ry，然後再進行重建；

(2)對於共中心點-偏移距域的5D數據D(t,cmpx,cmpy,hx,hy)做關於時間t的傅立葉變換到頻率域得到頻率域5D地震數據D(f,cmpx,cmpy,hx,hy)；

(3)提取第i個頻率切片數據D(fi,cmpx,cmpy,hx,hy)，設置迭代次數初始值k＝1，對該頻率切片數據進行迭代重建，fi是第i個頻率切片數據對應的頻率值；

(4)將由頻率切片數據D(fi,cmpx,cmpy,hx,hy)構建的四重塊Toeplitz矩陣記為A(4)，充分挖掘和利用四重塊Toeplitz矩陣沿對角線元素相等的結構特徵，只存儲矩陣A(4)的第一行和第一列來表示整個四重塊Toeplitz矩陣A(4)，以此可以達到壓縮數據存儲量的目的；

(5)對第一行和第一列表示的四重塊Toeplitz矩陣A(4)運用隨機QR分解進行降秩，得到頻率切片fi對應的重建數據。具體過程例如可為：設A(4)是一個秩為r(r<m≤n)，大小為m×n的四重塊Toeplitz矩陣，Ω是一個列向量互不相干且大小為n×r的隨機矩陣。將四重塊Toeplitz矩陣A(4)與隨機矩陣Ω相乘得到壓縮後的矩陣B，B＝A(4)Ω，此時壓縮後的矩陣B是一個維數較小的m×r型壓縮矩陣。在計算壓縮後的矩陣B時，可運用基於四維傅立葉快速變換的四重塊Toeplitz矩陣和向量的快速乘法來實現矩陣A(4)和Ω的相乘，而不直接將A(4)與Ω相乘，以此可提高計算效率。對壓縮後的矩陣B實施QR分解得B＝QR，Q為正交矩陣，R為上三角矩陣。將降秩後的低秩四重塊Toeplitz矩陣記為則取M＝Q*A(4)，則低秩四重塊Toeplitz矩陣需要說明的是此處只是一個符號標記，具體實施時不需要真正計算和存儲該低秩四重塊Toeplitz矩陣而只需計算和存儲子矩陣Q和M。在計算子矩陣M時可運用基於四維快速傅立葉變換的四重塊Toeplitz矩陣和向量的快速乘法技術來實現正交矩陣Q的共軛轉置矩陣Q*和四重塊Toeplitz矩陣A(4)相乘。該快速乘法可參照「A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions」(Geophysics,2013,78(1)：V21–V30)。對低秩四重塊Toeplitz矩陣可採用無展開快速求平均算法，沿低秩四重塊Toeplitz矩陣對角線元素求平均值，得到重建數據。無展開快速求平均算法例如可以表示為，

其中，變量i1,i2,i3和i4分別表示空間坐標cmpx,cmpy,hx和hy，設Nj表示地震數據D在第j維空間方向的長度，Lj和Kj滿足Lj+Kj-1＝Nj，符號表示取整數部分，j＝1,2,3,4；為一個四重塊Toeplitz矩陣，其首行元素由一個大小為K4×1的向量構成，其首列元素由一個大小為N4×1的向量構成，為一個大小為K4×1的列向量矩陣和向量中元素具體表達式為，

其中中的每個元素p4＝1,2,...,L4為一個三重塊Toeplitz矩陣，該三重塊Toeplitz矩陣中每個元素又可表示為p3＝1,2,...,L3。為一個二重塊Toeplitz矩陣，該二重塊Toeplitz矩陣中的每個元素可表示為p2＝1,2,...,L2。為一個普通Toeplitz矩陣，該Toeplitz矩陣中的每個元素可表示為p1＝1,2,...,L1。是位於正交矩陣Q中第p1行和第l列的元素。同樣，每個ν4＝1,2,...,K4都是一個K3×1的列向量。向量中每個元素ν3＝1,2,...,K3又是一個K2×1的列向量，中每個元素ν2＝1,2,...,K2是一個K1×1的列向量。向量中每個元素ν1＝1,2,...,K1取自子矩陣M的第l行和第ν1列元素。計算可運用基於四維快速傅立葉變換的四重塊Toeplitz矩陣和向量快速乘法技術來實現，該快速乘法的具體計算原理可參照「A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions」(Geophysics,2013,78(1)：V21–V30)。對D(fi,i1,i2,i3,i4)中的缺失點數據而言，做一次降秩分解重建後該點的振幅能量值並不能完全被恢復，需要不斷按照步驟(3)、(4)和(5)進行反覆迭代重建，直至滿足重建誤差條件。

一個實施例中，對頻率fi對應的地震頻率切片數據而言，第k次迭代重建的表達式可表示為，

D(k+1)(fi,i1,i2,i3,i4)＝αDobs(fi,i1,i2,i3,i4)+(I-αS)Dk(fi,i1,i2,i3,i4),k＝1,...,Niter.

其中，Niter為最大迭代次數，Dobs(fi,i1,i2,i3,i4)為原始輸入的觀測數據，α為權值因子，α∈(0,1]，S為採樣算子，其元素為0和1，其中0表示缺失點，1表示已知點；Dk(fi,i1,i2,i3,i4)表示第k次迭代的重建數據，D(k+1)(fi,i1,i2,i3,i4)表示第k+1迭代的重建數據，I表示元素值全為1的四維數組。

(6)對所有頻率切片數據D(fi,i1,i2,i3,i4)按步驟(2)-(5)進行重建。最後，對頻率f做傅立葉反變換到時間域，得到時間域重建數據，最終完成重建。

本發明實施例的方法，因對四重塊Toeplitz矩陣採用了隨機QR快速分解降秩方法，無展開求平均方法，基於四維傅立葉快速變換的四重塊Toeplitz矩陣和向量快速乘法計算技術以及依據四重塊Toeplitz矩陣結構特點只存儲矩陣的第一行和第一列來表示整個四重塊Toeplitz矩陣的有效存儲方法，使得矩陣降秩重建的計算效率得到明顯提升。在提升計算速度的同時還能保持重建質量。本發明實施例的方法是一種非SVD降秩重建方法，可以有效解決目前矩陣降秩重建方法因嚴重依賴SVD分解，計算量大而無法應用於大規模工業數據插值重建處理的瓶頸問題。

一個實施例中，設計不同大小的疊前5D無噪聲模型數據D(t,cmpx,cmpy,hx,hy)，分別用基於隨機SVD分解的多道奇異譜分析技術(Rsvds方法)，基於Lanczos雙對角分解的多道奇異譜分析技術(Lanczos方法)和本發明實施例的基於無展開隨機QR分解矩陣降秩技術進行重建計算效率比較。每個數據在時間方向含有nt個採樣點，例如nt＝256，採樣間隔dt＝0.004s。在空間方向的大小為N×N×N×N，共包含三個不同斜率的同相軸。可取N＝6,7,…,13，最大迭代次數Niter＝10，秩r＝6，加權因子α＝1.0，迭代停止誤差ε＝10-4，然後用三種方法進行重建。圖8是分別利用現有方法和本發明一實施例的利用隨機QR分解進行矩陣降秩的地震數據重建方法進行重建所耗費計算時間的對比示意圖。如圖8所示，本發明實施例的無展開隨機QR分解降秩重建方法要明顯快於現有的Rsvds方法和Lanczos方法兩種重建方法，並且隨著數據維數N的逐漸增大，三種重建方法的計算時間差異也更加明顯。取N＝10時的5D模型數據D(t,cmpx,cmpy,hx,hy)，然後隨機剔除10％,20％,…,80％的地震道形成不規則缺失道地震數據，分別運用三種方法進行重建並比較重建效果。定義重建質量因子Q，Drecon表示重建數據，Dtrue表示原始真實數據，F表示四維數組的斐波那契範數。圖9是分別利用現有方法和本發明一實施例的利用隨機QR分解進行矩陣降秩的地震數據重建方法進行重建的重建質量因子對比示意圖。如圖9所示，本發明實施例的無展開隨機QR分解降秩重建方法和現有的基於隨機QR分解的多道奇異譜分析方法具有相似的重建質量，但是前者比後者計算速度快，計算效率高。而且從圖9還可以看出現有的基於Lanczos雙對角分解的多道奇異譜分析方法的重建質量明顯低於本發明實施例的方法。

一個實施例中，合成一個大小為256×10×10×10×10的5D無噪聲模型數據，隨機剔除50％的地震道形成不規則缺失道地震數據，然後用本發明實施例的方法進行重建。該數據在時間方向含有256個時間採樣點，在空間方向包含10個CMP面元，10個偏移距，設置最大迭代次數Niter＝10次，秩r＝6，加權因子α＝1.0。圖10是本發明一實施例中原始完整數據的cmpx面元切片示意圖。如圖10所示，共中心點x坐標cmpx＝1，3，5，7，9，共中心點y坐標cmpy＝2，偏移距x坐標hx＝1：10，偏移距y坐標hy＝2。圖11是本發明一實施例中隨機缺失50％地震道數據的cmpx面元切片示意圖。圖12是本發明一實施例中重建數據的cmpx面元切片示意圖。圖13是圖10所示原始完整數據與圖12所示重建數據的cmpx面元切片差剖面，如圖13所示，從差剖面上可以看出缺失道數據的振幅值被有效恢復。

一個實施例中，將本發明實施例的方法應用於某油田陸上區塊疊前5D實際數據插值重建處理。該數據大小為301×15×15×13×13，在時間方向含有301個採樣點，在空間方向含有15個CMP面元，13個偏移距，數據整體缺失89％。設置最大迭代次數Niter＝100次，秩r＝10，加權因子α＝0.8，然後分別運用現有的基於Lanczos雙對角分解的多道奇異譜分析方法和本發明實施例的基於無展開隨機QR分解矩陣降秩重建方法對該共中心點-偏移距域不規則缺失道數據進行重建。圖14是本發明一實施例中原始實際採集的觀測數據的cmpy面元切片示意圖。如圖14所示，共中心點x坐標cmpx＝5，共中心點y坐標cmpy＝5，6，7，8，9，偏移距x坐標hx＝8和偏移距y坐標hy＝1:13。圖15是利用現有基於Lanczos雙對角分解的多道奇異譜分析方法的重建結果的cmpy面元切片示意圖。圖16是利用本發明實施例方法的重建結果的cmpy面元切片示意圖。從圖15和圖16的矩形窗中可以看出，在相同的重建條件下本發明實施例的方法在重建後缺失道數據得到有效恢復並且同相軸的振幅能量也更強光滑連續。

本發明實施例的基於隨機QR分解矩陣降秩的地震數據重建方法，通過僅存儲四重塊Toeplitz矩陣的部分行列元素來表示整個四重塊Toeplitz矩陣能夠降低地震數據在重建過程中的數據存儲量；利用隨機QR降秩分解算法可以實現對部分行列元素表示的四重塊Toeplitz矩陣進行降秩處理，以此可以減小降秩計算量，提高地震數據重建效率；利用無展開求平均算法對降秩處理後的四重塊Toeplitz矩陣的對角線元素求平均值，無需計算並存儲降秩處理後的四重塊Toeplitz矩陣的具體形式，無需存儲該大型低秩四重塊Toeplitz矩陣，可以有效減少數據存儲量，提高地震數據重建效率。本發明的方法能夠減少數據儲存量和計算量且不會降低重建質量。本發明實施例的方法不僅能夠用於五維地震數據的重建，還可以很好地用於N維地震數據的重建，N≥2。

在本說明書的描述中，參考術語「一個實施例」、「一個具體實施例」、「一些實施例」、「例如」、「示例」、「具體示例」、或「一些示例」等的描述意指結合該實施例或示例描述的具體特徵、結構、材料或者特點包含於本發明的至少一個實施例或示例中。在本說明書中，對上述術語的示意性表述不一定指的是相同的實施例或示例。而且，描述的具體特徵、結構、材料或者特點可以在任何的一個或多個實施例或示例中以合適的方式結合。各實施例中涉及的步驟順序用於示意性說明本發明的實施，其中的步驟順序不作限定，可根據需要作適當調整。

本領域內的技術人員應明白，本發明的實施例可提供為方法、系統、或電腦程式產品。因此，本發明可採用完全硬體實施例、完全軟體實施例、或結合軟體和硬體方面的實施例的形式。而且，本發明可採用在一個或多個其中包含有計算機可用程序代碼的計算機可用存儲介質(包括但不限於磁碟存儲器、CD-ROM、光學存儲器等)上實施的電腦程式產品的形式。

本發明是參照根據本發明實施例的方法、設備(系統)、和電腦程式產品的流程圖和/或方框圖來描述的。應理解可由電腦程式指令實現流程圖和/或方框圖中的每一流程和/或方框、以及流程圖和/或方框圖中的流程和/或方框的結合。可提供這些電腦程式指令到通用計算機、專用計算機、嵌入式處理機或其他可編程數據處理設備的處理器以產生一個機器，使得通過計算機或其他可編程數據處理設備的處理器執行的指令產生用於實現在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能的裝置。

這些電腦程式指令也可存儲在能引導計算機或其他可編程數據處理設備以特定方式工作的計算機可讀存儲器中，使得存儲在該計算機可讀存儲器中的指令產生包括指令裝置的製造品，該指令裝置實現在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能。

這些電腦程式指令也可裝載到計算機或其他可編程數據處理設備上，使得在計算機或其他可編程設備上執行一系列操作步驟以產生計算機實現的處理，從而在計算機或其他可編程設備上執行的指令提供用於實現在流程圖一個流程或多個流程和/或方框圖一個方框或多個方框中指定的功能的步驟。

以上所述的具體實施例，對本發明的目的、技術方案和有益效果進行了進一步詳細說明，所應理解的是，以上所述僅為本發明的具體實施例而已，並不用於限定本發明的保護範圍，凡在本發明的精神和原則之內，所做的任何修改、等同替換、改進等，均應包含在本發明的保護範圍之內。