李 響,王麗娜,張衛(wèi)東,楊 飛,楊振凱
1. 信息工程大學地理空間信息學院,河南 鄭州 450052; 2. 鄭州輕工業(yè)大學計算機與通信工程學院,河南 鄭州 450001; 3. 61512部隊,北京 100088
雙/多變量制圖(bivariate/multivariate mapping)是指能夠同時表示兩種或兩種以上地理現(xiàn)象的制圖方法[1]。目前的雙/多變量制圖方法主要分為兩類:
(1) 在專題地圖表示方法如等值區(qū)域法、分區(qū)統(tǒng)計圖表法等方法的基礎之上,運用多種視覺變量的組合或地圖符號的擴展的方法,從而達到表示兩種或者多種地理現(xiàn)象的目的。具體方法[2-4]包括使用不同色相的二維編碼(圖1(a))、Chemoff臉譜符號(圖1(b))、不同尺寸符號和明度組合(圖1(c))以及星形/雪花符號(圖1(d))等表達雙/多變量。
圖1 基于專題地圖表示方法的雙/多變量制圖方法示例Fig.1 Examples of traditional bivariate/multivariate cartographic methods
(2) 基于面域拓撲圖(area cartogram)表達雙/多變量。Cartogram是一種以地理對象面積或距離來表示其屬性特征的圖形表達方法,分為面cartogram和線cartogram。目前cartogram有拓撲地圖、統(tǒng)計地圖、示意地圖、屬性地圖或變形地圖等多種中文釋義,但這些中文釋義往往只強調(diào)了cartogram的某一方面[5],因此本文以下仍沿用英文cartogram這一術(shù)語。面cartogram主要包括4種類型[5-11]:簡單連續(xù)面cartogram(如矩形cartogram)、簡單離散面cartogram(如多林cartogram)、復雜離散面cartogram和復雜連續(xù)面cartogram。目前關(guān)于cartogram雙/多變量制圖主要是基于矩形cartogram和多林cartogram這類簡單圖形的面cartogram,在此基礎上通過符號擴展實現(xiàn)雙/多變量的表達。具體方法包括矩形cartogram和紋理、形狀視覺變量相結(jié)合[12](圖2(a));矩形cartogram與TreeMap相結(jié)合[13-15](圖2(b));多林cartogram與Chemoff臉譜符號相結(jié)合[16](圖2(c));多林cartogram基礎上通過對圓形符號的擴展[17](圖2(d))。
圖2 基于面cartogram表達雙/多變量制圖方法示例Fig.2 Examples of bivariate/multivariate cartographic methods based on area cartogram
雖然面cartogram中地理區(qū)域形狀扭曲變形,但其本身面積已經(jīng)表達了第1個變量,就為更多變量的表達以及不同變量空間分布規(guī)律的對比提供了更多可能。然而,目前基于面cartogram的雙/多變量制圖還存在以下兩個問題:
(1) 難以表達相鄰區(qū)域之間的基本狀況。目前基于面cartogram的雙/多變量制圖方法一般是以行政區(qū)域為單元進行變量統(tǒng)計,這樣能夠很好地反映各個獨立乃至整體區(qū)域的基本情況。但由于沒有考慮地理現(xiàn)象在區(qū)域中的位置,很難反映區(qū)域與區(qū)域之間,尤其是相鄰區(qū)域之間的基本情況。如圖3所示,無論是使用等值區(qū)域法或者多林cartogram的表示方法,區(qū)域單元線的略微偏移都會造成兩種迥然不同的結(jié)果。而且這兩種表示方法都無法反映兩個區(qū)域邊界地區(qū)學校密集、其他地區(qū)學校稀疏的地理分布特征。
圖3 學校分布案例Fig.3 The example of school distribution
(2) 不利于不同類型地理現(xiàn)象空間分布規(guī)律和差異特征的表達。由于矩形cartogram和多林cartogram難以保證整個區(qū)域的連續(xù)性,雖然圓形和矩形本身形狀簡單,面積易于估計和對比且易于分割,適合不同指標量化信息的傳遞,但難以有效表達不同地理現(xiàn)象的分布特征。
連續(xù)面cartogram通過變形來調(diào)整地理單元大小的同時保持整個區(qū)域的連續(xù)性,更能表現(xiàn)區(qū)域邊界處的地理要素分布狀況,同時也適用于表達不同地理現(xiàn)象的分布。因此,本文將連續(xù)面cartogram應用于雙/多變量的制圖表達中,首先對基于擴散模型的連續(xù)面cartogram生成算法進行優(yōu)化實現(xiàn),然后提出基于連續(xù)面cartogram的雙/多變量表達方法,最后通過兩個案例分別實現(xiàn)雙、多變量表達。研究方法的整體思路具體如圖4所示,大致可分為兩個步驟:
(1) 選擇基本變量并根據(jù)該變量值完成原地圖到連續(xù)面cartogram的轉(zhuǎn)換。在該轉(zhuǎn)換過程中,本文選擇最為經(jīng)典的基于擴散模型的連續(xù)面cartogram生成算法,并在具體實現(xiàn)中進行部分優(yōu)化改進,一定程度上提升了算法的效率。
(2) 雙/多變量在連續(xù)面cartogram中的表達。第2變量(如學校分布)通過空間內(nèi)插,重新計算其在連續(xù)面cartogram的位置,并與連續(xù)面cartogram進行疊加,從而實現(xiàn)面向雙變量的連續(xù)面cartogram可視化;其他變量(如學校能夠招收的學生數(shù)量及學生類型等),則通過擴展符號實現(xiàn)面向多變量的連續(xù)面cartogram可視化方法。
圖4 面向雙/多變量的連續(xù)面cartogram可視化方法的整體思路Fig.4 The general idea of visualization method of a continuous area cartogram for two or multiple variables
基于擴散模型的連續(xù)面cartogram生成算法來源于物理學中的擴散模型思想,模擬了液體從高密度向低密度擴散至密度均衡狀態(tài)的過程。在該轉(zhuǎn)換中,以人口密度地圖為例,允許人口從高密度地區(qū)流向低密度地區(qū),最終改變地理區(qū)域的形狀使得每個區(qū)域密度均衡,即首先通過格網(wǎng)密度計算每個格網(wǎng)邊界上的擴散速度,然后再將擴散速度進行積分從而求解格網(wǎng)邊界的位移,從而得到了變形后的格網(wǎng),再次求解格網(wǎng)密度,直至格網(wǎng)密度均衡,整個算法結(jié)束。由于該算法有效且快速,同時保留了區(qū)域間的連續(xù)性,在易讀性和地理準確性方面非常突出,此外在圖形扭曲度上保持了一定的靈活度,用戶可以在密度均衡程度和圖形扭曲程度之間作出調(diào)整,因此是應用最為廣泛的連續(xù)面cartogram生成算法[18-22]。該算法的基本思想和原理直觀明了,但在具體實現(xiàn)過程中,需要解決以下兩個基本問題:
(1) 面cartogram的邊界問題。在實際的面cartogram生成過程中,如果對面cartogram的邊界不作限制,人口數(shù)據(jù)會無限制地從高密度地區(qū)流入到低密度地區(qū),直至完全稀釋趨近于0;
(2) 任意時刻格網(wǎng)密度的求解問題。格網(wǎng)密度一直隨時間動態(tài)變化,且與位置、速度是相互作用的,因此難以求解任意時刻任意格網(wǎng)的密度。
文獻[18]中采用邊界條件[23]和傅里葉變換較好地解決了上述問題,并指出如果選定的邊界范圍足夠大于研究區(qū)域,那么不同類型的邊界條件對面cartogram變形的影響非常小。從算法實現(xiàn)和邏輯合理性上綜合考慮,采用了第2類邊界條件(諾依曼邊界條件[24])來限定擴散模型,即使用一個“墻”作為面cartogram的邊界,并認為在邊界“墻”處不存在數(shù)據(jù)的流動,速度為0。具體算法流程如下:
(1) 計算格網(wǎng)初始密度;
(2) 通過傅里葉變換(在諾依曼邊界條件下的傅里葉變換實際上是離散余弦變換)得到任意時刻的格網(wǎng)密度;
(3) 通過鄰近格網(wǎng)的密度差計算出該時刻格網(wǎng)邊界速度,由于限定了邊界條件,因此邊界處的格網(wǎng)邊界速度為0;
(4) 通過積分得出格網(wǎng)邊界位移;
(5) 判斷所有格網(wǎng)邊界位移是否發(fā)生變化,沒有變化則整個算法結(jié)束,否則回到步驟(2)。
擴散算法的基本原理及其實現(xiàn)如圖5所示。
圖5 基于擴散模型的等密度cartogram算法基本原理和實現(xiàn)Fig.5 Basic principles and implementation of diffusion-based cartogram
上述擴散算法在實現(xiàn)細節(jié)上仍存在以下兩個問題:
(1) 格網(wǎng)密度差異問題。當格網(wǎng)密度為0或者鄰接格網(wǎng)密度差異很大的時候,很可能會使得算法迭代多次都無法達到密度均衡的要求,從而導致算法失??;
(2) 積分循環(huán)次數(shù)過長問題。理論上只有時間趨近于無窮,才會生成等密度cartogram,同時為了確保積分的準確性,時間步長會設置的非常小,從而導致循環(huán)次數(shù)過長,運行效率低下。
本文提出通過格網(wǎng)密度補償和積分步長逐步試探的方法來解決上述兩個問題。
針對問題(1)本文將人口的數(shù)值做一個整體的正向偏移,具體偏移方法如式(1)所示
各格網(wǎng)正向偏移后的人口數(shù)=各格網(wǎng)人口數(shù)+
人口均值×偏移系數(shù)
(1)
偏移系數(shù)可以依據(jù)實際情況而定。這樣既可以避免零值出現(xiàn),也可以一定程度減小鄰接格網(wǎng)密度差異。
針對問題(2)本文提出了一種積分步長逐步試探法,能夠有效縮短循環(huán)次數(shù),提高運行效率。在計算機中進行積分運算,工程上通常采用最為廣泛的高精度單步算法-四階龍格-庫塔法[25],即已知時刻tn的位置rn,下一個時間間隔h的rn+1公式(2)如下
(2)
式中
(3)
式中,D為常量。為確保積分的準確性,本文限定了時間步長比率最大為4,如果依據(jù)式(3)計算出來的步長比率大于4,則步長按照4倍進行擴展,如果計算出來的步長比率小于4,則按照計算出來的步長比率進行擴展。算法流程如圖6所示,直至前后積分位移不再產(chǎn)生變化,算法結(jié)束。
圖6 積分步長逐步試探法流程Fig.6 Integral step size step-by-step trial method
以一個1 km×1 km的四格網(wǎng)為例,格網(wǎng)中的人口數(shù)如圖7所示,以左上角點為原點,建立二維平面直角坐標系,首先將人口數(shù)值依據(jù)式(1)做一個整體的正向偏移,該算例中將偏移系數(shù)設置為0.005,對格網(wǎng)進行了整體密度補償,并創(chuàng)建格網(wǎng)坐標,如圖7所示。
圖7 格網(wǎng)數(shù)值的總體正向偏移Fig.7 Overall forward offset of grid values
根據(jù)初始密度,可以計算出每個格網(wǎng)在水平方向和垂直方向的速度,如圖8所示。對于每個格網(wǎng)節(jié)點而言,數(shù)據(jù)從該格網(wǎng)節(jié)點所對應的左下方格網(wǎng)流入為正,流出為負。
圖8 依據(jù)初始密度計算水平和垂直方向的速度Fig.8 Calculate horizontal and vertical speeds based on initial density
表1 積分步長逐步試探法收斂過程
Tab.1 The convergence process of integral step size step-by-step method
序號時間/s步長位移誤差計算步長比率實際步長比率10 0.0026×10-12 80.37420.0020.0084.05×10-921.81430.010.0321.09×10-67.12440.0420.1284.783×10-53.343.3450.170.4281.371×10-42.712.7160.5981.1591.426×10-42.692.6971.7573.1161.469×10-42.672.6784.8748.3263.265×10-79.074913.20033.3041.611×10-15415.6841046.504133.2170——11179.721————
原始格網(wǎng)和變換后格網(wǎng)的結(jié)果對比如圖9所示。
雙/多變量在連續(xù)面cartogram中的表達方法主要有兩種:
(1) 面cartogram+點圖(dot map):可以將空間內(nèi)插后的第2變量的地理位置以點圖的形式疊加在面cartogram上,從而實現(xiàn)面向雙變量的面cartogram可視化。因此空間內(nèi)插是該表達方法的關(guān)鍵技術(shù),將在2.2節(jié)詳細討論。
圖9 原始格網(wǎng)和變換后格網(wǎng)對比圖Fig.9 The comparison of original grid and transformed grid
(2) 面cartogram+符號擴展:符號擴展的核心思想是靈活使用多種基本視覺變量,使得能夠表達兩種或者兩種以上的變量。通常在表示各種統(tǒng)計數(shù)據(jù)時,可以采用點狀符號的擴展,即點狀符號實際上是一個具有一定面積的二維平面,該平面的總面積表示事物的總量,用各分量占總量的百分比表示其內(nèi)部構(gòu)成。圓、環(huán)和扇面是最容易進行分割的圖形,也是人們應用最多的符號之一。將這些擴展的點狀符號疊加在面cartogram上,從而實現(xiàn)面向多變量的面cartogram可視化。
本文采用基于格網(wǎng)的雙線性插值方法來解決基于連續(xù)面cartogram的空間內(nèi)插,即在水平和垂直兩個方向上分別進行一次線性插值。如圖10所示,已知格網(wǎng)上的4個點分別在原地圖和面cartogram上的坐標。Q11、Q12、Q21和Q22表示4個點在地圖上的坐標,f(Q11)、f(Q12)、f(Q21)、f(Q22)表示4個點在面cartogram上的坐標。內(nèi)插格網(wǎng)中任意點P在面cartogram上的坐標f(P)。
首先在水平方向上進行一次線性內(nèi)插得到R1和R2兩個點,插值公式如下
R1(x,y1)
(4)
R2(x,y2)
(5)
圖10 二次線性內(nèi)插原理Fig.10 The principle of quadratic linear interpolation
然后在垂直方向上進行一次線性內(nèi)插,得到最終的結(jié)果P點,插值公式如下
(6)
雙線性插值與插值先后順序是無關(guān)的,無論先是從水平方向還是垂直方向上所得結(jié)果都相同。如果將4個已知點Q11、Q12、Q21和Q22作歸一化處理,將其分別簡化為坐標(0,0)、(0,1)、(1,0)和(1,1),則上述公式還可以進一步簡化為
f(x,y)=f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+
f(0,1)(1-x)y+f(1,1)xy
(7)
本文中基礎變量選用人口統(tǒng)計數(shù)據(jù)。由于國內(nèi)人口統(tǒng)計數(shù)據(jù)多按照省、市、自治區(qū)等行政區(qū)劃進行統(tǒng)計,人口按照格網(wǎng)進行統(tǒng)計的數(shù)據(jù)集比較少,精度也不太高,部分精度高的數(shù)據(jù)也較難以獲得。因此本文選用歐盟人口格網(wǎng)統(tǒng)計數(shù)據(jù)作為cartogram轉(zhuǎn)換的基礎數(shù)據(jù)。該數(shù)據(jù)集有兩個版本,最新版本數(shù)據(jù)為2016年,采用的是投影坐標系(ETRS89/LAEA),格網(wǎng)精度為1 km。
在試驗區(qū)域的選取上,考慮到便于實地查看和驗證,選取了兩個類型不同的德國南部城市,慕尼黑和奧格斯堡。慕尼黑是德國巴伐利亞州的首府,總面積達310 km2,2016年城市人口為146萬,是德國南部第一大城市,全德國第三大城市,僅次于柏林和漢堡。奧格斯堡位于德國南部巴伐利亞州西南,距離慕尼黑60多千米。人口截至2015年不足30萬人,面積146 km2。
通過編寫爬蟲程序從Google地圖上獲取了慕尼黑、奧格斯堡市的銀行、ATM、學校幼兒園等不同類型基礎設施的興趣點數(shù)據(jù)作為試驗中的第2變量。奧格斯堡的幼兒看護數(shù)據(jù)主要來自Google地圖以及奧格斯堡官方發(fā)布的奧格斯堡的幼兒看護(Kinderbetreuung in Augsburg,http:∥kinderbetreuung.augsburg.de/)本文對格斯堡市政府的全市幼兒園名冊進行了數(shù)字化,獲取了幼兒園的位置數(shù)據(jù)以及各幼兒園不同年齡段能夠接收的幼兒數(shù)量信息作為試驗研究的第3、第4變量。
3.2.1 可視化結(jié)果
按照統(tǒng)一的歐盟人口格網(wǎng)劃分,橫軸方向28個格網(wǎng),縱軸方向22個格網(wǎng)即可完全覆蓋慕尼黑城區(qū),每平方千米格網(wǎng)人口從0到2萬多不等,圖11分別是慕尼黑人口格網(wǎng)分布圖和根據(jù)人口數(shù)量生成的連續(xù)面cartogram。
圖11 慕尼黑人口格網(wǎng)分布圖和基于人口數(shù)量的連續(xù)面cartogram圖Fig.11 Population grid map of Munich and population grid cartogram of Munich
通過編寫爬蟲程序,從Google地圖上獲取了慕尼黑興趣點數(shù)據(jù)23 216條,以銀行(995條)和ATM機(454條)為例經(jīng)過空間內(nèi)插,分別疊加到基于人口數(shù)量的連續(xù)面cartogram圖上,其中紅的圓色點表示ATM機,紫的圓色點表示銀行,將它們分別疊加在地圖和連續(xù)面cartogram圖上,如圖12所示。
圖12 ATM和銀行分別疊加在連續(xù)面cartogram和人口格網(wǎng)圖Fig.12 The distribution of ATMs and banks in population grid cartogram and in population grid map
3.2.2 實例分析
通過兩幅圖對比并結(jié)合在Google Earth上觀察和實地走訪,得到如下分析結(jié)論:
(1) 在圖12(b)中能夠發(fā)現(xiàn),在人口稠密的中心城區(qū)ATM機和銀行分布較為密集,但是在連續(xù)面cartogram上還能夠進一步發(fā)現(xiàn),中心城區(qū)的1號區(qū)域比它周邊人口更為稠密區(qū)域的ATM機和銀行分布更多。通過在Google Earth上觀察以及實地走訪,1號區(qū)域正是從慕尼黑中心火車站通往瑪麗亞廣場的步行街地段,是慕尼黑重要的商業(yè)和旅游區(qū)域,如圖13所示,因此也就不難解釋為什么人口密度小于周邊區(qū)域,但是ATM機和銀行卻更多的原因。
圖13 Google Earth上觀察到的1號區(qū)域Fig.13 Area 1 on Google Earth
(2) 在圖12(b)中易于發(fā)現(xiàn),4號、5號和6號區(qū)域上存在較大的“洞”,即在這3個區(qū)域中幾乎沒有ATM機和銀行,雖然根據(jù)圖例能夠判斷這些區(qū)域人口較為稀少,但是無法確切獲知具體數(shù)量,因為對應的人口數(shù)量是一定范圍的(0—1179)。而在連續(xù)面cartogram上,這3個區(qū)域已經(jīng)急劇縮小,并未產(chǎn)生明顯的“洞”,因此并不引人注意。通過在Google Earth上觀察以及實地走訪,如圖14所示,這3個區(qū)域均存在非常廣袤的農(nóng)田,人口稀疏,因此ATM機和銀行數(shù)量較少。
圖14 Google Earth上觀察到的4號、5號和6號區(qū)域Fig.14 Area 4,5 and 6 on Google Earth
(3) 而在圖12(b)中難以發(fā)現(xiàn)到2號和3號區(qū)域上存在的“洞”,而在連續(xù)面cartogram上則非常明顯,這兩個區(qū)域人口都較為稠密,但是缺少ATM機和銀行的分布。
3.3.1 可視化結(jié)果
按照統(tǒng)一的歐盟人口格網(wǎng)劃分,橫軸方向17個格網(wǎng),縱軸方向25個格網(wǎng)即可完全覆蓋奧格斯堡城區(qū),如圖15所示。
在奧格斯堡幼兒看護按照年齡可以劃分為0~3歲(Krippe), 3~6歲(Kiga)學齡前兒童看護, 以及6歲以上(Hort)學齡兒童看護。圖16是幼兒看護地可接收不同年齡段的兒童數(shù)量在連續(xù)面cartogram和人口格網(wǎng)圖的分布情況。
圖15 奧格斯堡每平方千米人口分布圖Fig.15 Population grid map of Augsburg
3.3.2 實例分析
由于要表現(xiàn)幼兒看護地可接收的幼兒數(shù)量,需要使用到一定尺寸的分級圓符號,而不是單純的點符號。因此在人口格網(wǎng)圖中難以發(fā)現(xiàn)奧格斯堡的幼兒看護地有什么分布不均衡的地方,但是在連續(xù)面cartogram上可以看到在中心城區(qū)(1號區(qū)域)和中心城區(qū)下方(2號區(qū)域)出現(xiàn)了明顯的“洞”,如圖17所示。通過實地走訪和Google Earth可以驗證上述兩個區(qū)域確實沒有幼兒看護地,鄰近的也是在區(qū)域周邊。
圖16 幼兒看護地可接收不同年齡段兒童數(shù)量在連續(xù)面cartogram和人口格網(wǎng)圖的分布情況Fig.16 Child care locations which include the number of children of different ages in population grid cartogram and in population grid map
圖17 Google Earth上觀察到的1號和2號區(qū)域Fig.17 Area 1 and 2 on Google Earth
上述兩個案例中,由于均采用密度補償,因此試驗區(qū)域都沒有出現(xiàn)因密度差異較大而無法收斂的情況,且慕尼黑的邊界整體變形非常光滑。但由于奧格斯堡西南邊存在大量的森林(森林自然公園),人口極其稀少,即使進行了密度補償,仍然邊界出現(xiàn)一些急劇變形,這一點難以避免,如果對稀少人口進一步進行密度補償,則會影響cartogram本身的準確性。
此外,上述兩個案例由于都采用積分步長逐步試探法,因此都在迭代20次左右就達到收斂。從表2中可以看出,按照初始步長0.002 s開始循環(huán),整個慕尼黑地區(qū)的連續(xù)面cartogram生成需要經(jīng)歷20 712 s才能收斂,如果不采用積分步長逐步試探法,則需要循環(huán)10 356 000次才能收斂,采用積分步長逐步試探法循環(huán)迭代19次就達到收斂。整個奧格斯堡地區(qū)的連續(xù)面cartogram生成經(jīng)歷18 001.925 s,如果不采用積分步長逐步試探法,則需要循環(huán)9 000 963次才能收斂,采用積分步長逐步試探法循環(huán)迭代了24次就達到收斂。從表2中還能發(fā)現(xiàn)其收斂速度和地區(qū)大小并沒有相對關(guān)系,可能會同區(qū)域的人口分布存在一定關(guān)系,這一點還需要后續(xù)進一步驗證和研究。
表2 慕尼黑和奧格斯堡地區(qū)的連續(xù)面cartogram生成收斂過程
Tab.2 Convergence process of continuous area cartogram generation in Munich and Augsburg
序號慕尼黑區(qū)域(28×32格網(wǎng),單位 km)奧格斯堡區(qū)域(17×25格網(wǎng),單位 km)時間/s步長/s時間/s步長/s1 0 0.002 0 0.00220.0020.0080.0020.00830.010.0320.010.03240.0420.0850.0420.0250.1270.0910.0620.01460.2180.2800.0760.03170.4980.7440.1070.03781.2421.5270.1440.02492.7693.1370.1680.017105.9067.0360.1850.0691112.94216.7070.2540.2081229.64936.2530.4620.4671365.90277.1280.9291.08914143.030158.8852.0182.40715301.905384.0974.4255.48916686.012953.6199.91412.182
續(xù)表2
本文主要有兩點貢獻:①通過密度補償和積分步長逐步試探法改進和優(yōu)化基于擴散模型的等密度算法,有效縮短積分循環(huán)次數(shù),提高了算法效率;②將連續(xù)面cartogram應用到雙/多變量的制圖表達中,并通過試驗對比發(fā)現(xiàn),面向雙/多變量的連續(xù)面cartogram更能有效地發(fā)現(xiàn)不同變量空間分布差異明顯的區(qū)域,并且更容易察覺出變量之間差異的層次感,同時弱化差異不明顯的地方,有效引導人的視覺注意力,這為進一步分析內(nèi)部成因提供基礎,也能夠為城市規(guī)劃、政策制定等提供輔助參考決策。