引力模型最初源於物理學，用以描述物體間的相互作用力，後來被廣泛應用於社會科學領域，預測地理空間中的流動現象，例如貿易、通勤與人口遷徙。在網路科學中，此模型提供了一個強大的理論框架，用以解釋節點間的連結強度或流量大小，如何受到節點自身屬性（質量）與它們之間距離的共同影響。本章節的實證分析將此理論具象化，展示如何將抽象的「質量」與「距離」概念，轉化為航空網路中可量化的指標——機場的總旅客流量與地理距離。此過程不僅是數據處理的實踐，更是將一個經典的空間互動模型，應用於複雜網路結構，以揭示其內在組織原則與流動模式的關鍵步驟，特別是在參數估計的環節，我們將探索如何從觀測數據中反推出模型的基礎常數。
航空網路的引力模型實證分析與參數估計
本章節將深入探討如何將引力模型應用於實際的航空交通數據分析。我們將詳細闡述如何計算節點的「質量」（以機場總旅客流量為例），並利用 Haversine 公式計算機場間的距離，將這些資訊整合到網路結構中。更重要的是，我們將探討如何從數據中估計引力模型中的常數乘數，這是模型擬合的關鍵步驟，儘管我們將採用一種簡化的方法來處理。
距離與質量的整合
- Haversine 函數的實現與應用：
- 函數定義：
- 程式碼提供了
haversine(q, p)函數的完整實現，它接收兩個包含 (latitude, longitude) 的點q和p。 - 函數內部進行了度到弧度的轉換，並應用了 Haversine 公式來計算地球表面的大圓距離。
R_km = 6371：設定地球半徑為 6371 公里。d = R_km * c：最終計算出的距離d以公里為單位。
- 程式碼提供了
- 外部套件的考量：
- 文中也提及，有現成的 Python 套件（如
haversine和geopy）提供了更精確的距離計算方法，並且可能處理了地球非完美球體的因素。然而，為了完整性，這裡展示了手動實現。
- 文中也提及，有現成的 Python 套件（如
- 函數定義：
- 為邊添加距離屬性：
- 迭代邊：
for v, w in G_air.edges:：遍歷網路中的每一條邊。
- 獲取節點位置：
p_v = (G_air.nodes[v]['lat'], G_air.nodes[v]['long'])p_w = (G_air.nodes[w]['lat'], G_air.nodes[w]['long'])：從節點屬性中獲取兩個端點的緯度和經度。
- 計算並賦值距離：
G_air.edges[v, w]['distance'] = haversine(p_v, p_w)：調用haversine函數計算距離，並將結果作為'distance'屬性添加到邊(v, w)中。
- 迭代邊：
- 計算節點質量 (Mass)：
- 質量定義：
- 在航空交通的例子中，節點的「質量」被定義為該機場的總旅客流量。這代表了機場的規模和活躍程度。
- 文中也幽默地提到，根據應用場景，質量也可以是人口、甚至「平均松鼠大小」，強調了質量的抽象性。
- 計算方法：
- 總旅客流量就是機場的加權度數 (weighted degree)，其中權重是邊的旅客流量 (
count)。 degree = G_air.degree(weight='count')：使用 NetworkX 的degree方法，指定weight='count'來計算加權度數。nx.set_node_attributes(G_air, dict(degree), 'degree')：將計算出的加權度數（作為質量）設置為節點的'degree'屬性。
- 總旅客流量就是機場的加權度數 (weighted degree)，其中權重是邊的旅客流量 (
- 質量定義：
引力模型參數估計（簡化方法）
- 引力模型的假設：
- 流量 $I_{ij}$ 與節點質量乘積 $M_i M_j$ 成正比，與距離平方 $d_{ij}^2$ 成反比。
- $I_{ij} \propto \frac{M_i M_j}{d_{ij}^2}$
- 引入一個常數乘數 $k$：$I_{ij} = k \frac{M_i M_j}{d_{ij}^2}$。
- 常數乘數 $k$ 的估計：
- 挑戰：在實際應用中，準確估計這個常數 $k$ 是非常重要的，它代表了流量與質量、距離關係的整體尺度。
- 標準方法：通常會使用線性迴歸模型來擬合數據，將 $\log(I_{ij})$ 作為因變數，$\log(M_i) + \log(M_j) - 2\log(d_{ij})$ 作為自變數，來估計 $k$ 的對數。
- 簡化方法：
- 由於線性迴歸超出了本書的範圍，這裡採用一種更簡單的近似方法。
- 目標：對於每一條邊，計算一個「理論上」所需的常數 $k$，使得模型能夠「盡可能」地解釋觀察到的邊權重 (
count)。 - 程式碼邏輯：
g_list = []：創建一個列表來儲存計算出的常數。for v, w in G_air.edges():：迭代所有邊。if v >= w: continue：為了避免重複計算（例如，邊 (A, B) 和 (B, A)），只處理其中一條（例如，當節點名稱排序時，v < w）。try...except KeyError:：處理邊可能不存在或缺少關鍵屬性的情況。count = G_air.edges[v, w]['count']：獲取實際旅客流量。mass_v = G_air.nodes[v]['degree']：獲取節點 v 的質量（加權度數）。mass_w = G_air.nodes[w]['degree']：獲取節點 w 的質量。distance = G_air.edges[v, w]['distance']：獲取邊的距離。- 計算 $k$：
- 從公式 $I_{ij} = k \frac{M_i M_j}{d_{ij}^2}$，我們可以推導出 $k = I_{ij} \frac{d_{ij}^2}{M_i M_j}$。
k_val = count * (distance**2) / (mass_v * mass_w)：計算當前邊所需的常數 $k$。g_list.append(k_val)：將計算出的 $k$ 值添加到列表中。
- 對 $k$ 值的理解：
g_list中的值代表了每一條邊在獨立考慮時，為了使引力模型成立所需的 $k$ 值。- 這些值的平均值或中位數可以作為對整體常數 $k$ 的一個近似估計。
- 這些值的分佈也能提供關於模型擬合度的信息。如果 $k$ 的值非常分散，則說明引力模型對所有邊的解釋力不均勻。
程式碼實現細節
- Haversine 函數：
- 程式碼中提供了
haversine函數的實現，它接收兩個點的 (lat, long) 座標。
- 程式碼中提供了
- 添加距離屬性：
- 遍歷所有邊，獲取節點的地理屬性，調用
haversine計算距離，並將結果存入G_air.edges[v, w]['distance']。
- 遍歷所有邊，獲取節點的地理屬性，調用
- 計算節點質量：
- 使用
G_air.degree(weight='count')計算加權度數，並用nx.set_node_attributes設置為節點的'degree'屬性。
- 使用
- 計算常數乘數 $k$：
- 迭代邊，獲取
count,distance,mass_v,mass_w。 - 應用公式 $k = \text{count} \times (\text{distance}^2) / (\text{mass_v} \times \text{mass_w})$ 計算 $k$ 值。
- 將 $k$ 值存入
g_list。
- 迭代邊，獲取
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pathlib import Path
import math
# --- 假設 G_air_spatial 已載入並處理好節點屬性 (lat, long, count, degree)
---
# 為了程式碼的獨立性，這裡重新載入並處理部分數據
# 實際應用中，應接續前文的 G_air_spatial
# --- 重新創建模擬數據和網路 (如果前文未執行)
---
data_dir_sim = Path('./data_simulated')
carrier_data_path = data_dir_sim / 'BTS2018' / 'carrier.csv'
airport_db_path = data_dir_sim / 'partow' / 'GlobalAirportDatabase.txt'
# 模擬創建 carrier.csv
if not carrier_data_path.exists():
data_dir_sim.mkdir(parents=True, exist_ok=True)
mock_data_carrier = {
'ORIGIN_AIRPORT_ID': [10001, 10001, 10001, 10002, 10002, 10003, 10003, 10001, 10002, 10003, 10004, 10005, 10006, 10007],
'DEST_AIRPORT_ID': [10002, 10003, 10004, 10001, 10003, 10001, 10002, 10005, 10005, 10005, 10001, 10002, 10001, 10002],
'PASSENGERS': [1500, 800, 200, 1200, 900, 750, 600, 100, 150, 300, 50, 70, 1200, 1000],
'YEAR': [2018]*14, 'MONTH': [1]*14
}
df_mock_carrier = pd.DataFrame(mock_data_carrier)
df_mock_carrier.to_csv(carrier_data_path, index=False)
# 模擬創建 GlobalAirportDatabase.txt
if not airport_db_path.exists():
data_dir_partow = data_dir_sim / 'partow'
data_dir_partow.mkdir(parents=True, exist_ok=True)
mock_airport_db_content = """
1:JFK:John F. Kennedy International Airport:USA:40.639751:-73.778925
2:LAX:Los Angeles International Airport:USA:33.941556:-118.408530
3:ORD:Chicago O'Hare International Airport:USA:41.974209:-87.907321
4:ATL:Hartsfield-Jackson Atlanta International Airport:USA:33.640445:-84.427700
5:DEN:Denver International Airport:USA:39.856057:-104.673737
6:SFO:San Francisco International Airport:USA:37.619002:-122.371297
7:SEA:Seattle-Tacoma International Airport:USA:47.448981:-122.309310
8:MIA:Miami International Airport:USA:25.793199:-80.290594
9:DFW:Dallas/Fort Worth International Airport:USA:32.899811:-97.040416
10:LAS:McCarran International Airport:USA:36.084030:-115.153734
11:PHX:Phoenix Sky Harbor International Airport:USA:33.434250:-112.011590
12:IAH:George Bush Intercontinental Airport:USA:29.990163:-95.336770
13:BOS:Logan International Airport:USA:42.365599:-71.009667
14:MSP:Minneapolis–Saint Paul International Airport:USA:44.884794:-93.217711
15:PHL:Philadelphia International Airport:USA:39.874399:-75.242449
16:CLT:Charlotte Douglas International Airport:USA:35.214433:-80.947171
17:EWR:Newark Liberty International Airport:USA:40.689532:-74.174492
18:DTW:Detroit Metropolitan Airport:USA:42.211501:-83.353409
19:FLL:Fort Lauderdale–Hollywood International Airport:USA:26.074199:-80.150667
20:BWI:Baltimore/Washington International Airport:USA:39.177452:-76.668409
21:SLC:Salt Lake City International Airport:USA:40.789999:-111.979111
22:SAN:San Diego International Airport:USA:32.733819:-117.193751
23:TPA:Tampa International Airport:USA:27.977474:-82.531189
24:PDX:Portland International Airport:USA:45.589802:-122.597510
25:IAD:Washington Dulles International Airport:USA:38.953054:-77.456459
26:DCA:Ronald Reagan Washington National Airport:USA:38.851111:-77.037500
27:ANC:Ted Stevens Anchorage International Airport:USA:61.174355:-149.996333 # Alaska, will be removed
28:HNL:Daniel K. Inouye International Airport:USA:21.318694:-157.924694 # Hawaii, will be removed
"""
with open(airport_db_path, 'w') as f:
f.write(mock_airport_db_content.strip())
# --- 重新載入並處理數據
---
G_air_spatial = nx.Graph()
airport_locations = {}
try:
# 載入旅客流量
df_carrier = pd.read_csv(carrier_data_path)
df_carrier_2018 = df_carrier[df_carrier['YEAR'] == 2018]
for index, row in df_carrier_2018.iterrows():
origin = str(row['ORIGIN_AIRPORT_ID'])
dest = str(row['DEST_AIRPORT_ID'])
passengers = int(row['PASSENGERS'])
if origin == dest or passengers == 0: continue
if G_air_spatial.has_edge(origin, dest):
G_air_spatial.edges[origin, dest]['count'] += passengers
else:
G_air_spatial.add_edge(origin, dest, count=passengers)
if origin not in G_air_spatial: G_air_spatial.add_node(origin)
if dest not in G_air_spatial: G_air_spatial.add_node(dest)
# 載入地理位置
with open(airport_db_path) as f:
for row in f:
columns = row.strip().split(':')
if len(columns) >= 16:
code = columns[1]
country = columns[3]
if country == 'USA':
try:
lat = float(columns[14])
long = float(columns[15])
airport_locations[code] = (lat, long)
except ValueError: pass
# 添加地理屬性並篩選大陸機場
nodes_to_remove = []
for node_id in list(G_air_spatial.nodes()):
if node_id in airport_locations:
lat, long = airport_locations[node_id]
if 24 <= lat <= 49 and -125 <= long <= -66:
G_air_spatial.nodes[node_id]['lat'] = lat
G_air_spatial.nodes[node_id]['long'] = long
else:
nodes_to_remove.append(node_id)
else:
nodes_to_remove.append(node_id)
G_air_spatial.remove_nodes_from(nodes_to_remove)
print(f"網路準備就緒: {G_air_spatial.number_of_nodes()} 個節點, {G_air_spatial.number_of_edges()} 條邊。")
except Exception as e:
print(f"初始化網路時發生錯誤: {e}")
# --- Haversine 距離計算函數
---
def haversine(lat1, lon1, lat2, lon2):
"""Calculate the distance between two (lat, long) points using Haversine formula."""
R_km = 6371 # Earth radius in kilometers
lat1_rad, lon1_rad, lat2_rad, lon2_rad = map(math.radians, [lat1, lon1, lat2, lon2])
dlon = lon2_rad - lon1_rad
dlat = lat2_rad - lat1_rad
a = math.sin(dlat / 2)**2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2)**2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
distance = R_km * c
return distance
# --- 為邊添加距離屬性
---
for v, w in G_air_spatial.edges():
try:
p_v = (G_air_spatial.nodes[v]['lat'], G_air_spatial.nodes[v]['long'])
p_w = (G_air_spatial.nodes[w]['lat'], G_air_spatial.nodes[w]['long'])
G_air_spatial.edges[v, w]['distance'] = haversine(p_v[0], p_v[1], p_w[0], p_w[1])
except KeyError:
# 如果節點缺少 lat/long 屬性，則跳過此邊的距離計算
G_air_spatial.edges[v, w]['distance'] = None
print("已為邊添加距離屬性。")
# --- 計算節點質量 (加權度數)
---
# 使用 'count' 作為邊的權重來計算加權度數
degree_data = list(G_air_spatial.degree(weight='count'))
nx.set_node_attributes(G_air_spatial, dict(degree_data), 'mass') # 將質量命名為 'mass'
print("已計算節點質量 (總旅客流量)。")
# --- 簡化方法估計引力模型常數乘數 k
---
k_values = [] # 儲存計算出的 k 值
for v, w in G_air_spatial.edges():
# 為了避免重複計算 (例如，邊 (A,B) 和 (B,A))，我們只處理一次
# 這裡假設節點 ID 是字符串，可以進行比較
if v >= w:
continue
try:
count = G_air_spatial.edges[v, w]['count']
distance = G_air_spatial.edges[v, w]['distance']
mass_v = G_air_spatial.nodes[v]['mass']
mass_w = G_air_spatial.nodes[w]['mass']
# 確保所有值都有效，避免除以零或使用 None 值
if count is not None and distance is not None and distance > 0 and \
mass_v is not None and mass_w is not None and mass_v > 0 and mass_w > 0:
# 從公式 k = count * (distance^2) / (mass_v * mass_w) 計算 k
k_val = count * (distance**2) / (mass_v * mass_w)
k_values.append(k_val)
except KeyError:
# 如果邊缺少 'count', 'distance' 或節點缺少 'mass' 屬性，則跳過
continue
except Exception as e:
print(f"處理邊 ({v}, {w}) 時發生錯誤: {e}")
# --- 分析計算出的 k 值
---
if k_values:
k_values_np = np.array(k_values)
# 計算 k 值的統計數據
mean_k = np.mean(k_values_np)
median_k = np.median(k_values_np)
std_k = np.std(k_values_np)
print(f"\n估計的引力模型常數 k (簡化方法):")
print(f" - 平均值: {mean_k:.4f}")
print(f" - 中位數: {median_k:.4f}")
print(f" - 標準差: {std_k:.4f}")
# 繪製 k 值的直方圖，以了解其分佈
plt.figure(figsize=(10, 6))
plt.hist(k_values_np, bins=50, alpha=0.7, color='skyblue')
plt.axvline(mean_k, color='r', linestyle='dashed', linewidth=1, label=f'Mean k ({mean_k:.2f})')
plt.axvline(median_k, color='g', linestyle='dashed', linewidth=1, label=f'Median k ({median_k:.2f})')
plt.title("Distribution of Calculated Gravity Model Constant 'k'")
plt.xlabel("Value of 'k'")
plt.ylabel("Frequency")
plt.legend()
plt.grid(axis='y', alpha=0.5)
plt.show()
print("\nk 值分佈分析：")
print(" - 直方圖顯示了計算出的 k 值的分佈情況。")
print(" - 如果 k 值相對集中，表明引力模型對數據有較好的統一解釋力。")
print(" - 如果 k 值分佈廣泛，則意味著引力模型對不同邊的解釋力差異較大，可能需要更複雜的模型或考慮其他因素。")
else:
print("\n未能計算出任何引力模型常數 k 值。請檢查數據和網路屬性。")
@startuml
!define DISABLE_LINK
!define PLANTUML_FORMAT svg
!theme _none_
skinparam dpi auto
skinparam shadowing false
skinparam linetype ortho
skinparam roundcorner 5
skinparam defaultFontName "Microsoft JhengHei UI"
skinparam defaultFontSize 16
skinparam minClassWidth 100
start
:航空網路的引力模型實證分析與參數估計;:距離與質量的整合;
note right
Haversine 函數實現:
- 輸入: (lat, long) 座標對
- 計算: 大圓距離 (公里)
- 考量: 地球半徑, 度轉弧度
邊添加距離屬性:
- 遍歷邊 (v, w)
- 獲取節點 v, w 的 lat/long
- 調用 haversine 計算距離
- 存儲為 G_air.edges[v, w]['distance']
節點質量計算:
- 定義: 機場總旅客流量 (加權度數)
- 方法: G_air.degree(weight='count')
- 設置節點屬性: 'mass'
end note
:引力模型參數估計 (簡化方法);
note right
引力模型公式: I_ij = k * (M_i * M_j) / d_ij^2
- I: 實際流量 (count)
- M: 節點質量 (mass)
- d: 距離 (distance)
- k: 常數乘數 (待估計)
簡化估計 k:
- 對每條邊獨立計算 k 值
- k = count * (distance^2) / (mass_v * mass_w)
- 避免重複計算 (v < w)
- 處理缺失值和除以零
k 值列表 (k_values):
- 儲存所有計算出的 k 值
分析 k 值:
- 計算平均值、中位數、標準差
- 繪製 k 值分佈直方圖
- 評估模型對數據的統一解釋力
end note
stop
@enduml看圖說話：
此圖示總結了「航空網路的引力模型實證分析與參數估計」的內容，重點在於展示如何將計算出的距離和節點質量整合到網路中，並透過簡化的方法估計引力模型的核心參數 $k$。流程開頭首先聚焦於「距離與質量的整合」，說明了 Haversine 函數的應用、如何為邊添加距離屬性，以及如何計算節點質量（加權度數），接著詳細闡述了「引力模型參數估計 (簡化方法)」，解釋了引力模型公式，並描述了如何透過獨立計算每條邊所需的 $k$ 值，將這些值收集起來進行統計分析和視覺化，以評估模型的整體擬合情況。 縱觀引力模型在航空網路數據上的實證應用，我們不僅完成了從地理座標到網路距離、從旅客流量到節點質量的屬性整合，更重要的是，揭示了理論模型與真實世界數據之間的複雜互動。
透過對每條航線獨立估算常數 $k$ 的簡化方法，我們得到了一個分佈而非單一數值。此 $k$ 值分佈的離散程度，即成為評估模型解釋力的關鍵診斷指標。一個高度分散的結果直接表明，單純的「質量-距離」框架並不足以完整捕捉所有流量動態，這強烈暗示了票價結構、航空公司策略、樞紐機場效應等潛在變數的深遠影響。因此，這次參數估計的過程，其價值更接近於一次系統性的「壓力測試」，它所揭示的模型與現實的落差，為後續的特徵工程與模型優化提供了明確方向。
玄貓認為，對於數據分析師而言，此實證分析的核心價值，並非在於求得一個完美的常數以精確預測流量，而在於提供一個量化框架，讓我們得以系統性地識別並探索那些驅動複雜網路動態的、超越基礎引力法則的深層因素。
