Skip to content

Latest commit

 

History

History
218 lines (163 loc) · 11.7 KB

File metadata and controls

218 lines (163 loc) · 11.7 KB

pyTECTOR:方法細節

準則、Angelier 的參數化,以及由此而來的各項後果,比 README 更深入。 此處內容或出自已發表論文,或量測自原程式自身的輸出; 凡引用數字者,產生該數字的腳本都在 research/

回到 README


PSIDIR 什麼時候會把軸重新貼標籤

INVDIR 的軸來自未正規化的張量,因此它為 σ₁/σ₂/σ₃ 指派的標籤不保證正確。 PSIDIR 的職責就是檢查並在必要時修正;但它並非每次都會實際改動任何東西, 而一個 run 屬於哪一種情況,決定了最終應該把哪一組數字當成答案。

在凍結的 INVDIR 軸上,正規化張量的三個特徵值為 cos ψ、cos(ψ+120°)、cos(ψ+240°), 其中 ψ 取 PSIDIR 掃描收斂到的值。這三個值只有在 ψ 落在整圈的最後 60° (300° 到 360°)時才依大小排序,也就是與 INVDIR 原本設定的 σ₁ > σ₂ > σ₃ 順序一致。 其餘每個 60° 區間各對應另外五種排列之一。

  • ψ 落在 [300°, 360°) → AXES OK !。INVDIR 的 σ₁/σ₂/σ₃ 標籤維持不變。 Φ 一律會換成 PSIDIR 計算的值,但該 Φ 所屬的軸仍是 INVDIR 原本報告的那一組, 因此無論將此解歸給 INVDIR 或 PSIDIR,讀出來的都是同一組數字。
  • ψ 落在其他區間 → PERMUTATION。σ₁/σ₂/σ₃ 被重新指派到與 INVDIR 標籤 不同的凍結方向上。此時 INVDIR 印出的 Φ 與最終 PSIDIR 的 Φ 描述的是不同的 軸指派,而非同一組軸上的細微修正;這正是必須改讀 PSIDIR 數值、 不可再引用 INVDIR 數值的情況。pyTECTOR 遇到時一律標示:back-tilt 視窗的 摘要會印出 PSIDIR ...: sigma1/2/3 are NOT INVDIR's labels, 底層旗標為 permutation=Truepsidir_flag='PERMUTATION' (見 pytector.invdir.axis_order)。

以 archive 驗證:在 pyTECTOR 重現至 σ₁ 誤差 3° 以內的 56 個 run 中, 上述「300-360°」規則 56 次全部正確預測了檔案記錄的旗標。 其餘 6 個未命中的 run,本身連 INVDIR 解都未能重現, 因此那 6 筆檢驗的是重現程度,而非這條規則。

Angelier 本人的說明(Appendix IV)指出,這一步用於修正「資料方位變化不足」時, 未正規化的 INVDIR 階段會產生的「人為 σ₁/σ₃ 對調」,也就是斷層資料本身 對應力張量的約束不足、未正規化形式無法自行將軸置於正確位置的那些站點。 因此出現 PERMUTATION 並不代表解有問題,而正是 PSIDIR 被設計出來要處理的情形。

S4MINpytector.modern,代碼 S4MN)為同一個 S₄ 的精確最小值。 其採特徵分解參數化,故 λ 依定義恆為 √3/2,無須迭代,且搜尋為全域搜尋。 它在全部 92 個 archive 站均取得更低的 S₄,無一例外 (L12 為 0.2360 對 0.3018;0406-7 為 7.3201 對 7.6198)。 換言之,原程式並未達到其自身準則的最小值,因為 λ 停在未收斂之處。

λ 的定義:Angelier 的原始說明

以下內容並非逆向推測,而是他 1990 年論文 Section 4 與 Appendix IV 明確載明的。

λ 為該張量所能產生的最大剪應力。準則同時要求預測剪應力方向與實測滑動一致, 且其大小足以克服摩擦,λ 即為其目標值。

正規化的 A16 張量最大剪應力為常數,式 (14) 的張量則不然。原因他闡述得相當清楚: 對角線項帶有 ψ、非對角線項則無,因此轉動應力軸會改變應力大小。 他的結論是「軸的旋轉與應力大小在解析上不獨立」(1990, Section 4), 並指出若非如此,λ 便僅是一個常數,完全無須調整。

他的處理方式即為該迭代:執行數趟,每趟將 λ 更新為前一趟解的最大剪應力。 (NO k) 計數的正是此項。

Angelier 亦明確指出 A16 形式更為理想:採用該形式時 λ 的調整完全不必要、 λ 恆等於 √3/2;他未採用的原因是無法對該式求出解析解,同時註明 「沒有理由認為不可能」(Appendix IV)。 PSIDIR 這最後一步即是在收尾階段套用一次 A16 形式,他自述為「為求穩妥而增設」, 用以修正未正規化形式在資料方位變化不足時所產生的 σ₁/σ₃ 人為對調。

因此 S4MIN 並非將 Angelier 的方法「現代化」,而正是他描述過、 意圖採用但當年無法實現的那個形式,只是改以數值方法達成。 兩種跑法的差距源自 1990 年的解析限制,而非地質判斷上的分歧。

λ 迭代不會收斂,提早停止即為其設計要點

Angelier 僅提及「數趟連續的決定」,未指明趟數,亦未主張收斂。以 archive 驗證: 92 站中有 72 站會發散。 0406-7 即為其中之一:若持續執行, λ 由 0.866 → 1.009 → 1.115 → … 至第 200 趟達 1.1 × 10⁹, S₄ 則由高於最小值 4 % 惡化至 78 %。L12 則會收斂,停在 λ = 2.2404、高於最小值 1.5 %。

其機制為正回饋:λ 增大即要求更大的剪應力, 求解器的因應方式是撐大未正規化的張量,而撐大後的張量最大剪應力又更大, 成為下一輪的 λ。

因此採使用者自選趟數而非「迭代至收斂」並非簡化處理; 在多數站點上,這是唯一能使結果保持有限的做法。 archive 的紀錄亦印證此點:62 站用 NO 1、25 站用 NO 2、五站用 3 至 5。


從 INFO1 判讀收斂程度

INFO1 印出三個數字,容易誤讀其中之一:

SOLUTION INVDIR (NO 1)  LAMBDA= 0.68     <- INVDIR 實際採用的 λ
SOLUTION PSIDIR         AXES OK !
LAMBDA= 0.87            TAUMAX= 0.80     <- PSIDIR:λ 依定義即為 √3/2,
                                            故所有檔案皆印出 0.87

TAUMAX 亦非 √3/2。對特徵值為 cos(ψ + k·2π/3) 的正規化張量,其最大剪應力為

taumax = 3 / (4·√(Φ² − Φ + 1))

其值自 Φ = 0 或 1 時的 0.75 至 Φ = 0.5 時的 0.866。此式與 87 個 archive run 吻合至 0.005 以內,同時也驗證了整體理解無誤。

唯有第一個數字上升至與第三個數字相等,才代表收斂。 在具備這兩個數字的 88 個 run 中,差距中位數為 0.160,且無任何一筆落在 0.02 以內:

TAUMAX − LAMBDA 站數
≤ 0.02(收斂) 0
0.02 到 0.05 4
0.05 到 0.15 31
超過 0.15 53

此處並非批評當年的執行方式:迭代在多數站點上會發散,停在 NO 1 或 NO 2 是正確的判斷。 其意義在於記錄下來的 λ 是一個停止點,而非一個解, 而 archive LAMBDA 這項功能存在的目的即為重現該停止點。

由此可知 INVDIR 與 S4MIN 的差距並非一個固定數字, 而取決於站點以及當年所採用的趟數:

n ≥ 7(55 站) n ≥ 15(10 站)
受約束的軸 中位 8.9°、p90 20.8° 中位 4.8°、最大 12.5°
|ΔΦ| 中位 0.074 中位 0.065
S₄ 高於最小值 中位 27 % 中位 6 %

0406-7 的 4 % 屬於表現較佳的一端。 須注意 S₄ 的百分比在 n 較小時是不可靠的指標:張量有四個未知數, 因此僅有四至五筆資料時全域最小值會趨近於零,任何與之相比的比值都會失去意義。 應據以判讀的是軸與軸之間的夾角。

兩者的實際差距:僅檢視資料真正約束住的那根軸(Φ < 0.5 看 σ₁、Φ > 0.5 看 σ₃), 在 n ≥ 7 的 55 站中,中位數 8.8°、p90 為 16.7°、最大 28.4°。 n ≥ 15 時收斂至中位數 4.8°、最大 12.5°,可見此分歧主要為樣本數效應, 而非方法本身的缺陷。 對照基準:INVD 自身在無雜訊合成資料上的偏差約 4°, 而 Angelier 引用的野外擦痕觀測誤差為 ±5-15°。 兩種跑法的差距落在此雜訊水準之內。

實務建議:以 INVDIR 為主,以維持與既有 run 及 TENSOR 系文獻的一致性; 並以 S4MIN 作為穩健性檢驗併陳。


軸的位移

將實測的軸經同一旋轉轉換後,繪於回轉後的圖上,以空心圈表示,並以虛線弧連接至星形。 理由在於:單獨檢視一張回轉後的圖,無法看出「傾轉」本身, 軸僅是位於不同位置,圖面上並無任何資訊指出其原始方位。

空心圈與星形不必然重合,且是否重合取決於所用方法。

  • S4MIN 為精確等變。 S₄ 具旋轉不變性,因此其最小值隨資料一同旋轉。 由此可推論:回轉不可能改變 S4MIN 的 Φ 與 S₄, 傾轉檢驗的全部內容即在於「軸最終落於何處」。
  • INVDIR 非等變。 式 (14) 將張量的對角線固定於地理座標下的 cos ψ、cos(ψ+2π/3)、cos(ψ+4π/3)。資料一經旋轉,該四參數族即成為另一個族, 對回轉後的資料重新計算,所搜尋的是不同的解空間。

此為 Angelier 方法本身的性質,而非本重建版的缺陷。 以原程式自身執行的十四組 archive 回轉配對驗證(經 Kabsch 擬合 < 2° 確認為同一批資料): 將母站的軸經旋轉轉換後,與檔案中回轉版的計算結果相比, 中位差為 σ₁ 10.4°、σ₂ 24.3°、σ₃ 23.6°,最大達 88.7°。 差距最大的數筆多落在 Φ 接近 0 或 1、兩根軸近乎簡併的站點,但並非全部如此: 0214-5(13 筆、Φ 0.46 → 0.72)σ₁ 仍差 19.8°、σ₃ 差 22.8°。

對研究撰寫的意義:若僅採用 INVDIR,「回轉前後軸發生變化」有一部分反映的是 參數化的性質而非地質意義,因此「回轉後 σ₁₂₃ 回到水平/垂直」本身不足以作為證據。 較嚴謹的判準應採 S4MIN(其軸可證明僅為純旋轉),INVDIR 則保留作為與舊有 run 的連續性對照。視窗會同時列出兩者及其差距。相關測試見 tests/test_backtilt.py

「轉回水平」並非必然正確

該做法預設斷層形成早於傾轉。若斷層是在傾轉過程中形成的, 則僅有部分傾轉發生於其後,將全部傾轉轉回即構成過度回轉, 會得到一個從未存在過的應力張量。

Tilt test 將旋轉自 0 掃描至 125 %,每一步皆執行反演,並繪出兩項診斷量:

診斷量 說明
平均 ANG、RUP、S₄ 單一張量解釋資料的程度。最佳值落在 100 % 附近,表示斷層形成早於傾轉
Andersonian 失配 90° 減去最陡那根軸的傾沒,故 0 代表一軸垂直、兩軸水平。同時判定應力體制:σ₁ 垂直為正斷、σ₂ 為平移、σ₃ 為逆衝

最佳解落在遠低於完全回轉之處,即為同傾轉斷層的特徵,程式會明確標示。 若兩項判準的差異超過旋轉量的 20 %,程式亦會提示, 因為該分歧須先釐清,兩項判準才具參考價值。兩者皆非證明,而是診斷指標。

實例:0404-4C-2 對一個虛構的參考面,兩項判準均隨回轉而變差, 最佳值落在 10 %,Andersonian 失配自 40.8° 升至 49.2°。 將該面轉回水平,會在無警示的情況下得到比原始資料更差的結果。


效能

ψ 掃描原本每趟呼叫純量常式 4000 次,並疊加 120 步三分搜尋, 而 archive LAMBDA 又將整套流程重跑約九十次。實際上其中幾乎沒有一步與 ψ 相關: 僅對角線帶有 ψ,且為線性關係,因此 3×3 normal matrix 為資料集的常數, 每個 ψ 僅需對同一分解執行一次回代。PSIDIR 與 S4MIN 的起點掃描亦以特徵座標系 做同樣的化簡。

0406-7(29 筆)實測結果:

最佳化前 最佳化後
INVDIR 2 趟 0.975 s 0.004 s
S4MIN 400 起點 0.208 s 0.052 s
archive LAMBDA 6.775 s 0.063 s

與被取代的常式差異為 1e-13。