Angelier 古應力反演的 Python 重建版,對應 TENSOR 5.45 (jan91)
依已發表的方法重建,並以原程式自身的輸出檔案驗證
Angelier 的直接反演法從一組實測斷層面與擦痕線理,解出最能解釋它們的簡化應力張量,
也就是三個主應力方向與形狀比 Φ。他的程式 Tensor.exe 自 1991 年起執行這項計算,
大量已發表的古應力研究建立在其結果之上。
pyTECTOR 執行同一套運算、讀寫同一批檔案、繪製同一種圖, 並補上原程式未提供的部分:回轉(back-tilting)、傾轉檢驗, 以及將同一準則精確最小化的第二種跑法, 用以評估一個結果有多少來自方法本身、多少來自資料。
專案名稱取自 TECTOR,即 Angelier 為這套程式的構造方位資料庫所取的名稱, 見於它產生的每一份 INFO1。
| 兩種跑法並列 | 回轉與傾轉檢驗 |
|---|---|
![]() |
![]() |
以上圖片均由本 repo 所含的公開 fixture
tests/fixtures/L12-2/產生。 該站為合成資料而非野外資料,因此上述各圖皆可自行重製。
重建方式。 演算法已完整發表於 Angelier (1984, 1990),因此本專案依論文實作, 未對原始的 16 位元執行檔進行反組譯。
驗證。 以原程式產生的 92 個 run 作為回歸測試集。前向量 (SIGMN、TAU、TAUST、RUP、ANG)逐筆吻合至檔案本身的精度; 依各站記錄的 pass 數與 LAMBDA 重新執行反演後, 90 個可比對的站中有 85 站的三軸角度差在 3° 以內。
INVDIR 與 S4MIN 兩種跑法。 前者為 Angelier 原方法、亦即原程式的跑法; 後者為同一準則的精確最小值。兩者皆不應視為「真實應力」: υ 準則本身帶有系統性偏差,即使餵入零雜訊的理想合成資料, 結果仍與真實張量相差約 4°。並列呈現的目的在於顯示差異所在,而非在其中擇一。
回轉與傾轉檢驗。 此為原程式未提供的功能。旋轉角度以拉桿調整, σ₁、σ₂、σ₃ 即時重算。就 INVDIR 而言,回轉前後主軸的差距反映的是方法本身的性質 而非地質意義(14 組 archive 配對,實測 σ₁ 差異中位數為 10.4°)。
影響力診斷。 擬合殘差大的資料與實際主導結果的資料未必是同一批。 程式對每一筆資料執行 leave-one-out 重新反演,輸出剔除後殘差 ANG* 與 RUP*, 並將「全部資料」與「剔除該筆後」兩組結果並列寫入匯出的 INFO1。
Survey 視窗,一次處理整個 archive。 指向一個放 TENSOR run 的資料夾, 每一個 run 都會列進表格(三軸、Φ、ANG、RUP),旁邊是把應力軸畫在實測位置上的地圖, 以及各變形期的玫瑰圖。四個欄位可編輯,正是四個算不出來的欄位: 分期、斷層型式,以及經緯度。
原格式輸出。 INFO1 與 MOHR1 與原始檔案逐位元組相同; HPGL 匯出係重播原本的繪圖程序,而非另一套獨立實作。 HPGL 是向量格式,投影網在繪圖軟體裡是可編輯的曲線與文字而不是一張圖片, 見把圖帶進排版。 Survey 則匯出總表、斷層資料、地圖點位, 以及線幾何形式的應力軸 GeoJSON,GIS 打開就是一張應力方向圖。
Session 存檔。 全部工作狀態存為單一 JSON 檔。檔中僅保存張量, 其餘數值於載入時重新計算,因此存檔中的 Φ 不可能與存檔中的張量互相矛盾。
演算法已完整發表,因此本專案照論文重建,而非反譯 16 位元的原始執行檔: 編譯已丟棄名稱、型別與結構,分段定址也讓指標無法解析。工作集中於研讀 Angelier (1984, 1990, 1994),並量測原程式自身的輸出以補足論文未載明的部分。
名稱同樣出自 Angelier。他的論文從未為程式命名,但執行檔在它寫出的每一份
INFO1 上自報其名,而橫幅上兩個名稱之一即為 TECTOR,亦即它的構造方位
資料庫。「TENSOR」無法使用:在本領域它現指 Delvaux 的 Win-Tensor,
而 PyPI 上的 pytensor 是 PyMC 的張量庫。
兩個問題的完整說明、橫幅原文與參考文獻: docs/background.zh.md。
這是強烈建議,各平台皆然。pyTECTOR 在 python.org 的一般安裝上也能正常運作, 但 conda 可以排除掉大部分安裝失敗的三個原因:
- numpy、scipy、matplotlib 本來就有,或是以預先編譯好的二進位檔提供, 不需要現場編譯,網路慢或被過濾時也比較不會裝到一半就卡住;
- 在 Windows 上可以避開 Microsoft Store 的
python樁:它不是真正的直譯器, PyQt5 在它底下不能用。這是安裝失敗最常見的單一原因; - 在 Apple Silicon 的 Mac 上,Qt 會以 arm64 二進位檔提供, 而不是讓 pip 從原始碼編譯。
Miniconda 是兩者中較小的,而且完全夠用;已經裝了 Anaconda 的話同樣沒問題。 若是為機構而非個人安裝,請留意 Anaconda Inc. 的套件庫對較大型組織訂有授權條款, 而 Miniconda 搭配 conda-forge channel 不會有這個問題。
下載頁:https://www.anaconda.com/download/success。 安裝程式自己的預設選項就是對的,不需要系統管理員權限, 也不需要勾選「Add to PATH」。
下載整個 repo(Code → Download ZIP)、解壓縮,雙擊 install.bat。
它會一次把該裝的都裝好:
| 步驟 | |
|---|---|
| 尋找直譯器 | 先找 Anaconda 與 Miniconda,再找 PATH,最後找 py 啟動器。Store 樁會依路徑排除 |
| 安裝相依套件 | numpy、scipy、matplotlib、PyQt5;pip 裝不起來的部分改用 conda-forge |
| 確認四個都能 import | 而不是只看 pip 的離開碼 |
| 記下用的是哪一個 Python | 寫進 python-path.txt,讓 pyTECTOR.bat 啟動的是同一個,而不是 PATH 上的另一個 |
| 編譯檢查程式 | 並在桌面建立 pyTECTOR 捷徑 |
可以重複執行,也不需要系統管理員權限。
若機器上完全沒有 Python,安裝腳本會詢問是否代為從 repo.anaconda.com
下載並安裝 Miniconda:約 80 MB、兩三分鐘、只裝給目前這個使用者。
回答 n 則改為開啟下載頁。它會盡量裝到 C:\Miniconda3,
不行才退回使用者資料夾,因為 conda 與 Qt 在非純 ASCII 的家目錄底下都會出問題,
而中文或日文的 Windows 帳號正是這種情況。
函式庫本身沒有任何 Windows 專屬程式碼,因此反演、檔案讀取與各項匯出均可直接執行。
install.command 是 install.bat 的對應版本,步驟相同,
同樣會提供 Miniconda 安裝選項,在 Finder 中可直接雙擊。
若雙擊沒反應,表示它掉了執行位元:
chmod +x install.command
./install.command
介面本身尚未在 macOS 上實際測試:字體堆疊已納入 macOS 與 Linux 的字型 以維持定寬表格的對齊,若發現顯示異常歡迎回報。
不想跑安裝腳本的話:Python 3.8 或更新,加上那四個套件。
pip install .(或 pip install git+https://github.com/FormosaRes/pyTECTOR)
亦可,安裝後會提供 pytector 指令。
python -m pip install numpy scipy matplotlib PyQt5
在 Apple Silicon 的 Mac 上,PyQt5 需使用具備 arm64 wheel 的版本
(5.15.10 或更新);若 pip 開始從原始碼編譯 Qt,請改以 conda 安裝:
conda install -c conda-forge pyqt。
pyTECTOR.bat 桌面介面(Windows)
./pyTECTOR.command 桌面介面(macOS、Linux)
python demo_report.py [站檔] 反演一個舊站,印出 INFO1 + MOHR1
python run_batch.py [根目錄] [out.csv] 對整棵資料夾跑兩種方法
python make_survey.py [根目錄] [輸出夾] 產生表格、地圖點位與各期玫瑰圖
兩個啟動器都會使用安裝腳本記下的那個直譯器。
若要自己指定,設定 PYTECTOR_PYTHON 即可,它的優先序高於其他一切。
逐項功能的完整操作說明見 docs/manual.zh.md。
輸入分為四個欄位,輸入後即顯示於投影網上:
CS - 122 - 87W - 124
| | | |
| | | +-- rake +象限(62N),或直接給 trend(124)
| | +-------- 傾角+象限
| +-------------- 走向
+------------------- 信心度 C/P/S + 運動方式 I/N/S/D
亦可直接開啟舊有的 run:指向那個沒有副檔名、檔名等於站名的資料檔(例如 L12),
斷層資料與 INFO1 會一併載入。
- σ = T·n(式 3);τ = σ − (n·σ)n(式 4-5)
- υ² = λ² + |τ|² − 2λ(s·σ)(式 A1),最小化 S₄ = Συ²(式 13)
- RUP = 100·|υ|/(√3/2),範圍 0-200 %;ANG = s 與 τ 的夾角,0-180°
- 張量正規化採特徵值平方和 = 3/2(式 A16),並非固定 σ₁ − σ₃。 兩種定義僅在 Φ = 0.5 時相同,而此準則不具尺度不變性,故此差異會改變結果。
此準則本身帶有系統性偏差,並非任何程式的錯誤。 υ 同時要求「預測剪應力方向與實測滑動一致」與「剪應力大小接近 λ」, 因此即使餵入零雜訊的理想 Bott 合成資料,結果仍與真實張量相差約 4°, 並系統性地偏向使斷層承受高剪應力的方位。 純角度準則 F2 在同一批資料上的誤差為 0.00°。 這即是 TectonicsFP 等軟體與 Angelier 計算結果不一致的根本原因。
兩者以其實質內容命名,避免使用暗示優劣的代號。基於上述理由,兩者皆非「真應力」。
INVDIR(pytector.invdir,代碼 INVD)為 Angelier 原方法、原程式的跑法,
採用他自己的 (α, β, γ, ψ) 參數化(式 14 / A2)。此張量未經正規化,
其最大剪應力隨解的移動而改變,這正是 λ 必須逐趟迭代的原因,
也是 INFO1 所印出的 LAMBDA 小於 √3/2 的原因。流程分兩段:
- INVDIR 第 k 趟:在當前 λ 下最小化 S₄,隨後將 λ 更新為該解的 taumax。
INFO1 印出的
(NO k)為迭代趟數,而非「在兩個解之間擇一」。 - PSIDIR:將軸凍結於 INVDIR 的結果,改用正規化的 A16 形式、λ = √3/2, 對 ψ 整圈重新最小化。此步驟決定 Φ,並修正前一階段可能產生的 σ₁/σ₃ 人為對調。
兩項實作要點:固定 ψ 時 υ² 對 (α, β, γ) 為二次式,內層即精確的 3×3 線性求解, 無須轉抄附錄中的多項式(現有掃描檔的附錄辨識度不足,此處改以數值方式重新生成)。 另外 ψ 必須掃描整圈:若限制於 [0, π/3] 會得到 Φ = 0 的錯誤結果, 因為最小值常落在 ψ ≈ 336-353°。
S4MIN(pytector.modern,代碼 S4MN)為同一 S₄ 的精確最小值,
採特徵分解參數化,λ 依構造恆為 √3/2,無須調整迴圈,且搜尋為全域。
它在 92 個 archive 站上無一例外皆達到更低的 S₄:
| 站 | INVDIR S₄ | S4MIN S₄ |
|---|---|---|
| L12 | 0.3018 | 0.2360 |
| 0406-7 | 7.6198 | 7.3201 |
換言之,原程式並未達到其自身準則的最小值,原因是 λ 在收斂前即停止。
兩種跑法的完整處理、PSIDIR 何時重貼軸標籤、λ 究竟是什麼、迭代為何發散、 如何從 INFO1 判讀收斂程度,以及上面那段回轉警告背後的等變性實測數據, 全部在 docs/method.zh.md。
TENSOR 其餘四個方法(R4DT、R4DS、R2DT、R2DS)的準則罰的是什麼, 以及 S4MIN 為何與 INVD 同組而非與它們同組,在 docs/solvers.zh.md。
檔案格式、HPGL 畫風判準與專案結構在 docs/formats.zh.md。
回轉功能有獨立的視窗,由工具列開啟。主視窗僅呈現實測資料、不做任何旋轉, 因此該處的投影網無須額外標註目前所處的方位。
回轉視窗會將資料轉至指定方位,並對前後兩個狀態各執行一次反演,左右並列呈現: 左側為 as measured、右側為 back-tilted,下方列出兩者的數值。設定旋轉的方式有三種:
| 方式 | 輸入 | 作用 |
|---|---|---|
| 參考面 | 走向 / 傾角,或以其 pole 給定 trend / plunge | 將該面轉回水平所需的旋轉 |
| 旋轉軸 | trend / plunge / 角度 | 直接套用,依右手定則 |
| 部分回轉 | 0 到 125 % | 上述任一方式的任意比例 |
斷層面法向與滑動向量會一併旋轉,因此 rake 與運動感亦隨之改變。 此慣例並非推測,而是對 archive 驗證過的:解出原站與其回轉版之間的旋轉 (對法向施以 Kabsch 演算法),七組全部重現至 2° 以內。
參考面以虛線大圓及其 pole 繪出,並隨資料一同旋轉, 因此回轉是否正確可直接目視判斷:虛線圓會壓平至基準圓上,pole 則移至圓心。
旋轉角度並非計算而得。 該角度無解析解,須以試誤方式檢視結果,
這也是 archive 資料夾名稱直接記錄當年試用值的原因,例如 (backtilted 020 -20)。
本程式的作用僅在於加速試誤過程,並明確標示當前套用的旋轉;
參考面的選擇與旋轉角度的判斷仍由使用者決定。
單次反演回答的是一個站。一項研究真正要問的是一整批站說了什麼, 而在 0.4.0 之前,這個問題只能從命令列去問。工具列的 Survey 打開它。
指向一份成果表(一列一個解,.csv 或 .xlsx),每一列都會出現在表格裡:
站名、n、採用了哪個解法、有沒有回轉、Φ、ANG、RUP,以及三個軸。
輸入是一張表而不是一個放 run 的資料夾,這正是重點。資料夾裡放的是程式印出來 的東西,在整理過的資料上,那不等於研究採用的東西:玉里帶 47 站裡有 14 站改用 S4MIN 重解、5 站做了回轉,而這些事在資料夾的任何一個檔案裡都不存在。 去讀資料夾拿到的是印出來的那個答案,而且沒有任何東西會警告你。 表則把採用的三軸與產生它的決定並排放著;這個視窗不開任何測站檔, 所以它裡面的東西不可能跟測站檔打架。
必要欄位:site、s1_trend、s1_plunge、s3_trend、s3_plunge。
其餘(no、river、longitude、latitude、stage、type、solver、
backtilt、n、phi、ANG、RUP)都是選用,原樣帶過。
四個欄位可編輯,正是四個算不出來的欄位:
| 欄位 | 為什麼歸你決定 |
|---|---|
| phase | 一個解屬於哪一期是判斷。程式不會替你猜;沒有指定時玫瑰圖區塊會直接說明,而不是自行分組 |
| type | normal/thrust/strike-slip,決定該期讀哪一個軸 |
| longitude, latitude | 測站位置 |
這四欄有底色,讓表格讀起來像表單而非報表;其餘欄位是算出來的,維持原樣。 指定分期後,玫瑰圖與地圖會隨著你打字即時重畫。
要用 Excel 就開 xlsx,不要開 csv。 csv 沒有儲存格型別,Excel 會把 9-1 這種
no 重新讀成日期,再存成 cp950+TAB;同一張表已經因此掉過兩次那個欄位,而且不保證
救得回來,因為 10-1 和 1-10 轉出來是同一個日期。兩種格式都讀得進來,
所以 xlsx 拿來改、csv 留給腳本讀。被 Excel 動過的 csv 會直接被擋下並說明是 Excel
造成的,不會丟出解碼錯誤。讀 xlsx 需要 openpyxl,列為選用套件:
pip install .[excel]。
程式旁邊的 py_data/ 是放成果表與測站資料的地方,
整個資料夾都在 gitignore 裡,因為那是野外資料。
命令列的路徑仍然保留,做同樣的事:
python make_survey.py [根目錄] [輸出夾] --stages phases.csv --coords coords.csv
底下鋪 OpenStreetMap,第一次抓取後會快取。上面每一站沿其應力軸畫一個符號, 依分期上色。滾輪縮放、拖曳平移。
符號可以是穿過測站的線段,也可以是慣用的古應力箭頭: σ₁ 向內因為壓縮是推,σ₃ 向外因為張力是拉。 兩端都畫箭頭,理由和線段穿過測站而不是從測站射出去一樣: 軸沒有單一指向,只畫一個箭頭等於宣稱資料裡沒有的東西。 傾伏超過 45° 的軸不畫,因為那種軸的 trend 近乎任意。
圖層控制收納底圖、GeoTIFF 疊圖、符號選項,以及各期一個勾選框, 可以單獨檢視某一期 —— 25 站的期別畫在 3 站的期別上面會把它整個蓋掉。
trajectories 勾起來會畫 Lee 與 Angelier 的 LISSAGE 平滑應力軌跡, 一期一組線、用該期的顏色,畫在地圖上而不是另開一張圖 —— 平滑場的重點就是它在哪裡轉向,而一張軸標「距中心幾公里」的圖沒辦法跟河流對照。 軸是整期共用一根、由該期的斷層型態決定,所以伸張期畫的是 σ3, 期內那一兩個逆衝站不會把 σ1 灌進伸張場裡。
旁邊是搜尋半徑 radius,相信那些線之前要先看懂它。
它幾乎不影響方向:4 km 與 30 km 的中位殘差差不到一度,
散布在測站本身而不在內插。它完全決定畫到哪裡:
半徑內沒有任何測站的點根本不會被畫。這個值原本由各期自己的散布推算,
方向是反的 —— 四個散得很開的測站拿到 22 km 蓋滿全區,
六個分成相距 40 km 兩坨的測站只拿到 14 km,場在中間就斷了。
五期共用一個半徑、並寫在狀態列上,各期圖才能互相比較。
要用回舊規則就選 auto (per phase)。
GeoTIFF 以 Pillow 讀取,而非引入 rasterio 或 GDAL:那兩者都很大、 在 Windows 上難裝,而且會破壞安裝器「只裝四個套件」的承諾。 影像須為北方朝上,投影須為 EPSG 4326、3857、3826(TWD97 TM2) 或 32651(UTM 51N),包含投影寫在檔案 GeoKey 裡而非用代碼命名的常見情形。 其餘一律具名拒絕並提供手動選擇,而不是靠猜測擺上去: 一張悄悄偏移幾百公尺的底圖看起來是對的,這正是它比沒有底圖更糟的原因。
玫瑰圖採軸向而非方向統計,因為應力軸沒有箭頭:020 與 200 是同一條線, 故採倍角法,使兩端互相加強而非互相抵消。 唯有淺傾的軸其 trend 才視為方向;較陡的軸會被剔除, 且剔除了幾個會直接印在圖上而非默默略過。
一期該讀哪一個軸,取決於斷層型式,而非比較哪個軸的可用數較多。
比數量看似有原則,其實不然:正斷層情境下 σ₂ 與 σ₃ 都水平,
兩者可用數相同,勝負落到 R 值。
在一個 25 站的期別上,這讓 σ₂ 以 0.03 的差距勝過 σ₃,
報出來的方向與真正的張力方向差了 90°。平移期有完全相同的陷阱。
完整說明見 pytector/rose.py。
Save PNG 存出來的是圖片,Save HPGL 存出來的是向量。 要放進論文的圖用後者:每一個圓、刻度、標註都是可以重新上色、改字級、 自由縮放的物件,放大不會糊。
HPGL 是惠普 1970 年代的繪圖機語言。Angelier 的程式寫這個格式, 是因為當年的繪圖機吃這個;它比硬體活得久,正是因為夠簡單 —— 純文字、絕對座標、一行一個筆的指令。 本程式的匯出是重播螢幕上那支繪圖程序,所以畫出來的和存出來的不會走鐘。
CorelDRAW 可以直接匯入:檔案 ▸ 匯入,若檔案被當成未知格式,
在篩選器選 .plt。這是本匯出實際驗證過的路徑,優先用它。
這個格式夠老也夠簡單,其他向量與 CAD 軟體多半也讀得進去; 另外有一個成熟的免費轉檔工具 hp2xx, 可以把 HPGL 轉成 SVG、PostScript 或 PDF, 等於也打通了 Illustrator 與 Inkscape 這條路。 這些路徑我沒有實測過,如果哪一個把圖弄壞了,值得回報。
兩個實務上的提醒。有些匯入器認 .plt 而不認 .hpgl,改副檔名通常就解決了。
另外這種檔案本身不帶頁面尺寸,匯入後第一件事是把它縮放到圖框需要的大小。
十一個測試檔全數通過。有 archive 時即讀取,無 archive 時則 skip,不會 fail。
涵蓋範圍包括:整條流程分別對照原程式輸出與公開 fixture、INFO1 與 MOHR1 的
版面、打字輸入的解析、回轉慣例、等變性結果、影響力診斷、軸向統計,
以及 session 存讀一輪。各測試各自驗證什麼見 tests/。
0406-7(29 筆,傾角 42-89°)是確立演算法的關鍵站點:
前向模型 vs MOHR1 max |SIGMN| 0.001 |TAU| 0.001 |TAUST| 0.001
max |RUP| 0.099 |ANG| 0.113
INVDIR 流程 sigma1 0.047 度 sigma2 0.020 sigma3 0.032
Phi 0.138(檔案 0.138)
平均 ANG 20.898(20.900) 平均 RUP 54.097(54.100)
印出的 LAMBDA 0.682(0.680)
L12(六個近平行、近垂直的面)為簡併站點,重現程度較為寬鬆:軸約 1°、 平均 RUP 在 0.6 % 以內。容差係逐站設定以反映此差異,而非全域放寬。
測試中不會建立任何 Qt 物件:自動化 shell 啟動 QApplication 會彈出平台外掛 錯誤對話框並結束程序。測試改以驗證介面與函式庫之間的契約為主。
測試與推導腳本會讀取原程式的實際輸出。該資料屬未發表的野外資料, 因此路徑不寫入原始碼中:
set PYTECTOR_ARCHIVE=<放 TENSOR run 資料夾的那個目錄> REM Windows
export PYTECTOR_ARCHIVE=<放 TENSOR run 資料夾的那個目錄> # macOS、Linux
未設定時測試會 skip,不會 fail。資料本身不隨本 repo 發布。
本 repo 的程式碼以 MIT 授權釋出,全文見 LICENSE。 若本工具對已發表研究有所助益,歡迎在引用 Angelier 原始論文之外一併引用本 repo, 但並非必要。
以下三項不在本授權範圍內,因為它們並非本專案所能授權的對象:
- 方法本身為 Jacques Angelier 所提出,載於前文所引論文。本專案係依其發表的 論文,並輔以對其程式輸出檔案的量測所完成的獨立重寫,未自原執行檔取用任何 程式碼。
- 參考資料集屬未發表的野外資料,未隨本 repo 發布,亦不在本授權範圍內 (見「參考資料集」一節)。
- 開啟畫面採用 Angelier 繪製的台灣弧陸碰撞塊體圖。他以玉里附近的地震震源
機制作為「將此方法應用於地震學資料」的示範案例(1994, fig. 4.44)。該圖為
已發表的圖件,未納入本 repo,因此新 clone 的版本不會有開啟畫面,
其後的彩蛋亦無法開啟。將
Taiwan Tectonic Map.jpg置於專案根目錄即可恢復兩者, 僅限本機使用。
維護者:Chi-Hsiu Pang。
- 寫出 TENSOR 格式的資料檔(站頭欄位已由原程式實際執行解出, 見 docs/mesure_oracle.md)
- R4DT/R4DS/R2DT/R2DS 等 Angelier 的疊代搜尋法,刻意尚未著手 (TENSOR 自身的說明文件有相關記載,見 docs/mesure_oracle.md)
見 CHANGELOG.md。



