国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

蝸輪齒面刀痕打磨接觸力獲取方法研究

2023-10-07 12:58:22常旭趙康陳家榮
關(guān)鍵詞:蝸輪齒面卡爾曼濾波

常旭,趙康,陳家榮

(上海工程技術(shù)大學(xué)機(jī)械與汽車工程學(xué)院,上海 松江 201600)

目前,包絡(luò)環(huán)面桿副中的蝸輪主要通過滾齒機(jī)進(jìn)行加工,加工之后的蝸輪齒面不可避免地會產(chǎn)生刀痕,需要進(jìn)一步在不同負(fù)載(25%、50%、75%、滿負(fù)載)下進(jìn)行多次跑合[1-3],以消除蝸輪齒面刀痕并達(dá)到齒面嚙合精度,這種去除蝸輪齒面刀痕的方法工藝復(fù)雜,效率很低.為提高蝸輪加工效率,在蝸輪經(jīng)過滾切之后,使用工業(yè)機(jī)器人打磨蝸輪齒面刀痕,完成對蝸輪齒面的粗加工,并在滿負(fù)載狀態(tài)下進(jìn)行一次跑合即可完成蝸輪精加工.工業(yè)機(jī)器人打磨加工雖然具有效率高、成本低、一致性高等優(yōu)點(diǎn),但是尚需精確獲取工業(yè)機(jī)器人末端與蝸輪齒面之間的接觸力以實(shí)現(xiàn)柔順控制.

在接觸力精確獲取的研究上,文獻(xiàn)[4-5]建立的重力補(bǔ)償計算模型只考慮了機(jī)器人末端工具重力,沒有考慮機(jī)器人安裝傾角、傳感器零點(diǎn)和噪聲影響,而且僅進(jìn)行了仿真分析,沒有實(shí)驗(yàn)驗(yàn)證;文獻(xiàn)[6-7]提出了一種基于力傳感器重力補(bǔ)償?shù)臋C(jī)器人柔順控制方法,建立了重力補(bǔ)償數(shù)學(xué)模型,通過采集不同位姿下力傳感器數(shù)據(jù),計算機(jī)器人底座安裝傾角、力傳感器零點(diǎn)數(shù)據(jù)、末端工具重力及重心坐標(biāo),但是沒有考慮力傳感器采集的數(shù)據(jù)含有噪聲的影響;文獻(xiàn)[8-9]針對六維力傳感器的重力補(bǔ)償提出了一套完整的標(biāo)定與計算方法,綜合考慮了機(jī)器人安裝傾角、傳感器零點(diǎn)、工具重力和重心以及噪聲的影響,但是計算過程較為復(fù)雜;文獻(xiàn)[9]根據(jù)力等效原理,建立了零點(diǎn)校正和重力補(bǔ)償數(shù)學(xué)模型,通過至少3 個機(jī)器人的不同位姿計算出傳感器零點(diǎn)、安裝傾角以及重力補(bǔ)償參數(shù),但是缺少力跟蹤控制實(shí)驗(yàn),沒有驗(yàn)證實(shí)際應(yīng)用中力跟蹤效果.

本文主要對工業(yè)機(jī)器人與蝸輪齒面交互過程中真實(shí)力感知信息的獲取方法進(jìn)行了研究, 為實(shí)現(xiàn)工業(yè)機(jī)器人柔順打磨蝸輪齒面刀痕奠定基礎(chǔ).

1 接觸模型構(gòu)建

在工業(yè)機(jī)器人打磨蝸輪齒面刀痕過程中,末端安裝的電主軸始終與蝸輪齒面平行,末端與工件的坐標(biāo)系如圖1 所示,工業(yè)機(jī)器人傳感器坐標(biāo)系{Ts}姿態(tài)在打磨加工時始終與打磨臺坐標(biāo)系{Tb2}保持一致.工業(yè)機(jī)器人打磨蝸輪齒面刀痕過程較為復(fù)雜,直接對打磨進(jìn)行理論分析難度系數(shù)較大,并且加工傾角、打磨深度、進(jìn)給速度等外在因素的細(xì)微變化都會影響打磨過程中打磨力的變化.在工業(yè)機(jī)器人打磨加工過程中,為保護(hù)蝸輪齒面不被過切,在設(shè)置軌跡時要留有余量,設(shè)置的末端軸向打磨量遠(yuǎn)大于徑向打磨量,在這種情況下需要考慮蝸輪齒面受到的徑向力,即法向力.假設(shè)在打磨時所產(chǎn)生的力平行于高速電機(jī)主軸,那么在理論分析時就可以將打磨過程轉(zhuǎn)變?yōu)榇蚰デ娴拇蚰ミ^程,在打磨時需要控制垂直于打磨軌跡的法向力,即X 軸方向的力.

圖1 工業(yè)機(jī)器人打磨位姿分析Fig.1 Analysis of grinding pose of industrial robot

在蝸輪齒面刀痕打磨系統(tǒng)的笛卡爾操作空間中,對打磨工具坐標(biāo)系{Tt}進(jìn)行力/位分解,Xt方向進(jìn)行恒力控制,且始終保持水平,即打磨姿態(tài)為b1Rt=[1 0 0]T,Yt、Zt方向進(jìn)行位置控制.由運(yùn)動學(xué)模型計算得到打磨工具末端的位置為

式(1)中:Tb1表示工業(yè)機(jī)器人基坐標(biāo)系,Tf表示工業(yè)機(jī)器人末端坐標(biāo)系,Tt表示打磨工具末端坐標(biāo)系,Tb2表示打磨臺中心坐標(biāo)系,Ts表示六維力傳感器坐標(biāo)系,Tj表示工件坐標(biāo)系.

2 工業(yè)機(jī)器人末端與蝸輪齒面交互力預(yù)處理

2.1 力噪聲預(yù)處理

卡爾曼濾波器[10]是一種由線性系統(tǒng)狀態(tài)方程通過迭代更新最優(yōu)估計系統(tǒng)最佳狀態(tài)的濾波器,即使有諸多不確定信息,卡爾曼濾波器也能找到不易察覺的相關(guān)性.卡爾曼濾波器在對力感知信息進(jìn)行估算的過程中,只需要考慮系統(tǒng)的過程噪聲、測量噪聲以及系統(tǒng)當(dāng)前狀態(tài),不需要對整個過程的數(shù)據(jù)進(jìn)行采集、保存,占用內(nèi)存少,因此非常適用于力感知信息實(shí)時在線處理[11].

將卡爾曼濾波器的狀態(tài)變量設(shè)計為六維,分別是六維力傳感器采集的三維力信息:Fx、Fy、Fz與三維力變化信息:,觀測變量為三維力信息:Fx、Fy、Fz.六維力傳感器采集力感知信息的迭代方程可簡化為如下二階方程

式(2)可化簡為

系統(tǒng)的狀態(tài)方程為

式(3)中:Xt表示工業(yè)機(jī)器人打磨系統(tǒng)t 時刻狀態(tài)變量,A 表示工業(yè)機(jī)器人打磨系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣,T表示采樣周期,Ut表示工業(yè)機(jī)器人打磨系統(tǒng)過程噪聲.

工業(yè)機(jī)器人打磨系統(tǒng)的觀測方程可表示為

式(4)中:Zt+1 表示工業(yè)機(jī)器人打磨系統(tǒng)t+1 時刻的觀測變量,H 表示工業(yè)機(jī)器人打磨系統(tǒng)觀測矩陣,Vt表示工業(yè)機(jī)器人打磨系統(tǒng)的測量噪聲.卡爾曼濾波器以迭代的方式對系統(tǒng)進(jìn)行最優(yōu)在線估計,更新的過程包括時間更新和狀態(tài)更新,其中狀態(tài)更新包括模型預(yù)測和觀測兩個過程,迭代方程由式(5)表示.

2.2 重力補(bǔ)償

2.2.1 工業(yè)機(jī)器人安裝傾角和末端工具重力計算一般情況下,工業(yè)機(jī)器人在安裝時很難保證工業(yè)機(jī)器人處于水平平面,會與水平面存在一定的傾斜角度,因此在實(shí)際中工業(yè)機(jī)器人的基座標(biāo)系與實(shí)際坐標(biāo)系并不重合.如圖2 所示,建立工業(yè)機(jī)器人的坐標(biāo)系,定義Tw-X0Y0Z0為世界坐標(biāo)系,Z0軸豎直向上,T1-X1Y1Z1為機(jī)器人基座標(biāo)系,Ts-XsYsZs為六維力傳感器坐標(biāo)系,工業(yè)機(jī)器人基座標(biāo)系由世界坐標(biāo)系先繞X0軸旋轉(zhuǎn)u,再繞Y0軸旋轉(zhuǎn)v 得到.六維力傳感器實(shí)際測量數(shù)據(jù)記為:F=[f x,f y,fz]、M=[mx,my,mz],六維力傳感器的零點(diǎn)數(shù)據(jù)記為,打磨工具在力六維傳感上的作用力與力矩記為:G=[g,g,g]、M=[m,m,m].

工業(yè)機(jī)器人未與環(huán)境發(fā)生交互作業(yè)時,六維力傳感器實(shí)際上測量的數(shù)據(jù)是打磨工具在六維力傳感器上產(chǎn)生的力與力矩和六維力傳感器零點(diǎn)之和,所以

式(6)中:u、v 表示工業(yè)機(jī)器人的安裝傾角.

式(7)中:l 是3×3 單位矩陣.

測量N(N≥3)組工業(yè)機(jī)器人不同姿態(tài)下六維力傳感器數(shù)據(jù),代入式(7)得到

進(jìn)一步可化簡成

通過矩陣變換得到

由式(8)可計算出六維力傳感器的零點(diǎn)分力以及常數(shù)l、l、l.

同時可計算出工業(yè)機(jī)器人底座的傾斜角以及末端打磨工具重力為

2.2.2 六維力傳感器零點(diǎn)與末端工具中心坐標(biāo)計算末端打磨工具重心在六維力傳感器坐標(biāo)系Ts中的位置記為,在六維力傳感器上產(chǎn)生的力和力矩如圖3 所示,可求得打磨工具重力在六維力傳感器上的力矩計算公式

式中:P為塊段中鋯英石、或鈦鐵礦礦物量,t;V為塊段體積,m3);C為塊段鋯英石或鈦鐵礦的平均品位,kg·m3。

圖3 末端工具作用在六維力傳感器上的力矩示意圖Fig.3 Torque diagram of the end tool acting on the six-dimensional force sensor

聯(lián)立式(6)與式(9)可得

故式(10)可整理成

同樣選取N(N≥3)個工業(yè)機(jī)器人末端打磨工具不同姿態(tài)下的六維力傳感器數(shù)據(jù),代入式(11)可以求得工業(yè)機(jī)器人末端打磨工具重力在六維力傳感器上的坐標(biāo)p、p、p以及b、b、b.

由聯(lián)合式(10)和式(11)可求得

綜上公式計算,可以計算出六維力傳感器零點(diǎn)值:,工業(yè)機(jī)器人末端打磨工具重力G 和重心坐標(biāo)q=[q,q,q].

2.2.3 重力補(bǔ)償值以及負(fù)載重力計算由式(8)(11)可以求出打磨工具重力在六維力傳感器坐標(biāo)系上的力和力矩分量.故需要在六維力傳感器坐標(biāo)系上補(bǔ)償?shù)牧土貫?/p>

對六維力傳感器進(jìn)行重力補(bǔ)償之后,打磨工具末端接觸的真實(shí)力和力矩為

3 基于位置的阻抗控制

基于位置的阻抗控制[12]是一種適用于解決工業(yè)機(jī)器人與環(huán)境交互問題的控制策略,它建立了工業(yè)機(jī)器人與環(huán)境之間接觸力誤差和工業(yè)機(jī)器人運(yùn)動變量之間的函數(shù)關(guān)系,實(shí)現(xiàn)了工業(yè)機(jī)器人末端力控制.因?yàn)槠洳恍枰_的動力學(xué)模型,所以被廣泛地應(yīng)用于工業(yè)機(jī)器人控制領(lǐng)域[13,14].基于位置的阻抗控制架構(gòu)如圖4 所示,主要由位置控制內(nèi)環(huán)和力控制外環(huán)組成,其中Xd表示期望位置,Xe表示實(shí)際位置,Fe表示六維力傳感器采集的實(shí)際接觸力,即工業(yè)機(jī)器人與環(huán)境的接觸力,Fe'表示預(yù)處理后的實(shí)際接觸力,ΔF表示力誤差,Fd表示期望接觸力,ΔX 表示期望位置與實(shí)際位置的偏差,即對位置的阻抗補(bǔ)償量,力誤差ΔF 與位置補(bǔ)償量ΔX 的阻抗關(guān)系可用式(13)表示.

圖4 基于位置的阻抗控制架構(gòu)圖Fig.4 Architecture diagram of impedance control based on position

式(12)中:M 表示質(zhì)量參數(shù)矩陣,B 表示阻抗參數(shù)矩陣,K 表示剛度參數(shù)矩陣.

由于位置的阻抗控制中力控制和位置控制相互耦合以及六維力傳感器采集的力感知信息以一定的頻率傳送給連續(xù)的阻抗控制,因此為了使力控制與位置控制盡可能地分離,需要對連續(xù)的阻抗控制進(jìn)行時域離散化[15-17].另外,為了便于對工業(yè)機(jī)器人末端三個方向進(jìn)行柔順力控制以及方便調(diào)節(jié)阻抗參數(shù),本文簡化阻抗控制模型,即認(rèn)為機(jī)器人末端在X/Y/Z 三個方向的運(yùn)動相互解耦[18-19].以其中的X 方向?yàn)槔?由圖3 可得

式(13)中:e 表示期望位置與實(shí)際位置之差,ke表示環(huán)境剛度,表示機(jī)器人末端與環(huán)境的接觸力,xc表示未發(fā)生變形前的位置,xe表示環(huán)境當(dāng)前實(shí)際位置.

由式(13)可得:

將式(13)(14)進(jìn)行La-place 變換并化簡可得

所以系統(tǒng)穩(wěn)態(tài)力跟蹤誤差為

顯然要使力跟蹤穩(wěn)態(tài)誤差為零,滿足式(15)中任一條件即可,但在實(shí)際的工業(yè)機(jī)器人交互作業(yè)過程中,由于工業(yè)機(jī)器人末端剛度跟隨環(huán)境實(shí)時變化,因此實(shí)時、準(zhǔn)確地獲取剛度值很困難.為了簡化處理,本文將剛度參數(shù)kd設(shè)置為零,因此式(13)的阻抗控制函數(shù)可以化簡成

令x=v=,其中v 代表速度,則式(13)可化簡為

對式(16)求解可得

將式(17)離散化,假設(shè)系統(tǒng)采樣周期為Δt,在采樣周期內(nèi)慣性參數(shù)md、阻尼參數(shù)bd為常數(shù),令,可以得到工業(yè)機(jī)器人末端接觸力與機(jī)器人控制速度的表達(dá)式為

故位置補(bǔ)償量ΔX 的表達(dá)式為

4 實(shí)驗(yàn)驗(yàn)證與仿真分析

4.1 卡爾曼濾波器在線濾波實(shí)驗(yàn)

根據(jù)實(shí)驗(yàn)現(xiàn)場實(shí)際情況,卡爾曼增益初始值=1,為了保證卡爾曼濾波器收斂,先驗(yàn)估計協(xié)方差設(shè)定為:.將系統(tǒng)過程噪聲和系統(tǒng)采集觀測噪聲分別設(shè)定為

實(shí)驗(yàn)時,六維力傳感器實(shí)時采集打磨過程中的力感知信息,然后使用卡爾曼濾波器實(shí)時處理采集得到力感知信息.本文以X 方向的力信息濾波處理結(jié)果為例說明卡爾曼濾波器的處理效果.圖5、6 表示對平穩(wěn)力感知信息的處理效果,在未處理前力信息在±3 N 之間波動,卡爾曼濾波器處理后,波動值降到±1 N之間,使采集得到的力信息更加平穩(wěn),降低了外界因素對六維力傳感器采集力感知信息的干擾.圖3-7表示非平穩(wěn)力感知信息的處理效果,結(jié)果表明卡爾曼濾波器能夠適應(yīng)采樣值的變化,達(dá)到減少外界影響的目的.

4.2 重力補(bǔ)償實(shí)驗(yàn)

為驗(yàn)證本文重力補(bǔ)償算法的效果,根據(jù)傳感器標(biāo)定及重力補(bǔ)償算法進(jìn)行工業(yè)機(jī)器人在線重力補(bǔ)償與傳感器零點(diǎn)標(biāo)定實(shí)驗(yàn).首先測量工業(yè)機(jī)器人安裝傾角、末端工具重力、六維力傳感器零點(diǎn)以及末端工具中心.將工業(yè)機(jī)器人末端的位姿設(shè)置為如圖7 所示的6 種姿態(tài).為了減少未知因素帶來的誤差,每種姿態(tài)連續(xù)采集600 組數(shù)據(jù)求平均值, 并且每組數(shù)據(jù)采集時均經(jīng)過卡爾曼濾波器濾波,6 種姿態(tài)下六維力傳感器的力感知信息測量值如表1 所示.

圖7 機(jī)器人姿態(tài)示意Fig.7 Robot posture diagram

my/(N.m)0.982 3 0.000 0-1.013 6 10.012 5-0.114 5-9.821 4姿態(tài)123456重力方向-z-y+z+x+y-x fx/N-0.824 3-0.100 0 0.134 5 99.016 2 0.262 4-98.552 0 fy/N 0.734 5-97.175 9-0.127 5 0.000 0 98.163 2-0.624 3 fz/N 98.753 5 0.000 0 98.921 5 0.245 8-0.413 5-0.914 6 mx/(N.m)-1.154 3 9.813 2 1.234 6 0.000 0-9.811 4-0.213 5 mz/(N.m)0.000 0-0.981 2 0.137 8-0.963 9 0.983 8 1.247 9

根據(jù)表中數(shù)據(jù)以及本文所提及的重力補(bǔ)償算法,可以計算得出六維力傳感器的零點(diǎn)、安裝傾角、工具重力和打磨工具在傳感器坐標(biāo)系下質(zhì)心的位置

在工業(yè)機(jī)器人打磨蝸輪齒面刀痕的過程中,末端位姿會隨著打磨軌跡不斷變化,末端負(fù)載在六維力傳感器中的分量也會隨著末端位姿不斷改變,所以在標(biāo)定力傳感器零點(diǎn)和工具重力之后,也要驗(yàn)證標(biāo)定算法和重力補(bǔ)償算法在工業(yè)機(jī)器人運(yùn)動過程中的補(bǔ)償效果.具體實(shí)驗(yàn)方法是:使用示教器控制工業(yè)機(jī)器人夾持打磨工具在自由空間運(yùn)動任意一條軌跡,工業(yè)機(jī)器人運(yùn)行速度設(shè)置為7 mm/s,采集運(yùn)動過程中六維力傳感器讀數(shù),且采集時力感知信息已經(jīng)經(jīng)過卡爾曼濾波器處理.圖8 表示未經(jīng)過重力補(bǔ)償時六維力傳感器的讀數(shù),圖9 表示重力補(bǔ)償后六維力傳感器的讀數(shù).由圖8、9 中數(shù)據(jù)可以看出六維力傳感器未進(jìn)行補(bǔ)償時,傳感器讀數(shù)偏移較大且隨著工業(yè)機(jī)器人位姿變化而變化,經(jīng)過重力補(bǔ)償后,六維力傳感器數(shù)據(jù)的讀數(shù)誤差維持在±0.4N 之間,力矩誤差維持在±0.16 N·m 之間,其讀數(shù)均在零點(diǎn)附近波動,驗(yàn)證了零點(diǎn)標(biāo)定、重力補(bǔ)償算法的有效性.

4.3 機(jī)器人柔順打磨實(shí)驗(yàn)

工業(yè)機(jī)器人柔順打磨實(shí)驗(yàn),是在進(jìn)行濾波、重力補(bǔ)償?shù)幕A(chǔ)上,通過六維力傳感器測量打磨工具末端與蝸輪齒面的接觸力,完成對蝸輪齒面刀痕的柔順打磨,如圖10 所示.實(shí)驗(yàn)中設(shè)定恒定打磨力為20N,控制蝸輪齒面表面法向力為恒力,分別對兩個齒面進(jìn)行打磨,其中一個齒面利用沒有進(jìn)行力感知信息預(yù)處理的工業(yè)機(jī)器人系統(tǒng)進(jìn)行打磨,另一個齒面使用本文搭建的工業(yè)機(jī)器人打磨系統(tǒng)進(jìn)行打磨.打磨過程中法向受力如圖11 所示.相比于傳統(tǒng)工業(yè)機(jī)器人打磨系統(tǒng),本文搭建的打磨系統(tǒng),經(jīng)過力感知信息預(yù)處理后,恒力打磨時力的波動更小,更加穩(wěn)定,有效地減少了環(huán)境、誤差帶來的影響.

為進(jìn)一步驗(yàn)證本文力感知信息處理方法的效果,手持粗糙度測量儀對兩個齒面的表面粗糙度進(jìn)行測量,分別對不同區(qū)域進(jìn)行五次測量,且兩個齒面測量的區(qū)域相同,測量結(jié)果如表2 所示.可知在傳統(tǒng)打磨系統(tǒng),即無力感知信息預(yù)處理時,蝸輪齒面粗糙度為0.902 μm,而本文搭建的打磨系統(tǒng)的力感知信息預(yù)處理時,蝸輪齒面粗糙度為0.801 μm,打磨后蝸齒面輪表面質(zhì)量明顯改善.

處理方式無處理有處理測量一0.912 0.792測量二0.894 0.784測量三0.904 0.807測量四0.914 0.821測量五0.887 0.801均值0.902 0.801

5 小結(jié)

本文提出了一種蝸輪齒面刀痕打磨工業(yè)機(jī)器人接觸力獲取方法.針對工業(yè)機(jī)器人打磨過程中環(huán)境溫度、機(jī)器人振動等噪聲造成力感知信息波動的問題,建立了卡爾曼濾波數(shù)學(xué)模型;考慮到工業(yè)機(jī)器人底座安裝傾角、工具重力和六維力傳感器零點(diǎn)的影響,建立了六維力傳感器的零點(diǎn)校正和重力補(bǔ)償計算模型,完成了六維力傳感器的卡爾曼濾波和重力補(bǔ)償實(shí)驗(yàn)并進(jìn)行了數(shù)據(jù)分析.實(shí)驗(yàn)結(jié)果顯示,六維力傳感器數(shù)據(jù)的讀數(shù)誤差維持在±0.4 N 之間,力矩誤差維持在±0.16 N·m 之間,其讀數(shù)均在零點(diǎn)附近波動.與傳統(tǒng)的工業(yè)機(jī)器人打磨系統(tǒng)相比,經(jīng)過力感知信息預(yù)處理的工業(yè)機(jī)器人打磨系統(tǒng)更加穩(wěn)定、精確,并且加工出的齒面表面質(zhì)量明顯提高,為后續(xù)對工業(yè)機(jī)器人柔順打磨蝸輪齒面刀痕的研究奠定了基礎(chǔ).

猜你喜歡
蝸輪齒面卡爾曼濾波
新型材料的蝸輪蝸桿減速箱
基于NURBS理論的漸開線齒輪齒面修復(fù)
電梯蝸輪輪齒失效原因分析及警示
基于遞推更新卡爾曼濾波的磁偶極子目標(biāo)跟蹤
基于BP神經(jīng)網(wǎng)絡(luò)的面齒輪齒面粗糙度研究
高速動車組弧齒錐齒輪齒面疲勞點(diǎn)蝕失效分析
基于模糊卡爾曼濾波算法的動力電池SOC估計
17CrNiMo6齒輪軸齒面剝落原因分析
大型鑄鍛件(2015年1期)2016-01-12 06:33:01
基于擴(kuò)展卡爾曼濾波的PMSM無位置傳感器控制
基于Kisssoft與Workbench的汽車EPS蝸輪蝸桿優(yōu)化設(shè)計
大余县| 成都市| 德格县| 聊城市| 千阳县| 常德市| 宣城市| 苗栗市| 隆尧县| 清河县| 六枝特区| 突泉县| 定陶县| 屏山县| 宣武区| 伊金霍洛旗| 通海县| 赤水市| 正宁县| 宜君县| 南雄市| 灌南县| 沁水县| 读书| 丰宁| 长沙市| 宁德市| 鹤庆县| 景德镇市| 苏尼特右旗| 萨迦县| 商城县| 白水县| 湛江市| 弥勒县| 兴隆县| 潍坊市| 江北区| 延津县| 敦化市| 涿州市|