地圖上的各種等高線，是怎麼畫出來的？

簡單概括，畫等高線的步驟：如果有類似 DEM 這樣的柵格地形數據，就可以直接找 DEM 畫等高線。如果是一系列離散的高程數據點，那麼繪製過程就是：定位高程數據 -- 空白區空間插值 -- 勾畫等高線。

地理書上的各種等高線及其它相關的線，確實是經過大量的實地測量之後得到的。進行這些測量的人，大多數是有關部門（比如美國的 USGS、NOAA、Land Survey 之類的，估計是國土資源部？）專門負責搞測繪和地質調查的，而不是一般意義上的地理學家。而且實在不行的話，還有遙感呢……

海底的地形，可以用聲吶等設備來進行測繪。至於降水量、氣壓、氣溫之類的等值線的數據，全國各地的觀測點肯定比你想象的多，所以不用擔心，它們都是靠譜的哈。

不過對於降水量、氣壓、氣溫等數據，再多的觀測網點，確實也不可能覆蓋每一寸土地，難免有空白。如果你想問的是『在數據點之間的空白範圍，等值線都是如何畫出來的』，這裡就涉及到繪製等值線時會用到的給空白區域進行空間插值的演算法。常見的演算法有這麼幾個：

• 三角測量法（TIN）
• 反距離權重法（IDW）
• 克里金法（Kriging）
• 趨勢法（Trend）
• 最小曲率法（Minimum Curvature）

這幾種常見演算法是怎麼運作的？下面直接放以前我講課的課件截圖加以說明好了。

三角測量法 Triangulated Irregular Network - TIN：

核心思想：直截了當地假設兩點間的坡度不變

三角測量法是用得最多的一種方法，因為它在數學上最簡單，也是最易於手畫的一種方法。以地形為例。比如，我們現在有這麼一個區域，區域內有一堆海拔高程數據點（當年做課件的時候瞎編的數據，單位假設為米）：

每一個紅三角就是一個數據點，旁邊的數字是高程。我們現在用三角測量法來繪製這片區域的等高線圖。第一步，就是要建立 TIN 三角網路。也就是說，要把這些點用直線連接起來，組成一個三角形的網格。那麼具體哪些點之間需要連線？這裡有一個 Delaunay 法則：三角形網格中，任何三角形的外接圓，都不會含有別的數據點。也就是說，所有三角形的最小的那個內角都被最大化。

這裡的綠色虛線組成的三角網路就是符合 Delaunay 法則的網格，而淡紫色的那個三角形就是不符合要求的三角形。

有了 TIN 三角形網格，下一步就可以畫等高線了。這裡的主要思路就是根據兩個端點的高程，按比例分割每一條虛線。換句話說，我們需要假設，任何兩個高程點之間的坡度，都是不變的。比如，我決定每 50 米畫一條等高線，那麼，對於 1030 米和 694 米那兩個點之間的虛線，應該這麼分割：

照此方法，對每個三角形的每一條邊都做出劃分，那這個區域就大概會變成這樣（局部）：

好了，然後下一步就是把所有數值相同的點用盡量平滑的曲線給連起來（局部）：

這一步，有的教材說在三角測量法的體系下，應該用直線段來連接高程相同的標記點。那其實是做 3D 地形模型的時候採用的手段。在畫等高線的時候，我還是傾向於用平滑曲線。如果用直線段，做出來的效果會類似下面這張圖：

個人覺得這種鋸齒狀的『等高線』看起來很淡騰，所以我自己從來都是用的平滑曲線。畫好曲線之後，把三角網路的輔助線以及標記線擦除，做好標註，邊緣地區根據情況進行一下適當的調整，一張等高線地圖就畫好了：

上面這張圖是在 PPT 裏手畫的，所以並不太準確，曲線也並不太平滑。如果交給計算機（比如 ArcGIS）自動畫圖，那麼做出來的圖要更好，等高線之間的間隔會更平均，曲線也會更平滑。

以上是 TIN 的方法。其它的幾種演算法各有各的特色，但一般都沒法用手畫了（除非你有最強大腦那種超能力），得讓計算機來畫。其它幾種方法之間相似的是，它們都考慮了空間距離的權重，而非直接假設兩點間的斜率不變。

首先，剩下幾種畫法都涉及到尋找（或指定、定義）某個數據點的相鄰點（Neighboring points）。其實，TIN 的三角網格就是定義相鄰點的一種特殊方法。除此之外，還有區域搜索、半徑搜索、柵格搜索、方向搜索等很多方法，這裡不展開了。

當我們定義好了相鄰點之後，就可以運用其它幾種演算法來畫圖了。假設經過空間搜索之後，1030 和 694 那兩個點仍舊被定義為相鄰點，那剩下幾種演算法應該怎麼畫？

反距離權重法 Inverse Distance Weighting - IDW：

核心理念：近朱者赤，近墨者黑，距離誰近就更像誰。

這個方法是一種確定性插值法，思路是根據地理學第一定律來的：『離得近的事物之間的關聯比較大』。因此，一個數據點對周圍區域的海拔是有一定的影響力的。影響力多大呢？為此你要設定一個參數。這個參數是個冪參數，它用來增強某一個數據點對附近的區域的影響力。比如當這個參數是 2 的時候（比較常用的一個設定，參數越大、影響範圍越大）：

根據設定的權重參數，計算機會模擬出一條 AB 兩點間的海拔曲線。距離 A 點近的地方，海拔也接近 A；距離 A 越遠，A 的影響力越弱，海拔就變化越快，直到進入 B 的控制範圍，被 B 的海拔征服。這樣，就找到了 AB 間各條等高線的畫線位置，剩下的事兒，要麼像前面三角測量法 TIN 那麼手動畫，要麼就扔給計算機吧。

反距離權重法一般在數據點比較集中，且均勻分佈的時候使用，是一種經常用來定位異常值的方法。但是如果數據點太少或者太稀疏，或分佈不均勻，那這種方法就不太好用。如果你看到一幅等高線地圖上過於頻繁地出現下面的這種『牛眼』形狀，那麼這幅圖很可能就是用反距離權重法畫的。

克里金法 Kriging：

核心思想：有理有據的統計插值

克里金法是空間統計學和地球統計學中的常用方法，是一種高斯過程主導的插值法。發明這種方法的人叫丹尼·克里吉，他在南非金山大學讀碩士的時候創造出了這種方法，後來被巴黎高等礦業學校的教授馬瑟隆發揚光大。

克里金法是一種空間自相關的統計模型，它根據一系列相鄰點之間的統計關係（比如協方差）來模擬相鄰點之間空白區域的高程變化。換句話說，它就是把指定範圍內的所有數據點都放一起進行統計，然後根據每個距離範圍內的點之間高程半方差和距離，找出（預測出）這個範圍之內最佳的擬合地表。它被視作是最佳線性無偏差預測（Best Linear Unbiased Estimation，BLUE），因為在各種空間的線性預測裡面，它的偏差是最小的。

這麼說好像太抽象了，會暈菜。舉個栗子：

山坡上有一堆點。其中 A 是我選定的第一個參考點，綠色的點為定義的 A 的相鄰點。首先我找出 A 到每一個相鄰點之間的距離。然後，計算 A 和每一個綠點之間的高程的半方差（Semivariance）。一般來說，用克里金法的時候，數據點會很多（比圖上的多得多），因此我們要分組做。比如第一組是到 A 的距離大於 5 米的點，這裡只有 B 一個，那麼這一組的半方差就是（這裡只是舉例，現實中數據點會很多很多，不會出現這麼一個點的情況）：

SemiVariance_1 = 0.5 * (A 海拔 - B 海拔)^2

第二組是到 A 的距離在 3 米到 5 米之間的點，那麼這一組的半方差就是：

SemiVariance_2 = 0.5 * 1/4 * ((A 海拔 -C 海拔)^2 + (A 海拔 -D 海拔)^2 + (A 海拔 -E 海拔)^2 + (A 海拔 -F 海拔^2))

以此類推。這樣的距離叫做 Lag。Lag 的設定不能太大（也就是說，相鄰點的定義要嚴格），否則統計就會不顯著。

對 A 做完統計，得到各個 Lag 對應的半方差，之後如法炮製，對其它每一個點也都要先找出相鄰點，然後一一求得不同 Lag 的半方差。

假設我們算出了很多組這樣的半方差，每組半方差對應一組距離範圍，我們就可以把我們的結果作為散點，畫在一張 Variogram 圖表上，然後去用計算機擬合。擬合的方法至少有球面、圓弧、高斯和指數這幾種，根據需要自己選擇。比如：

通過這張 Variogram 圖表，我們可以得到距離數據點不同的距離，所對應的權重。根據權重，計算機就能得到兩個數據點之間海拔變化的曲線。比如，運用球面擬合，一般來說會得到形狀類似下面這樣的曲線（僅僅是示意圖，具體曲線根據權重來定）：

此後，再按照前面 TIN 的最後的步驟進行等高線勾繪就行了，或者繼續扔給計算機。

克里金法的計算量很大，一般適合用於數據點很多或者分佈不太均勻的情況（反正都是交給計算機做 666）。實際上，在畫地形圖的時候，克里金法使用得並不比其它幾種方法多。克里金法最受歡迎是在研究水文、礦產或者土壤的時候。比如一個區域的土壤，進行大規模採樣之後，繪製等 pH 值圖，等等。並且，用 ArcGIS 等軟體做 Kriging，最後會得到一個 Standard Error 的總結報告，告訴你等高線圖的靠譜程度。

趨勢法 Trend:

核心思想：嘗試用一張紙去擬合，只求個大概就行...

趨勢法比克里金法容易理解多了。比如我在山坡上有很多數據點，有高有低。然後我們去掉山坡，只剩下一堆點在空間里。我現在拿一張硬紙板，放到這堆點裡。至於怎麼放，有講究：放的時候肯定有的點在紙的上方，有的點在紙的下方，但紙的角度要讓空間里每一個點到硬紙板的距離的均方根最小。也就是說，每個點到紙面的距離進行平方，再取平均數，再開方，這個值要最小。

這是山坡的情況。那如果是個 U 形的峽谷呢？那就不用硬紙板了，用普通的紙，做成一個 U 形去擬合。也就是說，增加了代表紙面的多項式的階數，從一階多項式變成二階多項式。如果有更複雜的地形，還可以繼續增加階數，比如變成三階（轉兩個彎）。

之後呢，再對紙面繪製等高線。

這種方法的好處是便於理解、方便操作。但它的缺點也很多，比如大多數的點其實並不在紙面上，因此等高線無法真實反映每一個數據點的高程，只能反映一個大致的情況。這種方法一般在做初步研究的時候使用，可以提供一個區域的宏觀情況。真正研究的時候，一般還是會採用更為精確的方法。

最小曲率法 Minimum Curvature:

核心思想：嘗試用橡膠膜擬合，一個點也不能少...

前面的 Kriging 也好，TIN 也好，IDW 也好，兩點之間的所有地方，模擬出的高程都在兩點的高程之間。正常情況下，這沒問題。但如果我知道我畫的地方是一片丘陵，高低起伏不定，也許兩點間有地方會高於或低於端點的高程範圍，怎樣才能更好地模擬這種地形呢？這就要用最小曲率法。

最小曲率法是一種雙調和樣條函數插值法（Biharmonic Spline Interpolation）。所謂樣條函數，就是分段的光滑曲線函數，且各段的銜接處也光滑。如果說前面的趨勢法是用一張紙去擬合各點，最小曲率法就像是用有彈性的『橡膠膜』去擬合。並且，該方法和趨勢法最大的區別就是，該方法的『橡膠膜』必須要經過所有數據點，而趨勢法的紙只需要讓各點到紙面距離的均方根最小化，不要求通過所有點。

這個方法的具體要求是，每個點都照顧到，且要最小化『橡膠膜』表面的總曲率，找到對應的函數，然後得到空白區域的模擬地形。如何最小化總曲率？在每個點上，取『橡膠膜』對應函數的二階導數項，讓所有這些二階導數項的平方和最小，曲率就最小。這也是交給計算機去做的。

得到『橡膠膜』后，對『橡膠膜』畫等高線，就能得到最後的地圖。

最小曲率法畫的亮點是，畫出來的等值線平滑、漂亮，而且一般情況下可以很好地反應每一個數據點。但它的缺陷是，在坡度過於陡或者完全沒坡度的極端情況下，它的擬合不太靠譜。

這是 ArcGIS 官網的圖。這三幅圖用的同樣的數據，但從左到右分別是用反距離權重法、最小曲率法和克里金法畫的等高線，可以看出它們有明顯的區別：反距離權重法有很多『牛眼』，最小曲率法最平滑自然，克里金法對地形變化最敏感。

關於數據

俗話說，巧婦難為無米之炊。如果沒有靠譜的數據，再神的演算法，也只能畫出豬一樣的等高線地圖。一般來說，數據點分佈越均勻約密集，畫出來的圖就越好。如果數據分佈不均勻，有的地方數據密集，有的地方有大片的空白區域，那麼空白的地方的精準度肯定是遠遠不及數據密集的區域的。

最後吐個槽…

前面這些演算法都是理論上的，能讓大家知道如果把數據扔給計算機，比如 ArcGIS 等軟體，計算機是怎麼思考的，最後給你的等高線圖是怎麼得來的。

真實的操作中，ArcGIS 或者別的軟體比如 Petrel 畫出來的圖經常會很醜……而且太過於機械，不靈活。所以我們在做地質或地理的圖的時候，還是比較喜歡用手畫。但我們手畫的時候，也不完全按照 TIN 那麼來，畢竟那也太機械了。特別是地質的圖，在嚴格滿足數據點的前提下，對於空白處，經常要結合實際情況（比如當地的地形、沉積史、構造史、測井數據、地震數據等）進行自己的解讀和推測，從而畫出有地質意義的地圖。

比如畫某個地層的深度的等值線，用計算機畫，數學上說得過去，但有可能完全沒有地質意義。因為，地球上的事兒，啥奇怪的情況都可能發生。斷層（特別是逆斷層）、不整合等各種情況 ArcGIS 的默認演算法基本無能為力 。比如，如果我要畫下面圖裡的紅色的地層的上表面的等深度線：

這要讓 ArcGIS 來畫就不太好弄了。因此，那些演算法雖好，用的時候，還是得根據實際情況靈活調整的。比如左邊的逆斷層，就會有重疊區出現；右邊的正斷層，如果再多錯一點完全錯開，就會有空白區域出現。那如果工作量大，必須要用電腦呢？其實我們有 100 種方法調戲電腦，畫出靠譜的等高線，比如添加 pseudo-points 等等，只不過那就是另外一個話題了。

噗，寫完了才發現 Tag 話題裡面有個國中地理……好吧我說得太多了……

知乎專欄：地球的那些事兒（ShanYeTalking）

客官，這篇文章有意思嗎？

好玩！下載 App 接著看 (๑•ㅂ•) ✧

再逛逛吧 ˊ_>ˋ