一種基于sar數(shù)據(jù)的城區(qū)空氣動力學粗糙度反演方法
【技術領域】
[0001] 本發(fā)明涉及空氣動力學粗糙度計算技術領域，特別是涉及一種基于SAR數(shù)據(jù)的城 區(qū)空氣動力學粗糙度反演方法。
【背景技術】
[0002] 城市是人口密集、社會經(jīng)濟活動集中和資源相對聚集的空間地域?？焖僭黾拥慕?筑物、道路網(wǎng)等人工目標極大地改變了地表的粗糙度特性，下墊面粗糙特征的改變影響大 氣湍流和微氣象，改變湍流的傳輸、擴散和大氣污染物的沉降。因此，城市化進程導致了非 常復雜的城市邊界層結構，城市區(qū)域的地氣相互作用是當前研究的熱點之一?？諝鈩恿W 粗糙度（aerodynamic roughness) z。定義為地表上方風速為零的高度，是表征空氣動力學 特征的重要參數(shù)，其大小在一定程度上反映了地表與大氣的物質能量交換?？諝鈩恿W粗 糙度參數(shù)是城市氣象和大氣環(huán)境監(jiān)測研究的重要輸入?yún)?shù)，定量、準確獲取該參數(shù)對于正 確理解城市下墊面的空氣動力學特性具有重要意義。然而，城市是地表最粗糙的下墊面類 型，下墊面粗糙元高度在幾十米左右，且空間分布非常復雜，具有高度異質性，要在空氣動 力學粗糙度與粗糙元形態(tài)特征間建立準確的定量關系是很困難的。
[0003] 計算空氣動力學粗糙度的常規(guī)方法有氣象學方法和形態(tài)學方法。氣象學方法是根 據(jù)通量塔或者氣象站獲取的資料計算得到，但該方法通常需要風溫梯度數(shù)據(jù)或三維風速數(shù) 據(jù)，成本高，難度大，且只能獲取某個點上的空氣動力學粗糙度。形態(tài)學方法是利用基礎地 理信息數(shù)據(jù)和LiDAR數(shù)據(jù)獲取粗糙元的形態(tài)學參數(shù)，然后建立空氣動力學粗糙度與地面粗 糙元高度、粗糙元分布密度以及迎風面積等參數(shù)間的關系模型，該方法數(shù)據(jù)更新慢、成本高 的問題嚴重限制了其在城區(qū)下墊面粗糙特性研究中的應用。目前城區(qū)下墊面空氣動力學粗 糙度的研究仍處于初始階段，迫切需要針對城市下墊面發(fā)展適用的方法。
[0004] 合成孔徑雷達（Synthetic Aperture Radar, SAR)具有全天時全天候對地觀測能 力，對地物幾何結構和粗糙特性敏感，因此，為下墊面粗糙特性研究提供了重要手段?，F(xiàn)有 的利用SAR數(shù)據(jù)進行空氣動力學粗糙度反演的方法基本可以分為兩類。第一類方法是分類 查表法，即利用SAR圖像進行分類，利用查找表確定每種地物類型的空氣動力學粗糙度，進 而得出某個區(qū)域的空氣動力學粗糙度分布圖。這是目前最常用的方法，但這種方法通常用 于中尺度模型，當用于城市區(qū)域時要求非常精細的分類結果，這對于異質性很高的城市區(qū) 域來講是很困難的。第二類方法是在空氣動力學粗糙度與SAR圖像后向散射系數(shù)之間建立 關系模型，但現(xiàn)有研究未考慮到城市區(qū)域空氣動力學粗糙度的方向和尺度特性，忽略了城 市下墊面粗糙元空間分布的高度異質性，因此不能有效刻畫城區(qū)空氣動力學粗糙度的時空 分布特征。
[0005] 由此可見，上述現(xiàn)有的城區(qū)空氣動力學粗糙度計算方法顯然仍存在不足，如何能 創(chuàng)設一種充分考慮城市區(qū)域下墊面粗糙元分布具有高度異質性的特點，結合方向與尺度分 析技術的空氣動力學粗糙度反演方法，實屬當前重要研發(fā)課題之一。

【發(fā)明內(nèi)容】

[0006] 本發(fā)明要解決的技術問題是提供一種基于SAR數(shù)據(jù)的城區(qū)空氣動力學粗糙度反 演方法，使其充分考慮城市區(qū)域下墊面粗糙元分布的特點，借助方向與尺度分析技術，克服 現(xiàn)有的城區(qū)空氣動力學粗糙度計算方法的不足。
[0007] 為解決上述技術問題，本發(fā)明提供一種基于SAR數(shù)據(jù)的城區(qū)空氣動力學粗糙度反 演方法，其特征在于，包括如下步驟：
[0008] 第一步：對所述城區(qū)下墊面的SAR數(shù)據(jù)進行處理，提取出能表征所述下墊面粗糙 特性的極化散射參量；
[0009] 第二步：對利用氣象塔測量得到的多層風速風向梯度數(shù)據(jù)進行處理，計算出在不 同方向上的空氣動力學粗糙度；
[0010] 第三步：利用方向與尺度分析技術，通過分析所述極化散射參量與空氣動力學粗 糙度的相關性確定上風向扇形區(qū)域的最佳尺寸；
[0011] 第四步：從所述SAR圖像中提取所述最佳上風向扇形區(qū)域內(nèi)的極化散射參量，在 所述極化散射參量與城區(qū)空氣動力學粗糙度之間建立定量關系模型；
[0012] 第五步，利用第四步建立的定量關系模型反演得到城區(qū)空氣動力學粗糙度分布 圖。
[0013] 作為本發(fā)明的一種改進，所述第一步中對所述城區(qū)下墊面的SAR數(shù)據(jù)進行處理的 步驟包括：對SAR數(shù)據(jù)的輻射定標、濾波、幾何校正和配準處理，其中，全極化SAR數(shù)據(jù)處理 還包括極化分解和極化參數(shù)提取步驟，單/雙極化SAR數(shù)據(jù)處理還包括后向散射系數(shù)提取 步驟；以及分析城區(qū)各區(qū)域中地物的后向散射系數(shù)、極化方位角、同極化與交叉極化相位 差、面散射、二次散射、體散射極化散射參量的差異性；最終確定最能表征城區(qū)下墊面粗糙 度特性的極化散射參量。
[0014] 進一步改進，所述第二步中對多層風速風向梯度數(shù)據(jù)進行處理的步驟包括：
[0015] 首先剔除存在明顯錯誤的數(shù)據(jù)，然后以理查遜數(shù)確定大氣穩(wěn)定度，再選用滿足中 性層結條件的風速梯度數(shù)據(jù)，并根據(jù)式（1):
[0017] 其中，Uti為近地面層摩擦速度，k為Von Karman常數(shù)，z為離地面的實際高度，z。 為空氣動力學粗糙度，d。為零平面位移，u w為離地面實際高度z時的風速；
[0018] 利用擬合迭代法計算得到不同方向上的空氣動力學粗糙度.，其中i表示空氣 動力學粗糙度計算值的編號，Θ i表示第i個空氣動力學粗糙度對應的風向。
[0019] 進一步改進，根據(jù)式（1)和擬合迭代法計算z。和d。的具體步驟為：選取不同d。值 山，由ιι ω和ln(z-d J做線性回歸，直到％和ln(z-d J相關系數(shù)達到最大時山為d。的真 值，此時由線性斜率i^/k可得Uti，由截距UtlZk · ln(z。）可得到z。的真值。
[0020] 進一步改進，所述第三步中確定上風向扇形區(qū)域最佳尺寸的步驟包括：
[0021] (1)根據(jù)第二步獲取的不同方向上空氣動力學粗糙度；確定每個地面觀測結 果對應的時間和風向θ ι;
[0022] (2)根據(jù)空氣動力學粗糙度地面觀測結果對應的時間，選擇相應時相的SAR圖 像數(shù)據(jù)；
[0023] (3)定義上風向扇形區(qū)域半徑和夾角的取值范圍和步長；
[0024] (4)利用氣象觀測資料計算得到的空氣動力學粗糙度數(shù)據(jù)的時間及風向設定不同 大小的扇形區(qū)域，提取每個空氣動力學粗糙度數(shù)據(jù)對應的不同大小扇形區(qū)域內(nèi)的極化 音夂射參量;
[0025] (5)對所述每個空氣動力學粗糙度數(shù)據(jù)和與其對應的不同大小扇形區(qū)域內(nèi)的 極化散射參量μ進行相關性分析，計算兩者之間的相關系數(shù)，進而確定相關系數(shù)最大值 對應的尺寸為最佳上風向扇形區(qū)域半徑和夾角。
[0026] 進一步改進，所述步驟（3)中上風向扇形區(qū)域的半徑取值范圍為1000~7000m，步 長為500m，扇形區(qū)域的夾角取值為30度、45度或60度。
[0027] 進一步改進，所述極化散射參量V與城區(qū)空氣動力學粗糙度z。之間的指數(shù)關系模 型如下式：
[0029] 其中A和B表示通過擬合得到的模型參數(shù)。
[0030] 進一步改進，所述第五步中還包括利用多時相SAR數(shù)據(jù)分析所述城區(qū)空氣動力學 粗糙度的時空分布規(guī)律。
[0031] 采用上述的技術方案，本發(fā)明至少具有以下優(yōu)點：
[0032] 本發(fā)明提出了建立SAR數(shù)據(jù)與城區(qū)空氣動力學粗糙度定量關系模型的方法，可以 獲取大范圍的城區(qū)空氣動力學粗糙度的時空分布，為相關的氣象和氣候模型提供更精確的 輸入?yún)?shù)。
[0033] 本發(fā)明通過方向和尺度特性分析，可以根據(jù)下墊面特點確定空氣動力學粗糙度反 演計算所需的上風向扇形區(qū)域的最佳尺寸，能更精細地刻畫城市下墊面粗糙度特性。
【附圖說明】
[0034] 上述僅是本發(fā)明技術方案的概述，為了能夠更清楚了解本發(fā)明的技術手段，以下 結合附圖與【具體實施方式】對本發(fā)明作進一步的詳細說明。
[0035] 圖1是本發(fā)明基于SAR數(shù)據(jù)的城區(qū)空氣動力學粗糙度反演方法的流程圖。
【具體實施方式】
[0036] 針對現(xiàn)有的城區(qū)空氣動力學粗糙度反演方法存在的問題，本發(fā)明充分考慮城市下 墊面粗糙元異質性高的特點，借助方向和尺度分析技術，提出了一種基于SAR數(shù)據(jù)的城區(qū) 空氣動力學粗糙度反演方法。參照附圖1所示，該方法具體包括如下步驟：
[0037] 第一步：首先對SAR數(shù)據(jù)進行處理，包括SAR數(shù)據(jù)的輻射定標、濾波、幾何校正和配 準處理等，其中，全極化SAR數(shù)據(jù)處理還包括極化分解和極化參數(shù)提取步驟，單/雙極化SAR 數(shù)據(jù)處理還包括后向散射系數(shù)提取步驟；然后針對城市區(qū)域的特點，選取最具代表性的地 類，如住宅區(qū)、高層建筑區(qū)、中央商務區(qū)、水體、綠地、道路等類型，分析這些地物的后向散射 系數(shù)、極化方位角、同極化與交叉極化相位差、面散射、二次散射、體散射等極化散射參量的 差異性，最終確定最能反映下墊面粗糙度特性的極化散射參量。
[0038] 第二步：對利用氣象塔測量得到的多層風速風向梯度數(shù)據(jù)進行處理，首先剔除存 在明顯錯誤的數(shù)據(jù)，然后以理查遜數(shù)確定大氣穩(wěn)定度，選用滿足中性層結條件的風速梯度 數(shù)據(jù)，并根據(jù)莫寧一奧布霍夫相似理論，中性層結條件下近地面層水平風速分布滿足對數(shù) 律，即：
[0040] 其中Uti為近地面層摩擦速度，k為Von Karman常數(shù)，一般取為0. 4, z為離地面的 實際高度，為空氣動力學粗糙度，d。為零平面位移，u e)為離地面實際高度z時的風速。
[0041] 利用式⑴通過擬合迭代法計算z。和d。的具體步驟為：選取不同d。值d i，由?!ω 和In (Z-Cl1)做線性回歸，直到1^}和In (z-d J相關系數(shù)達到最大時山為d。的真值，由線性 斜率取i^/k可得Uti，由截距UtlZk · ln(z。）可得到z。的真值。